!   To compute the dynamics of 1D GPE under harmonic oscillator
!
!	-1/2 \p_xx  + 1/2[(x/eps)^2-1]/eps^2} u + \beta |u|^2 u,  0<t<T_0
!	
!	By change of function wu(\bx,t) = u(x,t), \bx =  x/eps, wu statisfies equation 
!
!             i u_t =1/eps/eps [-1/2 \p_xx  + 1/2 (x*x-1) ] u + \beta |u|^2 u,  0<t<T_0
!	
!	or im micro window tu(x,s) = wu(x,t) , s = t/eps
!
!             i u_s =1/2(-\p_xx + (x^2-1)) u + eps*eps \beta |u|^2 u, 0< s <T_0/eps/eps
!
!   Numerical method: 
!   marco window: Runge kutta fourth-order 
!   micro window: fourth-order time-splitting hermite -spectral method
!   approximation order: p = 2,4 or 6
!
!  Note: nz must be ODD !
!
!*************************************************************************** 
!  Variables:
!
! 'cc2': constant dimension for second order derivative, i.e.
! 'pi': constant \Pi
! 'dt':  time step
! 're': [0, re] solution interval 
! 'tmax': end time
!
	program main
	use constant
	implicit none
	complex(8) :: psi(0:nx),psi_gd(0:nx)
	real(8) :: fr(0:nx),fm(0:nx)
	real(8) :: fr1(0:nx),fm1(0:nx)
	complex(8) :: psiStored(0:nx,0:ntmacro)

	integer(4) :: i,j,k
	real(8) :: enkin,enint,entotal,time,mr
	complex(8) :: t

	call initial(psi) ! psi,mu,xg,gxw,ax
	
	do j = 0,nx
		psiStored(j,0) = psi(j)
	enddo

	time = 0.d0

	call normStat(psiStored(:,0),mr)
	write(*,*) time,mr

	call energy(enkin,enint,entotal,psiStored(:,0))
      write(36,40) time,entotal,enkin,enint	
      !write(*,40) time,entotal,enkin,enint	

	do k = 0,ntmacro-1

		time = time + dtMacro

		do j=0,nx
			psi(j)= exp(-im*beta*eps*2.d0*w1*dtmacro*abs(psi(j))**2.d0)*psi(j)
		enddo

		do j = 0,nx
			fr(j) = dreal(psi(j))
			fm(j) = dimag(psi(j))
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
		do j=0,nx
			t = exp(-im*2.d0*w2*dtmacro*mu(j))*dcmplx(fr1(j),fm1(j))
			fr(j) = dreal(t)
			fm(j) = dimag(t)
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
		do j = 0,nx
			psi(j) = dcmplx(fr1(j),fm1(j))
		enddo

		do j=0,nx
			psi(j) = exp(-im*beta*eps*2.d0*w3*dtmacro*abs(psi(j))**2.d0)*psi(j)
		enddo

		do j = 0,nx
			fr(j) = dreal(psi(j))
			fm(j) = dimag(psi(j))
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
		do j=0,nx
			t = exp(-im*2.d0*w4*dtmacro*mu(j))*dcmplx(fr1(j),fm1(j))
			fr(j) = dreal(t)
			fm(j) = dimag(t)
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
		do j = 0,nx
			psi(j) = dcmplx(fr1(j),fm1(j))
		enddo

		do j=0,nx
			psi(j) = exp(-im*beta*eps*2.d0*w3*dtmacro*abs(psi(j))**2.d0)*psi(j)
		enddo

		do j = 0,nx
			fr(j) = dreal(psi(j))
			fm(j) = dimag(psi(j))
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
		do j=0,nx
			t = exp(-im*2.d0*w2*dtmacro*mu(j))*dcmplx(fr1(j),fm1(j))
			fr(j) = dreal(t)
			fm(j) = dimag(t)
		enddo
		call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
		call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
		do j = 0,nx
			psi(j) = dcmplx(fr1(j),fm1(j))
		enddo

		do j=0,nx
			psi(j) = exp(-im*beta*eps*2.d0*w1*dtmacro*abs(psi(j))**2.d0)*psi(j)
		enddo

		do j=0,nx
			psiStored(j,k+1) =psi(j)
		enddo
		
		call energy(enkin,enint,entotal,psiStored(:,k+1))
		call  normStat(psiStored(:,k+1),mr)

      	write(36,40) time,entotal,enkin,enint	
      	write(*,40)  time,entotal,enkin,enint	
		write(*,*) time, mr

	enddo
40	format(t4,4(2x,f16.10))

	open(10,file='psi.bin',form='unformatted',access='direct',recl=4*(nx+1),Status='replace')
	open(20,file='psi_spec.bin',form='unformatted',access='direct',recl=4*(nx+1),Status='replace')
	kth = 1
	do k =0,ntmacro
		write(10,rec=kth) psiStored(:,k)	
		write(20,rec=kth) psiSpectralStored(:,k)	
		kth = kth +1
	enddo
	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)
	real(8),intent(out) :: enkin,enint,entotal

      complex(8) :: psit(0:nx)
      real(8) :: fr(0:nx),fm(0:nx)
      real(8) :: fr1(0:nx),fm1(0:nx)

	integer(4) :: i,j
	real(8) :: t1,t3

      enkin = 0.d0
      enint = 0.0d0

	psit = psi
	do j = 0,nx
		fr(j) = dreal(psi(j))
		fm(j) = dimag(psi(j))
	enddo

	call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
	call hetrfm1(nx+1,fm,fm1,ax,gxw,0)

	do j = 0,nx
	enkin = enkin + mu(j)*(fr1(j)*fr1(j)+fm1(j)*fm1(j))
	enddo

	do j = 0,nx
		fr(j) = abs(psi(j))**2.d0
	enddo
	call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
	do j = 0,nx
		enint = enint + 0.5d0*eps*beta*fr1(j)*fr1(j)
	enddo

      entotal = enkin + enint

      end


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

	real(8) :: fr(0:nx),fr1(0:nx)
	integer(4) :: j

	open(10,file='psi.dat',status='unknown')
	do j = 0,nx
		write(10,*) dreal(psi(j))
	enddo
	close(10)

	open(20,file='psiback.dat',status='unknown')
	do j = 0,nx
		fr(j) = dreal(psi(j))
	enddo
	call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
	call hetrfm1(nx+1,fr1,fr,ax,gxw,1)

	do j = 0,nx
		write(20,*) fr(j)
	enddo
	close(20)
	stop


	end subroutine testcode

	subroutine  normStat(psi,mr)
	use constant
	implicit none
	real(8) :: mr
	complex(8) :: psi(0:nx)

	integer(4) :: i,j
	real(8) :: fr(0:nx),fm(0:nx)
	real(8) :: fr1(0:nx),fm1(0:nx)

	do j = 0,nx
		fr(j) = dreal(psi(j))
		fm(j) = dimag(psi(j))
	enddo
	call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
	call hetrfm1(nx+1,fm,fm1,ax,gxw,0)

	mr = 0.d0
	do j = 0,nx
		mr = mr + fr1(j)**2.d0+fm1(j)**2.d0 
	enddo

	end subroutine normStat