MODULE fftw_mod
	IMPLICIT NONE
	REAL(8), PARAMETER :: h_planck=1.,pi=4*ATAN(1._8),hbar=h_planck/(2*pi)
	COMPLEX(8), PARAMETER :: ii=(0,1)
	TYPE fftw 
		COMPLEX(8), POINTER :: inq(:)
		COMPLEX(8), POINTER :: p(:)
		COMPLEX(8), POINTER :: outq(:)
		INTEGER(8) :: planf,planr
	END TYPE
	TYPE fftw2d 
		COMPLEX(8), POINTER :: inq(:,:)
		COMPLEX(8), POINTER :: p(:,:)
		COMPLEX(8), POINTER :: outq(:,:)
		INTEGER(8) :: planf,planr
	END TYPE
END MODULE fftw_mod


MODULE hamiltonians
	USE fftw_mod
	IMPLICIT NONE
	INTEGER(4), PARAMETER :: nq=1024,np=nq
	REAL(8), PARAMETER :: potential_period=1._8
	COMPLEX(8), PARAMETER :: hop=2.75,hop_even=2.,hop_odd=1.
	REAL(8), PARAMETER :: force=1E-2_8,site_even=1,site_odd=-1
	REAL(8), PARAMETER :: bar_height=0.,bar_width=5,bar_origin=512
	REAL(8), DIMENSION(nq) :: qn,pn
	REAL(8) :: dq,dp
CONTAINS	
	SUBROUTINE tight1d(psi,shift,scale) ! Assume nc=1
		IMPLICIT NONE
		REAL(8), INTENT(IN) :: shift,scale
		TYPE (fftw), INTENT(INOUT) :: psi
		
		psi%outq=conjg(hop/2)*cshift(psi%inq,1)+(hop/2)*cshift(psi%inq,-1)	&
			-qn*force*psi%inq+psi%inq*bar_height/COSH((qn-bar_origin)/bar_width)
		psi%outq=scale*(psi%outq+shift*psi%inq)
	END SUBROUTINE tight1d
	SUBROUTINE dp_tight1d(psi,shift,scale)
		IMPLICIT NONE
		REAL(8), INTENT(IN) :: shift,scale
		REAL(8), DIMENSION(nq/2) :: psi_even,psi_odd
		TYPE (fftw), INTENT(INOUT) :: psi
		
		psi_even=psi%inq(2:nq:2)
		psi_odd=psi%inq(1:nq-1:2)
		psi%outq(1:nq:2)=conjg(hop_odd)*psi_odd+(site_odd-qn(1:nq-1:2)*force)			&
			*psi_odd+hop_even*cshift(psi_odd,-1)
		psi%outq(2:nq:2)=conjg(hop_even)*cshift(psi_odd,1)+(site_even-qn(2:nq:2)*force)	&
			*psi_even+hop_odd*psi_odd
		psi%outq=scale*(psi%outq+shift*psi%inq)
	END SUBROUTINE dp_tight1d
END MODULE hamiltonians

MODULE hamiltonians2d
	USE fftw_mod
	IMPLICIT NONE
	INTEGER(4), PARAMETER :: nqx=300,nqy=256,npx=nqx,npy=nqy
	REAL(8), PARAMETER :: potential_periodx=1._8,potential_periody=1._8
	INTEGER(4), PARAMETER :: hoprangex=1,hoprangey=1		!  H=sum hop(ix,iy) cos(2pi(ix px+ iy py)+
	COMPLEX(8), DIMENSION (-hoprangex:hoprangex,-hoprangey:hoprangey) :: hop    	
	DATA hop(0,0),hop(1,0),hop(-1,0),hop(0,1),hop(0,-1)/0,2.,2.,2.,2./ !		+ qx forcex + qy forcey
	DATA hop(1,1),hop(1,-1),hop(-1,1),hop(-1,-1)/0.4,0.4,0.4,0.4/ 
	REAL(8), PARAMETER :: forcex=9E-5,forcey=0.
	REAL(8), PARAMETER :: bar_height=3.5,bar_width=5,bar_originx=128,bar_originy=128,bar_alpha=0	 !  + sech barrier
	REAL(8), DIMENSION(nqx,nqy) :: qx,px,qy,py
	REAL(8) :: dqx,dqy,dpx,dpy
CONTAINS	
	SUBROUTINE tight2d(psi,shift,scale) 
		IMPLICIT NONE
		INTEGER(4) :: ix,iy
		REAL(8), INTENT(IN) :: shift,scale
		TYPE (fftw2d), INTENT(INOUT) :: psi
				
		psi%outq=-(qx*forcex+qy*forcey)*psi%inq
		psi%outq=psi%outq+psi%inq*bar_height/	&
			COSH(((qx-bar_originx)*COS(bar_alpha)+(qy-bar_originy)*SIN(bar_alpha))/bar_width)
		do  iy=-hoprangey,hoprangey
			do ix=-hoprangex,hoprangex
				psi%outq=psi%outq+hop(ix,iy)/2*cshift(cshift(psi%inq,iy,2),ix,1)
			end do
		end do
		psi%outq=scale*(psi%outq+shift*psi%inq)
	END SUBROUTINE tight2d
END MODULE hamiltonians2d
