SUBROUTINE project_aph(tterm, phi_ir, ndim_ir, whichir,nchir,ndim_sector)

!==============================================================================
!   Written by: Jeff Crawford
!
! This routine obtains the vector of final state / propagated wave function overlaps at 
!    each time step
! These overlaps will be fourier transformed later to find the S matrix by solving the 
!    matrix equation AS=F.

USE Numeric_Kinds_Module
USE CommonInfo_Module
USE EnergyGrid_Module
USE APH_Module
USE SymGroup_Module
USE QuantumNumber_Module
USE SurfaceAPH_Module
USE Time_Module
USE Overlap_Module

IMPLICIT NONE

!==============================================================================
!   I N P U T / O U T P U T

INTEGER,INTENT(IN) :: tterm, ndim_ir, whichir, nchir, ndim_sector
COMPLEX(dp) :: phi_ir(ndim_ir)

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

CHARACTER(LEN=2)   :: l2
CHARACTER(LEN=100) :: filename
INTEGER     :: isf, lambdaval, sfindex, jcheck, neig, ndim_sec
REAL(dp)    :: probamp
COMPLEX(dp) :: overlap_element 

!==============================================================================
!   A L L O C A T A B L E
REAL(dp), ALLOCATABLE :: sfunction(:,:), senergy(:)
!COMPLEX(dp),ALLOCATABLE :: sector_wf(:)

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

CHARACTER(LEN=2) label2

!==============================================================================
! Obtain propagated sector wave function at rho_infty

  !nchir=nchi_ir(whichir)
  !ndim_sector=ntheta*nchir
  IF (tterm.eq.1) THEN
    ALLOCATE(sector_wf(ndim_sector))
  ENDIF

  CALL get_sector( tterm, phi_ir, ndim_ir, sector_wf, ndim_sector, whichir, nchir )

 sfindex=0

 DO lambdaval=0,jtotal

   l2=label2(lambdaval)
   ! Read in surface functions and energies
   filename='aph_surface_'//l2//'.bin'
   OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted')
     READ(bin_unit0) jcheck
     IF (jcheck.ne.jtotal) STOP 'wrong jtotal surface functions in project_aph'
     READ(bin_unit0) neig
     READ(bin_unit0) ndim_sec
     IF (ndim_sec.ne.ndim_sector) STOP 'wrong ndim_sec for surface functions in project_aph'
     ALLOCATE( sfunction(ndim_sector,neig), senergy(neig))
     READ(bin_unit0) senergy
     READ(bin_unit0) sfunction
   CLOSE(bin_unit0)

   DO isf=1,naphsf(lambdaval)

     CALL overlap_aph(sector_wf, sfunction(:,isf), ndim_sector, overlap_element, whichir, nchir)

     overlap_vec(isf+sfindex,tterm) = overlap_element

   ENDDO

   sfindex=sfindex+naphsf(lambdaval)

 ENDDO

 IF (tterm.eq.tstep) THEN
   DEALLOCATE(sector_wf)
 ENDIF

END SUBROUTINE project_aph