SUBROUTINE aph_surface(lambdaval, whichir)
!
!========================================================================================
! Written by: Jeff Crawford  : adapted from OneDim codes written by G. Parker
!
! This routine obtains the asymptotic aph eigenfunctions at a fixed value of rho.
! The eigenfunctions are solutions to equation 164 in the APH paper (JCP 87, 3888 (1987))
!
Use Numeric_Kinds_Module
USE CommonInfo_Module
USE APH_Module
USE SurfaceAph_Module
USE SurfaceJacobi_Module
USE QuantumNumber_Module
!USE Hamilt_Module
USE Aph_Ham_module
USE Masses_Module
USE SymGroup_Module
USE Convrsns_Module
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: lambdaval , whichir

!=========================================================================================
!   P A R A M E T E R S
LOGICAL,PARAMETER :: debug=.true.
CHARACTER(LEN=1),PARAMETER:: Jobz='V', UpLo='L', Range='V'
INTEGER, PARAMETER:: ndaf=2, mval=32, ldab=42
REAL(dp),PARAMETER:: AbsTol=1.d-12, tol=1.d-12

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

CHARACTER(LEN=2)   :: j2, l2, label2
CHARACTER(LEN=100) :: filename

INTEGER  :: ipt, indexchi, i, j, itheta, jtheta, ichi, jchi, sec_pt 
INTEGER  :: neig, lwork, liwork, info, nu, istat, nchir, ndim_sec, chishift 
REAL(dp) :: e, mag, x, y, upot, sum_chi, sum_theta, theta, chi, rhoinf2, factor, v0re,v0rc

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

INTEGER,ALLOCATABLE  :: iwork(:), ifail(:)
REAL(dp),ALLOCATABLE :: psi(:,:), psisq(:), eig(:), work(:), psi_temp(:,:), eig_temp(:)
REAL(dp),ALLOCATABLE :: kinetic(:,:), theta_ham(:,:), chi_ham(:,:), pot_ham(:,:),hamil(:,:)

!=========================================================================================
!   F U N C T I O N S
INTEGER :: aphindex_tc

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

 j2=label2(jtotal)
 l2=label2(lambdaval)

 indexchi = indexchi_array(whichir)
 nchir=nchi_ir(whichir)
 ndim_sec=ntheta*nchir
 chishift=chishift_ir(whichir)

!=========================================================================================
! Allocate arrays and initialize variables

 rhoinf2=rho_infty*rho_infty

 ALLOCATE( eig(ndim_sec), psi(ndim_sec,ndim_sec), psisq(ndim_sec) )
 ALLOCATE( ifail(ndim_sec) )
 ALLOCATE( kinetic(ndim_sec,ndim_sec) )

 liwork=MAX(1,5*ndim_sec)
 lwork=8*ndim_sec

 ALLOCATE( work(lwork), iwork(liwork) )

 psi=0.d0
 eig=0.d0
 ifail=0
 info=0
 neig=0
 mag=0.d0
 work=0.d0
 iwork=0

!=========================================================================================
! Create coordinate and diatomic potential grids
! Uses two-body potential

  IF (lambdaval.eq.0) THEN

    !ALLOCATE( aphpot(ndim_sec) )

    ALLOCATE( apbover2(ntheta) , lterm(ntheta) )
    IF (.not.ALLOCATED(naphsf)) ALLOCATE(naphsf(0:jtotal))

   ALLOCATE( pot_ham(ndim_sec,ndim_sec) )

   pot_ham=0.d0

    ! DEBUG
    ! Output potential data

    OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//'potential/aph_v_new.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
    IF (istat.ne.0) STOP 'OPEN failed - surface_aph'
    DO itheta=1,ntheta
      DO ichi=1,nchir
        i=aphindex_tc(nchir,itheta,ichi)
        theta=theta_val(itheta)
        chi=chi_val(ichi+chishift)
!        upot = potential_full(ind_rho_infty,itheta,ichi+chishift)+(15.d0/(8.d0*usys*rhoinf2))
        upot=v_pot(ind_rho_infty,itheta,ichi+chishift)+(15.d0/(8.d0*usys*rhoinf2))
        pot_ham(i,i)=upot
        x=TAN(0.5d0*theta)*COS(chi)
        y=TAN(0.5d0*theta)*SIN(chi)
        Write(output_unit0,*) real(x,sp),real(y,sp),real(upot,sp)
      ENDDO
    ENDDO
    CLOSE(output_unit0,IOSTAT=istat,STATUS='keep')
    IF (istat.ne.0) STOP 'CLOSE failed - surface_aph'
    ! Rotational terms
    factor=rhoinf2*usys
    IF (jtotal.eq.0) THEN
      apbover2=0.d0
      lterm=0.d0
    ELSE
      DO itheta=1,ntheta
        apbover2(itheta) = (a(itheta)+b(itheta))/(2.d0*factor)
        lterm(itheta) = (c(itheta)/factor) - apbover2(itheta)
      ENDDO
    ENDIF
 ENDIF

!=========================================================================================
! Obtain full 2D Hamiltonian: 

  ALLOCATE( theta_ham(ndim_sec,ndim_sec) , chi_ham(ndim_sec,ndim_sec) , hamil(ndim_sec,ndim_sec) ) 

  theta_ham = 0.d0
  chi_ham = 0.d0
  hamil = 0.d0

  ! Isolate correct portion of chi kinetic energy operator
  DO itheta=1,ntheta
    DO jchi=1,nchir
      DO ichi=1,nchir
        i = aphindex_tc(nchir, itheta, ichi)
        j = aphindex_tc(nchir, itheta, jchi)
        chi_ham(i,j)=2.d0*b(itheta)*t_chi(ichi+indexchi,jchi+indexchi)/rhoinf2
        !chi_ham(i,j)=2.d0*b(itheta)*kinetic_chi(ichi+indexchi,jchi+indexchi)/rhoinf2
      ENDDO
    ENDDO
  ENDDO

      ! Theta
      DO itheta = 1, ntheta
         DO jtheta = 1, ntheta
            DO ichi = 1, nchir
               i = aphindex_tc(nchir, itheta, ichi)
               j = aphindex_tc(nchir, jtheta, ichi)
               theta_ham(i,j) = t_theta (itheta, jtheta) / rhoinf2
               !theta_ham(i,j) = kinetic_theta (itheta, jtheta) / rhoinf2
            ENDDO
         ENDDO
      ENDDO

hamil = theta_ham + chi_ham + pot_ham


v0re=0.d0

v0rc=vasyj!+(15.d0/(8.d0*usys*rhoinf2))
!v0rc=0.05d0

PRINT*,'v0rc = ',v0rc,' au,   ',v0rc*autoev,' eV '

!=========================================================================================
!! Diagonalize with dsyevx

  CALL dsyevx(Jobz, Range, UpLo, ndim_sec, hamil, ndim_sec, v0re, v0rc, 1, 1, AbsTol, neig, eig, psi, ndim_sec, work, lwork, iwork, ifail, info)

!=========================================================================================
! Error message check for dsyevx

  DO i=1,neig
    IF (ifail(i).ne.0.or.info.ne.0) THEN 
      PRINT*,'ERROR in routine dsyevx'  
      PRINT*,'ifail =',ifail(i),'  info=',info   
      STOP 'aph_surface'
    ENDIF
  ENDDO

  DEALLOCATE( ifail, iwork, work )

 print*,'neig =',neig

!=========================================================================================
! Divide by sqrt(delta_td) to provide proper amplitude (for proper normalization)

  psi=psi/Sqrt(deltachi)

!=========================================================================================
! Save number of bound states for current lambda value

 naphsf(lambdaval)=neig

 IF (lambdaval.eq.0) THEN
   totnaphsf=neig
 ELSE
   totnaphsf=totnaphsf+neig
 ENDIF

!=========================================================================================
! Output eigenfunctions and energies to file:

! Energies
 IF (debug) THEN
   OPEN(UNIT=output_unit2,FILE=TRIM(outdir)//'surface_functions/aph/energy_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_aph'
   OPEN(UNIT=output_unit3,FILE=TRIM(outdir)//'surface_functions/aph/energy_long_'//l2//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
   IF (istat.ne.0) STOP 'OPEN failed - surface_aph'
   DO nu=1,neig
     e=eig(nu)
     WRITE(output_unit2,'(i4,e14.7)') nu, e
     WRITE(output_unit3,'(a,i4,a,3(es23.15,a))')   "nu=", nu, '  Energy=',e, "(Ha)  ", e*27.2117d0,"(eV)  ",e*219474.6314d0,"(cm-1)"
   ENDDO
   CLOSE(UNIT=output_unit2,IOSTAT=istat,STATUS='keep')
   IF (istat.ne.0) STOP 'CLOSE failed - surface_aph'
   CLOSE(UNIT=output_unit3,IOSTAT=istat,STATUS='keep')
   IF (istat.ne.0) STOP 'CLOSE failed - surface_aph'
 ENDIF
! Wave Functions

!DO i=1,neig
!     l3=label3(i)
!     OPEN(UNIT=output_unit4,FILE=TRIM(outdir)//'surface_functions/aph/function_'//l3//'.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
!       IF (istat.ne.0) STOP 'OPEN failed - surface_aph'
!     DO itheta=1,ntheta
!       DO ichi=1,nchir
!         ipt=aphindex_tc(nchir, itheta, ichi)
!         x=TAN(0.5d0*theta_val(itheta))*COS(chi_val(ichi+chishift))
!         y=TAN(0.5d0*theta_val(itheta))*SIN(chi_val(ichi+chishift))
!         IF (ABS(psi(ipt,i)).lt.1.d-12) THEN
!           WRITE(output_unit4,*) real(x,sp) , real(y,sp) , real(0.d0,sp)
!         ELSE
!           WRITE(output_unit4,*) real(x,sp) , real(y,sp) , real(psi(ipt,i),sp)
!         ENDIF
!       ENDDO
!     ENDDO
!     CLOSE(output_unit4,IOSTAT=istat,STATUS='keep')
!     IF (istat.ne.0) STOP 'CLOSE failed - surface_aph'
!ENDDO

!=========================================================================================
! Create binaries to be read in during analysis

  ALLOCATE( eig_temp(neig) , psi_temp(ndim_sec,neig) )
  DO ipt=1,neig
    eig_temp(ipt)=eig(ipt)
    psi_temp(:,ipt)=psi(:,ipt)
  ENDDO

   filename='aph_surface_'//l2//'.bin'
   OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
   WRITE(bin_unit0) jtotal ! Use this to make sure the correct data is read in later
   WRITE(bin_unit0) neig
   WRITE(bin_unit0) ndim_sec
   WRITE(bin_unit0) eig_temp
   WRITE(bin_unit0) psi_temp
   CLOSE(bin_unit0,STATUS='keep')

  DEALLOCATE( eig_temp , psi_temp )

!=========================================================================================
!  Check normalization of eigenfunctions

  vibloop: DO i=1,neig
    sum_theta=0.d0
    DO itheta=1,ntheta
      sum_chi=0.d0
      DO ichi=1,nchir
        sec_pt=aphindex_tc(nchir,itheta,ichi)
        sum_chi = sum_chi + psi(sec_pt,i)*psi(sec_pt,i)*deltachi
      ENDDO
      sum_theta=sum_theta+sum_chi!*weights(itheta)
    ENDDO

    IF (ABS(sum_theta-1.d0).gt.tol) THEN ! Check that it is normalized
      PRINT*,'NORMALIZATION ERROR: APH SURFACE FUNCTION:' 
      PRINT*,'EIG = ',i,'     1.0-NORM = ',1.d0-mag
    ENDIF
  ENDDO vibloop

!=========================================================================================
! Deallocate Arrays

  DEALLOCATE( psi, psisq, eig)
  DEALLOCATE(kinetic)  

  IF (lambdaval.eq.jtotal) THEN
    DEALLOCATE(apbover2,lterm)
    DEALLOCATE(pot_ham)
  ENDIF
 
END SUBROUTINE aph_surface