!====================================================================d
!   
!     GPE :  
!          Fast Oscillation in time t due to small epsilon
!     
!		i*u_t= H_x/epsilon + beta |u|^2 u
!               u(t=0,x) = u_0(x), 0 < t < TFinal
!
!      Or equivally can be written in slow change variable "s = t/epsilon "
!
!		i*u_s = H_x + epsilon  beta |u|^2 u
!               u(s=0,x) = u_0(x), 0< s < TFinal/epsilon
!
!  Operators :  H_x = -\prt_xx ;
!
!   File related description:
!       split.f90 contains 2nd and 4th order splitting method for slow changing  tau
!
!=====================================================================d
      program main
	use constant
	use fft_mkl
	implicit none 
      complex(8) :: psi(0:nx)
      complex(8) :: psi_gd(0:nx)
      complex(8) :: psit(0:nx)
      complex(8) :: k1(0:nx),k2(0:nx)
      complex(8) :: k3(0:nx),k4(0:nx)

	real(8) :: enkin,enint,entotal
	real(8) :: timef
	real(8) :: time,time_elapsed

	integer(4) :: inds,i,j,k,kth,inloop,ind_total 

      call initial ! x,umx, vtemp 
	call initial_1d(psi)
	call initial_fft_mkl
     
      time = 0.d0

	ind_total = nint(TFinal/dtMacro)
      open(36,file='energy.dat',status='unknown')
	call  energy(enkin,enint,entotal,psi)
      write(36,40) time,entotal,enkin,enint	

	do inds =1,ind_total
		time = time + dtMacro
		call gradient_1d(psi,psi_gd)
		k1 = dtMacro*psi_gd

		psit = psi+0.5d0*k1
		call gradient_1d(psit,psi_gd)
		k2 = dtMacro*psi_gd

		psit = psi+0.5d0*k2
		call gradient_1d(psit,psi_gd)
		k3 = dtMacro*psi_gd

		psit = psi+k3
		call  gradient_1d(psit,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,40) time,entotal,enkin,enint	
		
	enddo 
40	format(t4,4(2x,f16.10))

	kth = 1
	open(10,file='psi.bin',form='unformatted',access='direct',recl=4,Status='replace')
	do j = 0,nx-1
		write(10,rec=kth) dreal(psi(j)),dimag(psi(j))
		kth = kth +1
	enddo
	close(10)


	call free_fft_mkl
      !write(*,*)  'it calculated for',ind_total,'steps'

       end 

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

      complex(8) :: psit(0:nx)

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

	psit = psi
	Status = DftiComputeForward(dfft_x, psit(0:nx-1))
      hxl = (xmax-xmin)
      entotal = 0.0d0
      enkin = 0.d0
      enint = 0.0d0
      do i=0,nx-1
      	t1 = abs(psit(i))**2.d0*vtemp(i)*hxl
      	t3= 0.5d0*beta*eps*vext(i)*abs(psi(i))**4.d0*hx
      	enkin = enkin + t1
      	enint = enint + t3
	enddo
       entotal = enkin + enint

      end

!=====================================================================d

