Subroutine tcheby(irdim, phi_in, whichir, nchir)

Use Numeric_Kinds_Module
Use CommonInfo_Module
Use User_Module
Use APH_Module
Use NIP_Module
Use Time_Module
Use Cheby_Module
Use SymGroup_Module
Use Complex_Module
Use Phi_Module
Implicit None
Save

!=========================================================================================
! Tchebychev recurence relation.
!    packet is the initial wavefunction
!    T_i is a function of hnorm ==> T_i[hnorm]
! T_n[hnorm]|packet>=2*H*T_(n-1)[hnorm]|packet>-T_(n-2)[hnorm]|packet> 

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

Integer,Intent(In) :: irdim, whichir, nchir

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

Integer     :: ir_pt, irho, itheta, ichi, iterm, aphindex
Complex(dp) :: phi_in(irdim)
!Complex(dp) :: imag
!=========================================================================================
!     A L L O C A T A B L E

!REAL(dp),ALLOCATABLE    :: phi_real(:), phi_imag(:)
!COMPLEX(dp),ALLOCATABLE :: phi_n(:), phi_nm1(:), phi_nm2(:)

!=========================================================================================
! ALLOCATE( phi_real(irdim), phi_imag(irdim) )
! ALLOCATE( phi_n(irdim), phi_nm1(irdim), phi_nm2(irdim) )

!=========================================================================================
!  Initialize Chebychev summation
!     chi_1 = |chi> = chi_nm2      initial wave function
!     chi_2 = hnorm|chi> = chi_nm1

!CALL findtime(.true.,'init')

  phi_nm2=phi_in
!  phi_nm1=phi_in
!  phi_real=real(phi_nm1)
!  phi_imag=aimag(phi_nm1)
  phi_real=real(phi_in,dp)
  phi_imag=aimag(phi_in)

  Call h_apply( phi_real, irdim, whichir )
  Call h_apply( phi_imag, irdim, whichir )

  phi_nm1=Cmplx(phi_real,phi_imag,dp)
  phi_nm1=-imag*phi_nm1

!=========================================================================================
!  Accumulate first two terms into sum.

!   phi_in=0.d0  ! added to reset the input vector
   phi_in=phi_nm2*bessel_func(1)
   phi_in=phi_in+(phi_nm1*2.d0*bessel_func(2))
!CALL findtime(.false.,'init')

!CALL findtime(.true.,'loop')

!=========================================================================================
!  Start Chebychev summation loop.
!  Method for each iteration:
!  (1) Operate with the Hamiltonian on phi_nm1, store result in phi_n.
!  (2) Multipy phi_n by 2i.
!  (3) Add phi_nm2 to phi_n.
!  (4) Accumulate terms in the sum along with bessel function coeffs .
!  (5) Move phi_nm1 to phi_nm2.
!  (6) Move phi_n   to phi_nm1.
   Do iterm=3,ncheby

!     phi_n=phi_nm1
!     phi_real=real(phi_n,dp)
!     phi_imag=aimag(phi_n)
      phi_real=real(phi_nm1,dp)
      phi_imag=aimag(phi_nm1)
     Call h_apply( phi_real, irdim, whichir )
     Call h_apply( phi_imag, irdim, whichir )
!     phi_n=2.d0*(-imag)*cmplx(phi_real,phi_imag,dp)+phi_nm2
     phi_n=cmplx(phi_real,phi_imag,dp)
     phi_n=phi_n*2.d0*(-imag)
     phi_n=phi_n+phi_nm2
     phi_in=phi_in+(2.d0*bessel_func(iterm)*phi_n)
     phi_nm2=phi_nm1
     phi_nm1=phi_n

   End Do

!CALL findtime(.false.,'loop')
!=========================================================================================
!  Multiply by leading phase factor

!CALL findtime(.true.,'phase')

   phi_in=Exp(((e_range/2.d0)+eigmin)*deltat*-imag)*phi_in

!CALL findtime(.false.,'phase')
!=========================================================================================
! Apply NIP

! Call findtime(.true.,'nip')

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

        ir_pt=aphindex(ntheta, nchir, irho, itheta, ichi, .false.)

         phi_in(ir_pt) = phi_in(ir_pt)*Exp(-vnip(irho)*deltat)

       ENDDO
     ENDDO
   ENDDO

! Call findtime(.false.,'nip')

! DEALLOCATE( phi_real, phi_imag )
! DEALLOCATE( phi_n, phi_nm1, phi_nm2 )

Return
End Subroutine tcheby