MODULE chebyshev2d
	USE fftw_mod	
!	USE hamiltonians, ONLY : h
	IMPLICIT NONE
	PRIVATE :: estimate_emax

CONTAINS
!	SUBROUTINE cheb(psi,t_even,t_odd,time_t,h,n_iter,hbar,calc_emax)
	SUBROUTINE cheb(psi,t_even,t_odd,time_t,h,n_iter,hbar,emin,emax)
		IMPLICIT NONE
		INTERFACE
			SUBROUTINE h(psi,shift,scale)
				USE fftw_mod
				IMPLICIT NONE
				REAL(8), INTENT(IN) :: shift,scale
				TYPE (fftw2d), INTENT(INOUT) :: psi		
			END SUBROUTINE h
		END INTERFACE
		INTEGER(4), INTENT(IN) :: n_iter
		REAL(8), INTENT(IN) :: time_t,hbar
		COMPLEX(8), DIMENSION(:,:), INTENT(INOUT) :: psi
		TYPE (fftw2d), INTENT(INOUT) :: t_even,t_odd		
!		LOGICAL(4), INTENT(IN)  :: calc_emax
		INTEGER(4)  :: n2,n,sgn,c
		REAL(8), INTENT(IN)  :: emin,emax
		REAL(8)  :: shift,scale,k
!		REAL(8), SAVE  :: emin,emax,shift,scale,k
		REAL(8), PARAMETER :: etol=1e-3_8,ctol=1e-18_8

!		IF (calc_emax) THEN
!			CALL estimate_emax(t_even,etol,h,emin,emax)		!  (t_even%inq overwritten)
			shift=(emax+emin)/2
			scale=(emax-emin)/2
!		ENDIF
		
		t_even%inq=psi
		CALL h(t_even,-shift,1/scale)
		t_odd%inq=t_even%outq
		n=1
		sgn=-1
		k=scale*time_t/hbar
		psi=BESJ0(k)*t_even%inq-(0,2)*BESJ1(k)*t_odd%inq

		DO n2=1,n_iter/2
			n=n+1
			CALL h(t_odd,-shift,1/scale)
			t_even%inq=2*t_odd%outq-t_even%inq
			psi=psi+sgn*2*BESJN(n,k)*t_even%inq
!			print*,'n=',n,' bes=',BESJN(n,k)
		
			n=n+1
			CALL h(t_even,-shift,1/scale)
			t_odd%inq=2*t_even%outq-t_odd%inq
			psi=psi-sgn*(0,2)*BESJN(n,k)*t_odd%inq
			
!			if (modulo(n2,100)==0) print*,n,besjn(n,k)
			if (n>k .and. ABS(BESJN(n,k))<ctol) exit
			
			sgn=-sgn
		ENDDO
		write(0,*),'n=',n
		psi=psi*exp((0,-1)*shift*time_t/hbar)
	END SUBROUTINE cheb
	
	SUBROUTINE estimate_emax(psi,etol,h,emin,emax) 
		IMPLICIT NONE
		INTERFACE
			SUBROUTINE h(psi,shift,scale)
				USE fftw_mod
				IMPLICIT NONE
				REAL(8), INTENT(IN) :: shift,scale
				TYPE (fftw2d), INTENT(INOUT) :: psi		
			END SUBROUTINE h
		END INTERFACE
		TYPE(fftw2d), INTENT(INOUT)  :: psi
		REAL(8), INTENT(IN) :: etol
		REAL(8), INTENT(OUT) :: emin,emax
		
		REAL(8) :: shift,norm
		
		norm=estimate_norm(0._8)
		shift=norm
		emax=estimate_norm(shift)-shift
		shift=-emax
		emin=-estimate_norm(shift)-shift
		RETURN
				
	CONTAINS
		REAL(8) FUNCTION estimate_norm(shift)
			REAL(8), INTENT(IN) :: shift
			REAL(8) :: last_norm2,norm2
			
			psi%inq=0
			psi%inq(1,1)=1
			last_norm2=0
			DO
				norm2=dot_product(pack(psi%inq,.true.),pack(psi%inq,.true.))
!				write(0,*) 'norm2=',norm2
				if (abs(last_norm2-norm2)<norm2*etol) then
					estimate_norm=sqrt(norm2)*(1+sqrt(etol))
					return
				end if
				last_norm2=norm2
				psi%inq=psi%inq/sqrt(norm2)
				CALL h(psi,shift,1._8)
				psi%inq=psi%outq
			END DO
		END FUNCTION estimate_norm
	END SUBROUTINE estimate_emax
END MODULE chebyshev2d

