SUBROUTINE zanalysis(whichir)

!=========================================================================================
! Written by: Jeff Crawford
!
! This routine performs the state-to-state time dependent analysis using Delves' 
!    coordinates. The propagated wavefunction is projected onto a surface of constant rho.
!    Then, overlaps between this projection and Delves' theta eigenfunctions are 
!    calculated and Fourier transformed (F vector). A coefficent matrix is generated that
!    contains the overlaps between the Delves' theat eigenfunctions and the mass-scaled
!    Jacobi eigenfunctions in S and s (A matrix).
! These are used to solve the matrix equation AS=F for the S matrix column S.

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Jacobi_Module
USE APH_Module
USE Delves_Module
USE Overlap_Module
USE EnergyGrid_Module
USE QuantumNumber_Module
USE SurfaceJacobi_Module
USE SurfaceAPH_Module
USE SymGroup_Module
USE Masses_Module
use SMat_Module
USE Complex_Module
IMPLICIT NONE
SAVE

!=========================================================================================
!   I N P U T

INTEGER,INTENT(IN) :: whichir

!=========================================================================================
!   I N T E R N A L S

CHARACTER(LEN=1)   :: ir1
CHARACTER(LEN=100) :: filename

LOGICAL :: anyopen

INTEGER     :: i, j, ipt, jstart,  jbfdim, jsfdim, info, ieng, jeng
REAL(dp)    :: eval
!COMPLEX(dp) :: imag

!=========================================================================================
!   A L L O C A T A B L E

INTEGER,ALLOCATABLE     :: ipiv(:), njvib_open(:,:), njvib_max(:,:),njbf(:)
INTEGER,ALLOCATABLE     :: maxopen(:), j_vector(:)
COMPLEX(dp),ALLOCATABLE :: f_vec(:), a_mat(:,:), ata(:,:), atf(:)

!=========================================================================================
!   F U N C T I O N S

CHARACTER(LEN=1) :: label1

!=========================================================================================
!   N A M E L I S T S

!=========================================================================================

 print*,'----------------------------------------------------------------------------'
 print*,'Starting final analysis'
 print*,'----------------------------------------------------------------------------'

 ir1=label1(whichir)

! imag=Cmplx(0.d0,1.d0,dp)

 ALLOCATE(njvib_open(0:jmax,req_chanls))
 ALLOCATE(njbf(nchanls))
 ALLOCATE(maxopen(req_chanls))

 firstcall=.true.
 njvib_open=0
 maxopen=0
!========================================================================================
! Open binary file to store the s matrices at each energy
!---------------------------------------------------------------------------------------- 
   IF (whichir.le.4) THEN
     filename='smatrix_'//ir1//'.bin'
     OPEN(UNIT=smat_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
   ENDIF
!========================================================================================
! Open file to store the irreducible representation s matrices at each energy
!----------------------------------------------------------------------------------------
!    filename='smatrix_ir'//ir1//'.bin'
!    OPEN(UNIT=smatrix_unit19,FILE=TRIM(outdir)//TRIM(filename))
!=========================================================================================
! Begin smatrix calculation at each energy 
!-----------------------------------------------------------------------------------------
 DO ieng=initial_ept,edim ! Energy loop

   !=====================================================================================
   ! Determine the number of open jacobi  states
   !-------------------------------------------------------------------------------------

   CALL open_states( ieng, njvib_open, maxopen, anyopen, whichir)

   !=====================================================================================
   ! Obtain symmetrized jacobi eigenfunctions
   !-------------------------------------------------------------------------------------

   IF (anyopen) THEN

     CALL zget_sym_jacobi(ieng, njvib_open, njbf,  whichir)

     IF (nirrep.eq.8) THEN ! C6v
       jsfdim = njsf_ir
       jbfdim = njbf(2)
     ELSEIF (nirrep.eq.4) THEN ! C2v
       jsfdim = njsf_a+njsf_ir
       jbfdim = njbf(1)+njbf(2)
     ELSE
       jsfdim = njbf(1)+njbf(2)+njbf(3) !C2
       jbfdim = njbf(1)+njbf(2)+njbf(3)
     ENDIF

     ALLOCATE(j_vector(jbfdim))

     CALL combine_abc(jbfdim, j_vector, njbf, whichir)

     ALLOCATE( a_mat(totnaphsf,jsfdim) )
     CALL za_matrix(ieng, jbfdim, a_mat, whichir)

     ALLOCATE( f_vec(totnaphsf) )
     CALL f_vector( ieng, f_vec )

     ALLOCATE( ata(jbfdim,jbfdim), atf(jbfdim), ipiv(jbfdim) )
     ata=Matmul(Transpose(a_mat),a_mat)
     atf=Matmul(Transpose(a_mat),f_vec) 


    CALL zgesv( jbfdim , 1 , ata , jbfdim , ipiv , atf , jbfdim , info)
    !zgesv(N=dim,NRHS=col of F,A=coef matrix, LDA=dim of A , IPIV = interchange indicator , B=RHS matrix , LDB , INFO)

    IF (info.ne.0) THEN
      PRINT*,'Error : U factorization element i is zero : i = ',info
      PRINT*,       'for energy ieng = ',ieng
      PRINT*,' '
      STOP 'analysis'
    ENDIF

    !=====================================================================================
    ! Output smatrix info
    !-------------------------------------------------------------------------------------

    ! Output irreducible representation smatrix
    !eval=e_vals(ieng)*autoev
    !WRITE(smatrix_unit19,'(1X,e14.7,1000(1X,e14.7))') eval, (atf(ipt),ipt=1,jsfdim)   

    IF (nirrep.eq.8) THEN ! IF C6v
      IF (whichir.le.4) THEN ! IF A1,A2,B1,B2

        write(smat_unit1) jbfdim
        write(smat_unit1) j_vector
        write(smat_unit1) atf

        IF (ieng.eq.edim) THEN
          react_dim=maxopen(req_chanls)
          nonreact_dim=njsf_ir
          ALLOCATE(njvib_max(0:jmax,req_chanls))
          njvib_max=njvib_open

          OPEN(unit=Output_Unit0,file=TRIM(outdir)//"jr_vector.txt")
          OPEN(unit=Output_Unit1,file=TRIM(outdir)//"jnr_vector.txt")
          DO ipt=0,jmax
            IF (njvib_max(ipt,1).gt.0) THEN
              WRITE(Output_Unit0,'(i3)') njvib_max(ipt,1)
              IF (mod(j_in,2).eq.mod(ipt,2)) WRITE(Output_Unit1,'(i3)') njvib_max(ipt,1)
            ENDIF
          ENDDO
          CLOSE(Output_Unit0)
          CLOSE(Output_Unit1)

          jstart=Mod(j_in,2)

          WRITE(Out_Unit,*) 'List of Open States: ( j ) highest open vib for j '
          WRITE(Out_Unit,*) 'NONREACTIVE'
          WRITE(Out_Unit,'(1000(a,i3,a,i3))')('(',ipt,')',njvib_max(ipt,1),ipt=jstart,jmax,2)

        ENDIF

      ELSE ! IF E reps

        CALL zsmatrix(ieng, jbfdim, j_vector, atf, njvib_max)

        IF (ieng.eq.edim) THEN
          WRITE(Out_Unit,*) 'REACTIVE'
          WRITE(Out_Unit,'(1000(a,i3,a,i3))')('(',ipt,')',njvib_max(ipt,1),ipt=0,jmax)
        ENDIF
      ENDIF

    ELSEIF (nirrep.eq.4) THEN ! IF C2v
  
      write(smat_unit1) njsf_a
      write(smat_unit1) njsf_ir
      write(smat_unit1) jbfdim
      write(smat_unit1) j_vector
      write(smat_unit1) atf   

      IF (ieng.eq.edim) THEN
        CLOSE(smat_unit1)
        react_dim=maxopen(req_chanls)
        nonreact_dim=njsf_a
        ALLOCATE(njvib_max(0:jmax,req_chanls))
        njvib_max=njvib_open

        OPEN(unit=Output_Unit0,file=TRIM(outdir)//"jr_vector.txt")
        OPEN(unit=Output_Unit1,file=TRIM(outdir)//"jnr_vector.txt")
        DO ipt=0,jmax
          IF (njvib_max(ipt,2).gt.0) WRITE(Output_Unit0,'(i3)') njvib_max(ipt,2)
          IF (njvib_max(ipt,1).gt.0.and.mod(j_in,2).eq.mod(ipt,2)) THEN 
            WRITE(Output_Unit1,'(i3)') njvib_max(ipt,1)
          ENDIF
        ENDDO
        CLOSE(Output_Unit0)
        CLOSE(Output_Unit1)

        jstart=Mod(j_in,2)

        write(Out_Unit,*) 'List of Open States: ( j ) highest open vib for j '
        write(Out_Unit,*) ' NONREACTIVE  '
        write(Out_Unit,'(1000(a,i3,a,i3))') ('(',ipt,')',njvib_max(ipt,1), ipt=jstart,jmax,2)
        write(Out_Unit,*) ' REACTIVE  '
        write(Out_Unit,'(1000(a,i3,a,i3))') ('(',ipt,')',njvib_max(ipt,2), ipt=0,jmax)

        DO jeng=initial_ept,edim
          CALL szsmatrix(jeng, njvib_max,  whichir)
        ENDDO

      ENDIF
   
    ELSE ! C2

      write(smat_unit1) njbf(1)
      write(smat_unit1) njbf(2)
      write(smat_unit1) njbf(3)
      write(smat_unit1) jbfdim
      write(smat_unit1) j_vector
      write(smat_unit1) atf

      IF (ieng.eq.edim) THEN
        CLOSE(smat_unit1)
        ALLOCATE(njvib_max(0:jmax,req_chanls))
        njvib_max=njvib_open

        OPEN(unit=Output_Unit0,file=TRIM(outdir)//"ja_vector.txt")
        OPEN(unit=Output_Unit1,file=TRIM(outdir)//"jb_vector.txt")
        OPEN(unit=Output_Unit2,file=TRIM(outdir)//"jc_vector.txt")
        DO ipt=0,jmax
          IF (njvib_max(ipt,1).gt.0) WRITE(Output_Unit0,'(i3)') njvib_max(ipt,1)
          IF (njvib_max(ipt,2).gt.0) WRITE(Output_Unit1,'(i3)') njvib_max(ipt,2)
          IF (njvib_max(ipt,3).gt.0) WRITE(Output_Unit2,'(i3)') njvib_max(ipt,3)
        ENDDO
        CLOSE(Output_Unit0)
        CLOSE(Output_Unit1)
        CLOSE(Output_Unit2)

        write(Out_Unit,*) 'List of Open States: ( j ) highest open vib for j '
        write(Out_Unit,*) ' CHANNEL A  '
        write(Out_Unit,'(1000(a,i3,a,i3))') ( '(',ipt,')',njvib_max(ipt,1),ipt=0,jmax )
        write(Out_Unit,*) ' CHANNEL B  '
        write(Out_Unit,'(1000(a,i3,a,i3))') ( '(',ipt,')',njvib_max(ipt,2),ipt=0,jmax )
        write(Out_Unit,*) ' CHANNEL C  '
        write(Out_Unit,'(1000(a,i3,a,i3))') ( '(',ipt,')',njvib_max(ipt,3),ipt=0,jmax )

        DO jeng=initial_ept,edim
          CALL smatrixABC(jeng, njvib_max, maxopen, whichir)
        ENDDO

      ENDIF

    ENDIF

    DEALLOCATE( a_mat, f_vec, ata, atf, ipiv, j_vector )

  ENDIF ! Open states loop

ENDDO ! Energy loop

CLOSE(smat_unit1)
CLOSE(smatrix_unit19)
DEALLOCATE(njvib_open, maxopen, njbf)

END SUBROUTINE zanalysis