PROGRAM bloch2d
	USE chebyshev2d
	USE hamiltonians2d
	IMPLICIT NONE
	INCLUDE "fftw3.f"
	REAL(8) :: qx_period,px_period,qy_period,py_period	
	INTEGER(4) :: n,nx,ny
		
	COMPLEX(8), DIMENSION(nqx,nqy) :: psi
	COMPLEX(8), DIMENSION(nqx,nqy), TARGET :: source_even,source_odd,transform,outcome
	TYPE (fftw2d) :: even,odd	
	REAL(8), PARAMETER :: timestep=3.
	REAL(8) :: time=0._8
	INTEGER(4), PARAMETER :: n_iter=4000,n_steps=5
	REAl(8), PARAMETER :: qx_init=240._8,px_init=0.202526,deltaqx_init=20._8 ! 0<=q,q_init<q_period
	REAl(8), PARAMETER :: qy_init=110._8,py_init=0.,deltaqy_init=20._8 ! 0<=q,q_init<q_period
	REAL(8) :: emin,emax
	
	PRINT*,'2d tight-binding'
	PRINT*,'hop=',hop(-hoprangex:hoprangex,-hoprangey:hoprangey)
	PRINT*,'time=',time,'timestep=',timestep
	PRINT*,'force=',forcex,forcey
	PRINT*,'q_init=',qx_init,qy_init
	PRINT*,'p_init=',px_init,py_init
	PRINT*,'deltaq_init=',deltaqx_init,deltaqy_init
	PRINT*,'barrier height=',bar_height,' width=',bar_width,		&
		' origin = ',bar_originx,bar_originy,' tilt =',bar_alpha
		
	emax=sum(abs(hop))/2+bar_height ! assume height>0
	emin=-sum(abs(hop))/2-forcex*qx_period-forcey*qy_period ! assume force>0

	dqx=potential_periodx
	qx_period=dqx*nqx
	px_period=h_planck*nqx/qx_period
	dpx=px_period/npx
	dqy=potential_periody
	qy_period=dqy*nqy
	py_period=h_planck*nqy/qy_period
	dpy=py_period/npy
		
	FORALL (nx=1:nqx, ny=1:nqy)
		qx(nx,ny)=dqx*(nx-1)
		qy(nx,ny)=dqy*(ny-1)
	END FORALL
	FORALL (nx=1:npx, ny=1:npy)
		px(nx,ny)=dpx*(nx-1)
		py(nx,ny)=dpy*(ny-1)
	END FORALL
	
	even%inq=>source_even
	even%p=>transform
	even%outq=>outcome
	odd%inq=>source_odd
	odd%p=>transform
	odd%outq=>outcome
	
	psi=exp(-torus_dist(qx,qx_init,qx_period)**2/(2*deltaqx_init**2)+ii*2*pi*px_init*qx	&
		-torus_dist(qy,qy_init,qy_period)**2/(2*deltaqy_init**2)+ii*2*pi*py_init*qy)
	psi=psi/(sqrt(norm2(pack(psi,.true.)*dqx*dqy)))
	
!	psi=0
!	psi(2,2)=1.
!	even%inq=psi
!	CALL print_psi()
!	CALL tight2d(even,0.,1.)
!	time=1.
!	even%inq=even%outq
!	even%inq=cshift(cshift(psi,0,2),-1,1)
!	CALL print_psi()
!	stop
!	time=2.
!	even%inq=cshift(cshift(psi,-1,2),0,1)
!	CALL print_psi()
!	stop
	
	even%inq=psi
	CALL print_psi()
	DO n=1,n_steps
		write(0,*) 'time=',time
!		CALL cheb(psi,even,odd,timestep,tight2d,n_iter,hbar,n.eq.1)
		CALL cheb(psi,even,odd,timestep,tight2d,n_iter,hbar,emin,emax)
		time=time+timestep
		even%inq=psi
		CALL print_psi()
	ENDDO
	
	CONTAINS
		REAL(8) FUNCTION norm2(vec)
			IMPLICIT NONE
			COMPLEX(8), INTENT(IN) :: vec(:)
			norm2=DOT_PRODUCT(vec,vec)
		END FUNCTION norm2
		SUBROUTINE print_psi()
			IMPLICIT NONE
			INTEGER(4) n
			
			DO ny=1,nqy
				DO nx=1,nqx
					print "('q_psi(',3f22.12,' ) =',2e25.16e3)",time,qx(nx,ny),qy(nx,ny),even%inq(nx,ny)
				ENDDO
			ENDDO
			PRINT*
		END SUBROUTINE print_psi
		FUNCTION torus_dist(ar,sc,period)
			IMPLICIT NONE
			REAL(8), INTENT(IN) :: sc,period
			REAL(8), DIMENSION(nqx,nqy), INTENT(IN) :: ar
			REAL(8), DIMENSION(nqx,nqy) :: torus_dist
						
			torus_dist=min(ABS(ar-sc),ABS(ar+period-sc),ABS(ar-sc-period))
		END FUNCTION torus_dist
END PROGRAM bloch2d
