Subroutine zztcheby(irdim, 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)
!=========================================================================================
!     A L L O C A T A B L E

!REAL(dp),ALLOCATABLE  :: phi_real_temp(:), phi_imag_temp(:)
!REAL(dp),ALLOCATABLE  :: phi_real_n(:), phi_real_nm1(:), phi_real_nm2(:)
!REAL(dp),ALLOCATABLE  :: phi_imag_n(:), phi_imag_nm1(:), phi_imag_nm2(:)
!REAL(dp),ALLOCATABLE  :: phi_real_sum(:), phi_imag_sum(:)

!=========================================================================================

! ALLOCATE( phi_real_temp(irdim), phi_imag_temp(irdim) )
! ALLOCATE( phi_real_sum(irdim), phi_imag_sum(irdim))
! ALLOCATE( phi_real_n(irdim), phi_real_nm1(irdim), phi_real_nm2(irdim) )
! ALLOCATE( phi_imag_n(irdim), phi_imag_nm1(irdim), phi_imag_nm2(irdim) )

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

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

  !phi_real_n=real(phi_in)
  !phi_imag_n=aimag(phi_in)

  phi_real_nm2=phi_real_n
  phi_imag_nm2=phi_imag_n


  phi_real_sum=phi_real_nm2*bessel_func(1)
  phi_imag_sum=phi_imag_nm2*bessel_func(1)

  Call h_apply( phi_real_n, irdim, whichir )
  Call h_apply( phi_imag_n, irdim, whichir )

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

   phi_real_nm1=phi_imag_n
   phi_imag_nm1=-1.d0*phi_real_n

   phi_real_sum=phi_real_sum+(phi_real_nm1*2.d0*bessel_func(2))
   phi_imag_sum=phi_imag_sum+(phi_imag_nm1*2.d0*bessel_func(2))

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

!=========================================================================================
!  Start Chebychev summation loop.

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

   Do iterm=3,ncheby

     phi_real_n=phi_real_nm2
     phi_imag_n=phi_imag_nm2
     phi_real_nm2=phi_real_nm1
     phi_imag_nm2=phi_imag_nm1
     Call h_apply( phi_real_nm1, irdim, whichir )
     Call h_apply( phi_imag_nm1, irdim, whichir )
     phi_real_n=phi_real_n+phi_imag_nm1*2.d0
     phi_imag_n=phi_imag_n-phi_real_nm1*2.d0

     phi_real_sum=phi_real_sum+(2.d0*bessel_func(iterm)*phi_real_n)
     phi_imag_sum=phi_imag_sum+(2.d0*bessel_func(iterm)*phi_imag_n)

     phi_real_nm1=phi_real_n
     phi_imag_nm1=phi_imag_n

   End Do

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

!=========================================================================================
!  Multiply by leading phase factor

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

 phi_real_n=(COS(((e_range/2.d0)+eigmin)*deltat)*phi_real_sum) &
           +(SIN(((e_range/2.d0)+eigmin)*deltat)*phi_imag_sum)
 phi_imag_n=(COS(((e_range/2.d0)+eigmin)*deltat)*phi_imag_sum) &
           -(SIN(((e_range/2.d0)+eigmin)*deltat)*phi_real_sum)

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

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

 phi_real_n=phi_real_n*vnip_ir
 phi_imag_n=phi_imag_n*vnip_ir

 !phi_in=Cmplx(phi_real_n,phi_imag_n,dp)


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

Return
End Subroutine zztcheby