module constant
implicit none
real(8),parameter :: pi = 4.d0*datan(1.d0)
complex(8),parameter :: im=dcmplx(0.d0,1.d0)

integer(4),parameter :: nx = 81 ! nx should be odd
integer(4),parameter :: ny = nx ! ny should be odd

real(8) ::  xg(0:nx),gxw(0:nx),ax(0:nx,0:nx) ! Gauss-Hermite 
real(8) ::  yg(0:ny),gyw(0:ny),ay(0:ny,0:ny) ! Gauss-Hermite 

real(8) ::  mux(0:nx),cx(0:nx)
real(8) ::  muy(0:ny),cy(0:ny)
real(8) ::  vtemp(0:nx,0:nx)

real(8),parameter :: beta = 5.d0
real(8),parameter :: eps = 1.d-3
real(8),parameter :: Per = 2.d0*pi

real(8),parameter :: TFinal = Per ! 
real(8),parameter :: dtMacro = Per/1.d3  
real(8),parameter :: dtMicro = Per/(1.d2) 
integer(4),parameter :: order = 4 

integer(4),parameter :: ntmicro = nint(per/dtmicro)! 
integer(4),parameter :: ntmacro = nint(tfinal/dtmacro)! 


real(8),parameter :: w1 = 0.33780179798991440851d0
real(8),parameter :: w2 = 0.67560359597982881702d0
real(8),parameter :: w3 =-0.08780179798991440851d0
real(8),parameter :: w4 =-0.85120719795965763405d0


end module constant

	subroutine initial(psi)
	use constant
	implicit none
	complex(8) :: psi(0:nx,0:ny)

	integer(4) :: i,j

	call zehega(nx+1,xg,cx) ! gx is the Gauss-Hermite points
	call heinit(nx+1,xg,ax,gxw) ! gxw  wights

	call zehega(ny+1,yg,cy) ! gx is the Gauss-Hermite points
	call heinit(ny+1,yg,ay,gyw) ! gxw  wights

	do i= 0,nx
		mux(i) = dble(i)
	enddo
	do j = 0,ny
		muy(j) = dble(j)
	enddo

	do j = 0,ny
		do i = 0,nx

		psi(i,j) = ay(j,0)*ax(i,2)

		!psi(i,j) = ax(i,0)*(ay(j,0)+1.0/2.d0*ay(j,1)+1.d0/4.d0*ay(j,2))+ &
		!ay(j,0)*(ax(i,0)+1.0/2.d0*ax(i,1)+1.d0/4.d0*ax(i,2))
!
		vtemp(i,j) = eps*mux(i) + muy(j)

		enddo
	enddo

	end subroutine initial
