SUBROUTINE f_vector( ieng, f_vec )

!=========================================================================================
! Written by: Jeff Crawford
!
! This routine fourier transforms the time-dependent overlap integrals
!
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE PiFactrs_Module
USE EnergyGrid_Module
USE Jacobi_Module
USE SurfaceAPH_Module
USE Delves_Module
USE APH_Module
USE Overlap_Module
USE Print_Module
USE Time_Module
USE Masses_Module
USE QuantumNumber_Module
USE SymGroup_Module
USE Complex_Module
IMPLICIT NONE

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

INTEGER,INTENT(IN) :: ieng

!=========================================================================================
!   O U T P U T

COMPLEX(dp),INTENT(OUT) :: f_vec(totnaphsf)

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

INTEGER     :: tpt, tppt, istate, lambdaval, sfindex
REAL(dp)    :: intout_real, intout_imag
COMPLEX(dp) :: ft_int, coeff, num
!COMPLEX(dp) :: imag
!=========================================================================================
!   A L L O C A T A B L E

REAL(dp),ALLOCATABLE :: integrand_real(:),integrand_imag(:)

!=========================================================================================
!   F U N C T I O N S


!=========================================================================================
! Define i

! imag=Cmplx(0.d0,1.d0,dp)

!=========================================================================================
! Fourier transform all the collected Delves' theta "vibrational" function / propagated 
!    wave function overlap integrals at a specific energy

 sfindex=0

!ALLOCATE(integrand_real(tstep),integrand_imag(tstep))
 ALLOCATE(integrand_real(tstep-tdelaystep+1),integrand_imag(tstep-tdelaystep+1))

 ! Multiply by all factors here, such as 4/rho^5/2

 num = eta_e(ieng)*2.d0*pi*Sqrt(rho_infty**5)
 coeff = 4.d0/num

DO lambdaval=0,jtotal
 DO istate=1,naphsf(lambdaval)

!   DO tpt=1,tstep

    Do tpt=1,tstep-tdelaystep+1

        tppt=tpt+tdelaystep-1

        ft_int=overlap_vec(istate+sfindex,tppt)*Exp(imag*e_vals(ieng)*tppt*deltat)

        integrand_real(tpt)=Real(ft_int)
        integrand_imag(tpt)=AImag(ft_int) 

        ft_int=0.d0

   ENDDO


   CALL SimpExtInt(integrand_real,tstep-tdelaystep+1,deltat,intout_real)
   CALL SimpExtInt(integrand_imag,tstep-tdelaystep+1,deltat,intout_imag)

   f_vec(istate+sfindex)=coeff*Cmplx(intout_real,intout_imag,dp) 

 ENDDO

 sfindex=sfindex+naphsf(lambdaval)

ENDDO

DEALLOCATE(integrand_real,integrand_imag)

END SUBROUTINE f_vector