SUBROUTINE analysis(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     :: ipt, 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(:,:)
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(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)

   !IF(e_vals(ieng).lt.vib_energy_in) anyopen=.false.

   !: write these to a binary because they won't change for each irred. rep.(can just read them in)

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

   !print*,'maxopen   = ',maxopen

   IF (anyopen) THEN

     CALL zget_sym_jacobi(ieng, njvib_open,  whichir)

     IF (nirrep.eq.8) THEN
       jsfdim = njsf_ir
     ELSEIF (nirrep.eq.4) THEN
       jsfdim = njsf_a+njsf_ir
     ELSE
       jsfdim = totnjsf
     ENDIF

     ALLOCATE(j_vector(jsfdim))

     CALL combine_abc(jsfdim,j_vector,whichir)

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

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

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


    CALL zgesv( jsfdim , 1 , ata , jsfdim , ipiv , atf , jsfdim , 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) jsfdim
        write(smat_unit1) j_vector
        write(smat_unit1) atf

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

! Fix this to output as c2v does
          OPEN(unit=output_unit0,file=TRIM(outdir)//"j_vector_"//ir1//".txt")
            write(output_unit0,'(1000(a,i3,a,i3))') ('(',ipt,')',j_vector(ipt), ipt=1,nonreact_dim)
          CLOSE(output_unit0)

        ENDIF

      ELSE ! IF E reps

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

      ENDIF
! fix this too
      IF (ieng.eq.edim) THEN
        OPEN(unit=output_unit0,file=TRIM(outdir)//"j_vector_"//ir1//".txt")
            write(output_unit0,'(1000(a,i3,a,i3))') ('(',ipt,')',j_vector(ipt), ipt=1,react_dim)
          CLOSE(output_unit0)
      ENDIF

    ELSE ! IF C2v or C2
  
      write(smat_unit1) njsf_a
      write(smat_unit1) njsf_ir
      write(smat_unit1) jsfdim
      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
          print*,'nonreact_dim =',nonreact_dim
          print*,'react_dim =',react_dim
          print*,'njsf_ir =',njsf_ir
          ALLOCATE(njvib_max(0:jmax,req_chanls))
          njvib_max=njvib_open

         !OPEN(unit=output_unit0,file=TRIM(outdir)//"open_states.txt")
            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=0,jmax)
            write(Out_Unit,*) ' REACTIVE  '
            write(Out_Unit,'(1000(a,i3,a,i3))') ('(',ipt,')',njvib_max(ipt,2), ipt=0,jmax)
          !CLOSE(output_unit0)


          DO jeng=initial_ept,edim
            CALL szsmatrix(jeng, njvib_max,  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)

END SUBROUTINE analysis