SUBROUTINE S_wavepacket(spacket)

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Convrsns_Module
USE PiFactrs_Module
USE Jacobi_Module
USE SurfaceJacobi_Module
USE Gaussian_Module
USE APH_Module
USE Masses_Module

IMPLICIT NONE
SAVE

!=========================================================================================
!     O U T P U T
COMPLEX(dp),INTENT(OUT) :: spacket(ndim)

!=========================================================================================
!     P A R A M E T E R S
LOGICAL,PARAMETER :: debug=.true.
!=========================================================================================
!     I N T E R N A L S

INTEGER  :: i, irho, itheta, ichi, ispt, istat, array(3)
REAL(dp) :: gcoord, gaussian, fac, exparg, realpart, imagpart, mag
REAL(dp) :: fws, fwk, fwe, roismin, roikmin, roiemin, roismax, roikmax, roiemax

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

!===============================================================================
! Calculate the gaussian wave packet as a function of capital S (Jacobi Mass-Scaled)
! Note that initial wavepacket is given a negative momentum.

IF (debug) THEN
  OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//"iwf/gaussian.txt",IOSTAT=istat,STATUS="replace",ACTION="write")
  IF (istat.ne.0) STOP 'S_wavepacket - OPEN failed'
ENDIF

  fac=1.d0/(2.d0*pi*sigma*sigma)**0.25d0

  DO irho=1,nrho 
    DO itheta=1,ntheta 
      DO ichi=1,nchi

        ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)

        gcoord=larges(ispt)
        exparg=(gcoord-s0)*(gcoord-s0)
        gaussian=fac*Exp(-exparg/(4.d0*sigma*sigma))
        realpart=Cos(-k0*gcoord)*gaussian
        imagpart=Sin(-k0*gcoord)*gaussian
        spacket(ispt)=Cmplx(realpart,imagpart,dp)
        mag=(Real(spacket(ispt)))**2+(Aimag(spacket(ispt)))**2
        IF (debug) write(output_unit0,*) real(gcoord,sp),real(mag,sp)

      ENDDO
    ENDDO 
  ENDDO 

  IF (debug) THEN
    CLOSE(UNIT=output_unit0,IOSTAT=istat,STATUS="KEEP")
    IF (istat.ne.0) STOP 'S_wavepacket - CLOSE failed'
  ENDIF

!=========================================================================================
! Output initial wave packet widths and regions of interest for Jacobi large S,
! wave number, and energy. 
! Put this in getinfo ???

 array=(/2,10,100/)

 WRITE(Out_Unit,31) 'Wave Packet Info:', 'Full Width', 'S', 'k, au', 'E, au', 'E, eV' 

 31 FORMAT(/1X,a,/1X,a,T21,a,T26,a,T34,a,T43,a,/1X,50('-'))
   
 DO i=1,3
   mag=Real(array(i),dp)
   fws=4.d0*sigma*Sqrt(Log(mag)) 
   fwk=2.d0*Sqrt(Log(mag))/sigma
   fwe=(fwk**2)/(usys2)
   WRITE(Out_Unit,32) array(i), fws, fwk, fwe, fwe*autoev
 ENDDO

 32 FORMAT(1X,I3,12('.'),T18,F5.2,T26,F5.2,T34,F6.3,T42,F5.2)  

  WRITE(Out_Unit,33) 'Wave Packet Info:', 'Region of Interest:', 'S', 'k, au', 'E, au','E, eV'

 33 FORMAT(/1X,a,/1X,a,T21,a,T26,a,T34,a,T43,a,/1X,50('-'))

 DO i=2,3
   mag=Real(array(i),dp)
   roismin=s0-2.d0*sigma*Sqrt(Log(mag))
   roismax=s0+2.d0*sigma*Sqrt(Log(mag))
   roikmin=k0-Sqrt(Log(mag))/sigma
   roikmax=k0+Sqrt(Log(mag))/sigma
   roiemin=(roikmin**2)/usys2+vib_energy_in
   IF (roikmin.lt.0.d0) roiemin=0.d0
   roiemax=(roikmax**2)/usys2+vib_energy_in
   WRITE(Out_Unit,34) 'min',array(i), roismin, roikmin, roiemin, roiemin*autoev
   WRITE(Out_Unit,34) 'max',array(i), roismax, roikmax, roiemax, roiemax*autoev
 ENDDO

 34 FORMAT(1X,a,1X,I3,8('.'),T18,F5.2,T26,F5.2,T34,F6.3,T42,F5.2)

Return
End Subroutine S_wavepacket