NLS Equation using split step method in fortran

[Afroja_A_]

Hi Everyone,

I am trying to solve nonlinear schrodinger equation by spilt step method in fortran95. The equation is : ∂U/∂ξ=- i s/2  ∂^2 U/∂τ^ 2+ i N^ 2 |U |^ 2 U 

here is my code in below:

program nls
use mkl_dfti !, FORGET=>DFTI_DOUBLE, DFTI_DOUBLE=>DFTI_DOUBLE_R

implicit none

real*8 :: length , dispersion, N , m , chirp, nt, Tmax,Input
real*8:: step_num, deltaz, dtau
real*8, parameter :: pi= 3.14159265358979323d0
integer, parameter :: j=1024, k=1024
integer :: rti, o, i
integer :: status
real*8, dimension(1024)  :: omega
real*8, dimension(1024):: tau
complex*16, dimension(1024):: disp ,  temp, f_temp,  Out,Output, F_D_I, F_D_O, u
complex*16 :: hhz, A,B=(0,1)

type(DFTI_DESCRIPTOR), POINTER :: hand

 Status = DftiCreateDescriptor(hand,DFTI_double, DFTI_complex, 1, j)
 Status = DftiSetValue (hand,DFTI_FORWARD_SCALE,1.0d0 )
 Status = DftiSetValue (hand,DFTI_BACKWARD_SCALE, 1.0d0/dble(j))
 Status = DftiCommitDescriptor(hand)

!print*, 'Enter fiber length ='
!read (*,*) length
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
!read (*,*) dispersion
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
!read (*,*) N
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian='
!read (*,*) m
m =0.0d0
 chirp =0.0d0
 
! print*, 'Enter Sampling point nt='
! read(*,*) nt
 nt =1024.0d0
! print*, 'Enter window size Tmax='
! read(*,*) Tmax
Tmax= 32.0d0
 
 
 
 step_num= 1000.0d0*length*(N**2)
 deltaz= length/step_num
 dtau= (2.0d0*Tmax)/nt
 
 
  rti= int(nt)
 
 
 
  !!=============
 
do i= 1,int(nt)
  tau(i)=dble(i-nt/2-1)*dtau

  if (i<= 512) then
   
  omega(i)=dble(i-1)* pi/Tmax
 
  else
  omega(i)=dble(i-nt -1)* pi/Tmax

end if

 !u(i)= 1d0/dcosh(tau(i))
 
end do
 
 u = 1.0d0/dcosh(tau)
 
 hhz= (N**2)*deltaz*B
 A= hhz/2.0d0
 
 do i = 1, 1024
 
 temp(i)=  u(i)* exp(abs(u(i))**2 * A)
 
 end do

 do i =1, 1024
 
 disp(i) = exp(0.50d0*B*dispersion*omega(i)**2*deltaz)
 
 end do
 
 
do i= 1,1024
 open(1, file='Input_time_domain.dat')
 write(1,*) tau(i), abs(u(i))**2
 end do
 
 !!=========================================!!
 

 
!!-==========Main Program===============!!
 

 
 do o = 1, int(step_num)

 Status = DftiComputeBackward (hand,temp)
 
 do i = 1,1024
 
 f_temp(i) = temp(i)*disp(i)
 
 end do
 
 Status = DftiComputeForward (hand,f_temp)
 
 
 
 do i = 1, 1024
 
 Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)
 
 end do
 

 end do
 
 

 
 
 do i = 1, 1024
 
 Output(i)= Out(i)*exp(-abs(f_temp(i))**2*A)
 

 
   open(2, file='Output_time_domain.dat')
   write(2,*) tau(i), abs(Output(i))**2
 
end do

 
  !!========== FRequency domain input and Output file======!!
 
 Status = DftiComputeBackward (hand, u)
 
 
 do i = 1, 1024
 F_D_I(i)  = u(i)*((nt*dtau)/dsqrt(2.0d0*pi))
 end do
 
 
 Status = DftiComputeBackward (hand, Output)
 
 do i = 1,1024
 F_D_O(i)= Output(i)*((nt*dtau)/sqrt(2.0d0*pi))
 end do
 
 do i =1, 1024
 
 open(3, file='Input_Frequency_domain.dat')
 write(3,*) omega(i)/2.0d0*pi, abs(F_D_I(i))**2
 
 
 open(4, file='Output_Frequency_domain.dat')
 write(4,*)omega(i)/(2d0*pi), abs(F_D_O(i))**2

end do

 
 Status = DftiFreeDescriptor (hand)
 
 
 end program nls

===============================

I face problem in propagation calculation portion, where I suppose to go from 1 to step_num, and get dispersion effect. but in my program

as in:  do o = 1, int(step_num)

Status = DftiComputeBackward (hand,temp)
 
 do i = 1,1024
 
 f_temp(i) = temp(i)*disp(i)
 
 end do
 
 Status = DftiComputeForward (hand,f_temp)
 
 
 
 do i = 1, 1024
 
 Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)
 
 end do
 

 end do
 

==== this do loop is i calculate from 1 to 1 in do loop as highlited

program nls
use mkl_dfti !, FORGET=>DFTI_DOUBLE, DFTI_DOUBLE=>DFTI_DOUBLE_R

implicit none

real*8 :: length , dispersion, N , m , chirp, nt, Tmax,Input
real*8:: step_num, deltaz, dtau
real*8, parameter :: pi= 3.14159265358979323d0
integer, parameter :: j=1024, k=1024
integer :: rti, o, i
integer :: status
real*8, dimension(1024)  :: omega
real*8, dimension(1024):: tau
complex*16, dimension(1024):: disp ,  temp, f_temp,  Out,Output, F_D_I, F_D_O, u
complex*16 :: hhz, A,B=(0,1)


type(DFTI_DESCRIPTOR), POINTER :: hand 

 Status = DftiCreateDescriptor(hand,DFTI_double, DFTI_complex, 1, j)
 Status = DftiSetValue (hand,DFTI_FORWARD_SCALE,1.0d0 )
 Status = DftiSetValue (hand,DFTI_BACKWARD_SCALE, 1.0d0/dble(j))
 Status = DftiCommitDescriptor(hand)

!print*, 'Enter fiber length ='
!read (*,*) length
length = 100.0d0
!print*, 'dispersion 1 for normal, -1 for anomalous'
!read (*,*) dispersion
dispersion = 1.0d0
!print*, 'Nonlinear parameter N ='
!read (*,*) N
N= 1.0d0
!print*, 'm = 0 for sech , m>0 for super-gaussian=' 
!read (*,*) m
m =0.0d0
 chirp =0.0d0
 
! print*, 'Enter Sampling point nt='
! read(*,*) nt
 nt =1024.0d0
! print*, 'Enter window size Tmax='
! read(*,*) Tmax
Tmax= 32.0d0
 
 
 
 step_num= 1000.0d0*length*(N**2)
 deltaz= length/step_num
 dtau= (2.0d0*Tmax)/nt
 
 
  rti= int(nt)
 
  
  
  !!=============
  
do i= 1,int(nt)
  tau(i)=dble(i-nt/2-1)*dtau

  if (i<= 512) then
   
  omega(i)=dble(i-1)* pi/Tmax 
  
  else 
  omega(i)=dble(i-nt -1)* pi/Tmax

end if 

 !u(i)= 1d0/dcosh(tau(i)) 
  
end do
 
 u = 1.0d0/dcosh(tau) 
 
 hhz= (N**2)*deltaz*B
 A= hhz/2.0d0
 
 do i = 1, 1024
 
 temp(i)=  u(i)* exp(abs(u(i))**2 * A)
 
 end do

 do i =1, 1024
 
 disp(i) = exp(0.50d0*B*dispersion*omega(i)**2*deltaz)
 
 end do
 
 
do i= 1,1024
 open(1, file='Input_time_domain.dat')
 write(1,*) tau(i), abs(u(i))**2 
 end do
 
 !!=========================================!!
 

 
!!-==========Main Program===============!!
  

 
 do o = 1, int(step_num)

 Status = DftiComputeBackward (hand,temp)
 
 do i = 1,1024
 
 f_temp(i) = temp(i)*disp(i)
 
 end do
 
 Status = DftiComputeForward (hand,f_temp)
 
 
 
 do i = 1, 1024
 
 Out(i) = f_temp(i) *exp(abs(f_temp(i))**2*hhz)
 
 end do
 

 end do
 
 

 
 
 do i = 1, 1024 
 
 Output(i)= Out(i)*exp(-abs(f_temp(i))**2*A)
 

 
   open(2, file='Output_time_domain.dat')
   write(2,*) tau(i), abs(Output(i))**2
 
end do

 

 
  !!========== FRequency domain input and Output file======!!
  
 Status = DftiComputeBackward (hand, u)
 
 
 do i = 1, 1024
 F_D_I(i)  = u(i)*((nt*dtau)/dsqrt(2.0d0*pi))
 end do
 
 
 Status = DftiComputeBackward (hand, Output)
 
 do i = 1,1024
 F_D_O(i)= Output(i)*((nt*dtau)/sqrt(2.0d0*pi))
 end do 
 
 do i =1, 1024
 
 open(3, file='Input_Frequency_domain.dat')
 write(3,*) omega(i)/2.0d0*pi, abs(F_D_I(i))**2 
  
 
 open(4, file='Output_Frequency_domain.dat')
 write(4,*)omega(i)/(2d0*pi), abs(F_D_O(i))**2

end do


 
 Status = DftiFreeDescriptor (hand)
 
 
 end program nls

, then data is ok, but if it increase from 1 to 2 and above it does not give any data. plesae help.

Thanks in advance