SUBROUTINE basis_functions
!
!=========================================================================================
! Written By: J. Crawford
!
! This routine obtains the asymptotic diatomic eigenfunctions as a function of
! Jacobi mass-scaled s. It also maps the functions to the APH grid and stores
! the functions in a binary file to be used later.
!
! Variables:
! jac_function = array storing Jacobi basis functions
! j_in = diatomic angular momentum j in initial wave packet
! jmax = max value of diatomic angular momentum j
! jtotal = value of total angular momentum
! l_in = space fixed orbital angular momentum of initial wave packet
! nvib_in = diatomic vibrational state in initial wave packet
! parity = parity of initial wave packet
! req_chanls = number of unique arrangement channels
! sdim = dimension of small s Jacobi grid (read in from namelist)
! s_dim = dimension of small s Jacobi grid (read in from binary)
! s_max = max small s Jacobi grid value
! s_min = min small s JAcobi grid value
! s_vals = array of small s Jacobi grid values
!=========================================================================================
Use Numeric_Kinds_Module
USE CommonInfo_Module
USE SurfaceJacobi_Module
USE SurfaceDelves_Module
USE SurfaceAph_Module
USE QuantumNumber_Module
USE Print_Module
USE SymGroup_Module
USE APH_Module
IMPLICIT NONE
!========================================================================================
! I N P U T
!=========================================================================================
! I N T E R N A L S
CHARACTER(LEN=1) :: c1
CHARACTER(LEN=2) :: j2
CHARACTER(LEN=100):: filename
INTEGER :: i, ichanl, ivib, ijrot, istat, jcheck, ccheck, neig, sdim, s_dim, tcdim
REAL(dp) :: s_min, s_max
!=========================================================================================
! A L L O C A T A B L E
REAL(dp),ALLOCATABLE :: jac_function(:,:), s_vals(:)
REAL(dp),ALLOCATABLE :: wfrot(:), wfvib(:)
!COMPLEX(dp),ALLOCATABLE :: wfrot(:), wfvib(:)
!=========================================================================================
! F U N C T I O N S
CHARACTER(LEN=1) :: label1
CHARACTER(LEN=2) :: label2
INTEGER :: aphindex, aphindex_tc
save
!=========================================================================================
! N A M E L I S T S
NAMELIST / surface_grid / sdim, td_dim
!=========================================================================================
! Read in NAMELIST parameters
OPEN(UNIT=Input_Unit,FILE=Trim(inputfile1),IOSTAT=istat,STATUS='old',ACTION='read')
IF (istat.ne.0) STOP 'OPEN failed - surface_function'
READ(Input_Unit,surface_grid)
CLOSE(UNIT=Input_Unit,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - surface_function'
!=========================================================================================
! Check input parameter validity
IF (j_in.gt.jmax) jmax=j_in
!=========================================================================================
! Output pertinent info to screen
WRITE(*,21) 'jtotal ',jtotal
WRITE(*,21) 'jmax ',jmax
WRITE(*,21) 'nvib_in',nvib_in
WRITE(*,21) 'j_in ',j_in
WRITE(*,21) 'l_in ',l_in
WRITE(*,21) 'parity ',parity
WRITE(*,22) 'Calculating Jacobi Basis Functions'
21 FORMAT(1X,a,8('.'),I3)
22 FORMAT(1X,a)
!=========================================================================================
! Determine number of unique arrangement channels
nchanls=3
req_chanls=nchanls
IF (system.eq.'AAA') req_chanls=1
IF (system.eq.'ABB') req_chanls=Min(2,nchanls)
!==========================================================================================
! Calculate Jacobi rovib functions
tcdim=ntheta*nchi
firstcallj = .true.
firstcalld = .true.
!=========================================================================================
! Obtain rovibrational eigenfunctions in Jacobi and Delves coordinates
ALLOCATE( wfrot(tcdim), wfvib(tcdim) )
DO ichanl=1,nchanls ! DO (1) : channel loop
c1=label1(ichanl)
CALL jacobi_grid(ichanl)
DO ijrot = 0, jmax ! DO (2) : j loop
j2=label2(ijrot)
IF (ichanl.le.req_chanls) THEN ! IF (1) : check for req_chanls
IF (.not.use_pvs.or.ijrot.eq.j_in) THEN ! IF (2) : check for use_pvs or j_in
!==================================================================================
! Obtain vibrational basis functions for j=ijrot in arrangement channel=ichanl
CALL jacobi_basis( ijrot, ichanl, sdim )
!CALL delves_surface( ijrot, ichanl)
ENDIF ! IF (2) : use_pvs or j_in
ENDIF ! IF (1) : req_chanls
!===================================================================================
! Read in rovibrational small s rovib functions calculated in jacobi_surface routine
IF (ichanl.gt.req_chanls) c1=label1(req_chanls)
filename='jacobi_surface_c'//c1//'_'//j2//'.bin'
OPEN(UNIT=bin_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
READ(bin_unit1) ccheck
READ(bin_unit1) jcheck
IF (jcheck.ne.ijrot) STOP 'wrong jjval readin : surface_functions'
READ(bin_unit1) neig
READ(bin_unit1) s_dim
!IF (sdim_read.ne.sdim) STOP 'wrong sdim readin : surface_functions'
IF(ALLOCATED(s_vals)) DEALLOCATE(s_vals)
ALLOCATE( s_vals(s_dim))
READ(bin_unit1) s_vals
READ(bin_unit1) s_min
READ(bin_unit1) s_max
IF(ALLOCATED(jac_function)) DEALLOCATE(jac_function)
ALLOCATE( jac_function(s_dim,neig))
READ(bin_unit1) jac_function
! IF (ichanl.eq.nchanls) THEN
! IF (ichanl.ne.1.and.ijrot.ne.j_in) THEN
! CLOSE(UNIT=bin_unit1, STATUS='delete')
! ELSE
! CLOSE(bin_unit1)
! ENDIF
! ELSE
CLOSE(bin_unit1)
! ENDIF
IF (.not.use_pvs) THEN ! IF (3) : check for use_pvs
!===========================================================================================
! Open output file to store Jacobi small s cap theta rovib functions until needed in analysis
c1=label1(ichanl)
filename='jacobi_totsurface_c'//c1//'_'//j2//'.bin'
OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
!============================================================================
! Obtain rotational component of asymptotic basis function mapped to APH grid
IF (ntheta.eq.1) THEN
wfrot=0.d0
!wfrot=Cmplx(1.d0,0.d0,dp)
ELSE
CALL rotational(ijrot,0,tcdim, wfrot,'O','SF') ! Obtain cap theta legendre polynomial
ENDIF
DO ivib=1,neig ! DO (3) : vibrational loop
!======================================================
! Obtain vibrational basis function mapped to APH grid
CALL vibrational(s_dim, s_vals, s_min, s_max, jac_function(:,ivib),tcdim,wfvib,'O')
wfvib=wfvib*wfrot
WRITE(bin_unit0) wfvib
ENDDO ! DO (3) : vib loop
CLOSE(bin_unit0)
ENDIF ! IF (3) use_pvs
ENDDO ! DO (2) : j loop
firstcallj=.true.
ENDDO ! DO (1) : channel loop
!=========================================================================================
! Deallocate arrays
DEALLOCATE(wfrot,wfvib)
END SUBROUTINE basis_functions