	program main
	use constant
	implicit none
	complex(8) :: psi(0:nx,0:ny),psi_gd(0:nx,0:ny)
	complex(8) :: k1(0:nx,0:ny),k2(0:nx,0:ny),k3(0:nx,0:ny),k4(0:nx,0:ny)

	integer(4) :: i,j,k,kth
	real(8) :: enkin,enint,entotal,time

	call initial(psi) ! psi,mu,xg,gxw,ax
	
	time = 0.d0

	open(36,file='energy.bin',form='unformatted',access='direct',recl=8,Status='replace')
	open(10,file='modex.bin',form='unformatted',access='direct',recl=64,Status='replace')
	open(20,file='modey.bin',form='unformatted',access='direct',recl=64,Status='replace')

	kth = 1
	call energy(enkin,enint,entotal,psi)
	write(36,rec=kth) time,entotal,enkin,enint

	call hetrfm2_complex(psi,psi_gd,0)
	write(10,rec=kth) psi_gd(0:15,0)
	write(20,rec=kth) psi_gd(2,0:15)
	kth = kth+1



	do k = 0,ntmacro-1

		time = time + dtMacro

		call gradient(psi,psi_gd)
   		k1 = dtmacro*psi_gd

		call gradient(psi+0.5*k1,psi_gd)
		k2 = dtmacro*psi_gd

   		call gradient(psi+0.5*k2,psi_gd)
   		k3 = dtmacro*psi_gd

   		call  gradient(psi+k3,psi_gd)
   		k4 = dtmacro*psi_gd

   		psi= psi  + (k1+2.d0*k2+2.d0*k3+k4)/6.d0
		
		call energy(enkin,enint,entotal,psi)
		write(36,rec=kth) time,entotal,enkin,enint
      	!write(36,40) time,entotal,enkin,enint	
      	!write(*,40) time,entotal,enkin,enint	
		!write(*,*) time
!
		call hetrfm2_complex(psi,psi_gd,0)
		write(10,rec=kth) psi_gd(0:15,0)
		write(20,rec=kth) psi_gd(2,0:15)
		kth = kth+1
	enddo

40	format(t4,4(2x,f16.10))
	!close(36)
	close(10)
	close(20)



	end program main

      subroutine  energy(enkin,enint,entotal,psi)
!======================================================================d
	use constant
      implicit none
      complex(8),intent(in) :: psi(0:nx,0:ny)
	real(8),intent(out) :: enkin,enint,entotal

      complex(8) :: psit(0:nx,0:ny)
      complex(8) :: psis(0:nx,0:ny)

	integer(4) :: i,j

      enkin = 0.d0
      enint = 0.0d0

	call hetrfm2_complex(psi,psit,0)
	do j = 0,ny
		do i = 0,nx
			enkin = enkin+ abs(psit(i,j))**2.d0*(mux(i)+muy(j)/eps)
		enddo
	enddo

	do j = 0,ny
		do i = 0,nx
		psit(i,j) = abs(psi(i,j))**2.d0
		enddo
	enddo
	call hetrfm2_complex(psit,psis,0)
	do j = 0,ny
		do i = 0,nx
		enint = enint + 0.5d0*beta*abs(psis(i,j))**2.d0
		enddo
	enddo

      entotal = enkin + enint

      end

