SUBROUTINE eta

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE PiFactrs_Module
USE Convrsns_Module
USE EnergyGrid_Module
USE Gaussian_Module
USE SurfaceJacobi_Module
USE APH_Module 
USE Masses_Module
IMPLICIT NONE
SAVE

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

LOGICAL :: firstpt
INTEGER :: ipt, istat
REAL(dp) :: k_in, gaussian, exparg, fac, realpart, imagpart
REAL(dp) :: etaemag(edim)
REAL(dp) :: etamag, normfac
COMPLEX(dp) :: etain(edim), etaout(edim)
!=========================================================================================

ALLOCATE(eta_e(edim))

firstpt=.true.

!=========================================================================================
! Obtain the initial gaussian k and energy distribution coefficients
! Use analytic expresion: See Mathematica notebook 'initial_E_extraction.nb'
IF (debug) THEN
OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//'iwf/initial_k_dist_in.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - eta'
OPEN(UNIT=output_unit1,FILE=TRIM(outdir)//'iwf/initial_e_dist.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - eta'
OPEN(UNIT=output_unit2,FILE=TRIM(outdir)//'iwf/initial_e_dist_eV.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - eta'
OPEN(UNIT=output_unit3,FILE=TRIM(outdir)//'iwf/initial_e_dist_kin.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - eta'
OPEN(UNIT=output_unit4,FILE=TRIM(outdir)//'iwf/initial_k_dist_out.txt',IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - eta'
ENDIF

! Incoming part 
fac=Sqrt(sigma)*(2.d0/pi)**0.25d0

eta_e=0.d0

DO ipt=1,edim

  IF (e_vals(ipt).ge.vib_energy_in) THEN

    IF(firstpt.eq..true.) initial_ept=ipt
    firstpt=.false.

    k_in=Sqrt(usys2*(e_vals(ipt)-vib_energy_in))

    exparg=(k_in-k0)*(k_in-k0)
    gaussian=fac*Exp(-exparg*sigma*sigma)
    realpart=Cos(s0*(k_in-k0))*gaussian
    imagpart=Sin(s0*(k_in-k0))*gaussian
!    realpart=Cos(s0*(k_in))*gaussian
!    imagpart=Sin(s0*(k_in))*gaussian
    etain(ipt)=Cmplx(realpart,imagpart,dp)
    etamag=Real(etain(ipt))**2+Aimag(etain(ipt))**2
    IF (debug) WRITE(output_unit0,*) k_in, etamag

!    exparg=(-k_in-k0)*(-k_in-k0)
!    gaussian=fac*Exp(-exparg*sigma*sigma)
!    realpart=Cos(s0*(-k_in-k0))*gaussian
!    imagpart=Sin(s0*(-k_in-k0))*gaussian
!    etaout(ipt)=Cmplx(realpart,imagpart,dp)
!    etamag=Real(etaout(ipt))**2+Aimag(etaout(ipt))**2
!    WRITE(output_unit4,*) k_in, etamag

    exparg=(k_in+k0)*(k_in+k0)
    gaussian=fac*Exp(-exparg*sigma*sigma)
    realpart=Cos(-s0*(k_in+k0))*gaussian
    imagpart=Sin(-s0*(k_in+k0))*gaussian
!    realpart=Cos(-s0*(k_in))*gaussian
!    imagpart=Sin(-s0*(k_in))*gaussian
    etaout(ipt)=Cmplx(realpart,imagpart,dp)
    etamag=Real(etaout(ipt))**2+Aimag(etaout(ipt))**2
    IF (debug) WRITE(output_unit4,*) k_in, etamag

    normfac=Sqrt(usys/Abs(k_in))
    IF (k0.gt.0.d0) eta_e(ipt)=normfac*etain(ipt) 
    IF (k0.lt.0.d0 )eta_e(ipt)=normfac*etaout(ipt)
    etaemag(ipt)=Sqrt(Real(eta_e(ipt))**2+Aimag(eta_e(ipt))**2)
    IF (debug) THEN
      WRITE(output_unit1,*) e_vals(ipt),etaemag(ipt)
      WRITE(output_unit2,*) e_vals(ipt)*autoev,etaemag(ipt)
      WRITE(output_unit3,*) e_vals(ipt)-vib_energy_in,etaemag(ipt)
    ENDIF
  ELSE

    !eta_e(ipt)=0.d0 ! could possibly cause errors
    !eta_cd(ipt,:)=0.d0
    etaemag(ipt)=Sqrt(Real(eta_e(ipt))**2+Aimag(eta_e(ipt))**2)
    IF (debug) THEN
      WRITE(output_unit1,*) e_vals(ipt),etaemag(ipt)
      WRITE(output_unit2,*) e_vals(ipt)*autoev,etaemag(ipt)
      WRITE(output_unit3,*) e_vals(ipt)-vib_energy_in,etaemag(ipt)
    ENDIF
  ENDIF

ENDDO

IF (debug) THEN
CLOSE(UNIT=output_unit0,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - eta'
CLOSE(UNIT=output_unit1,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - eta'
CLOSE(UNIT=output_unit2,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - eta'
CLOSE(UNIT=output_unit3,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - eta'
CLOSE(UNIT=output_unit4,IOSTAT=istat,STATUS='keep')
IF (istat.ne.0) STOP 'CLOSE failed - eta'
ENDIF

Return
End Subroutine eta