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 = 79 ! nx should be odd
real(8) ::  xg(0:nx),gxw(0:nx),ax(0:nx,0:nx) ! Gauss-Hermite 
real(8) ::  mu(0:nx),cx(0:nx)

real(8),parameter :: beta = 1.d0
real(8),parameter :: eps_para = 2.d0**dble(-3)
real(8),parameter :: Per = 2.d0*pi

real(8),parameter :: TFinal = per  
real(8),parameter :: dtMacro = pi/2.d0**dble(12) ! in tau variable
real(8),parameter :: dtMicro = pi/2.d0**dble(5)   ! in tau variable
integer(4),parameter :: order = 4 ! gradient approximation order

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

real(8),parameter :: w1 = 1.d0/4.d0/(2.d0-2.d0**(1.d0/3.d0))
real(8),parameter :: w2 =  1.d0/2.d0/(2.d0-2.d0**(1.d0/3.d0))
real(8),parameter :: w3 = (1.d0-2.d0**(1.d0/3.d0))/4.d0/(2.d0-2.d0**(1.d0/3.d0))
real(8),parameter :: w4 = -2.d0**(1.d0/3.d0)/(2.d0-2.d0**(1.d0/3.d0))/2.d0

end module constant

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

	integer(4) :: j

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

	do j= 0,nx
		psi(j) = ax(j,0) + ax(j,1) 
		mu(j) = dble(j)
	enddo

	end subroutine initial
