%
% This code solves the NLS using a fourth order sympletic integrator.
%       i u_t + (1/2) u_xx - alpha*|u|^2 u = 0
%  alpha = 1 is defocusing
%  alpha =-1 is focusing
%  
% Ex:  NLS_instability(n,dt,tf,L,alpha,dispframes,saveframes)
%
% npts=2^n, dt:time step, tf:t_final, domain:[-L,L]
%
%Examples:
%global allt allx allu ran; ran=1
%n=8;dt=0.0025;tf=15;L=20;alpha=-1;dispframes=50;saveframes=200;
%NLS_instability(n,dt,tf,L,alpha,dispframes,saveframes)
%wave_plotsurf(allt,allx,abs(allu),L,1,length(allt),'no')
%wave_plotsurf(allt,allx,imag(allu),L,1,length(allt),'no')
%wave_plotsurf(allt,allx,angle(allu),L,1,length(allt),'no')
%
% DON'T FORGET: 
%global allt allx allu ran; ran=1

function [] = NLS(n,dt,maxtime,L,alpha,dispframes,saveframes)

 global allt allx allu;
 global ran;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking if variables are global
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 if(isempty(ran)==1)
   fprintf('\n');
   error('Error: first do: global allt allx allu ran; ran=1')
 end

 maxstep=maxtime/dt;
% if(exist('allt'))
%  if(isglobal(allt)~=1)
%   fprintf('\n');
%   error('Error: first do: global allt allx allu ran; ran=1')
%  end
% end

 dispscreen= fix(maxstep/100); 

 if(dispframes>0) stopdisp= fix(maxstep/dispframes); 
 else stopdisp=9e+99;
 end
 if(saveframes>0)stopsave= fix(maxstep/saveframes);
 else stopsave=9e+99;
 end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial computational domain defined
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 N=2^n;
 dx0=2*pi/N; x=-pi+(0:N-1)*dx0;
 M = N/2; kvec = fftshift(-M:M-1);
 kvec=kvec.';
 x=x.'; 
 x=x*L/pi;
 allx = x;
 allt = 0:stopsave*dt:maxtime;
 dx=x(2)-x(1);

 plotrealimag=1;
 plotsquare=0;
 u=0; 

 strlength=0;
 kk=0;
 t0=0;
 lw=2;
 fs=16;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% one bright soliton (alpha=-1)
 A = 0.25;
 x0 = -L+5/A;
 phi0=1;
 c = 1.1;
 u1 = A*sech(A*(x-x0)).*exp(i*(c*x+phi0));

% another bright soliton (alpha=-1)
 A = 0.25;
 x0 = -L+12/A;
 phi0=0;
 c2 = 0.75;
 u2 = A*sech(A*(x-x0)).*exp(i*(c2*x+phi0));

% a dark soliton (alpha=+1)
 A = 0.1;
 x03 = -L/2;
 phi0=0;
 k=20;
 c3 = pi*(k+0.5)/L;
 u3 = -A*tanh(A*(x-x03)).*exp(i*(c3*x+phi0));

% u=u1;
% u=u1+u2;
% u=u3;c=c3;

size(u1)
u=u1*0+1+0.001*2*((rand(size(u1))-0.5)+sqrt(-1)*(rand(size(u1))-0.5));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('\n');
 fprintf('L \t\t= %f\n',L);
 fprintf('alpha \t\t= %f\n',alpha);
 fprintf('A \t\t= %f\n',A);
 fprintf('c \t\t= %f\n',c);
 fprintf('\hi_0 \t\t= %f\n',phi0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(1);clf;hold on;
 if(plotsquare)
  plot(x,abs(u).^2,'k')
  high = max(max(abs(u).^2));
  low  = min(min(abs(u).^2));
 else
  plot(x,abs(u),'k','LineWidth',lw)
  high = max(max(abs(u)));
  low  = min(min(abs(u)));
 end
 dhl = (high-low)*0.6+1.0*high;
 title('IC: u(x,t) (black)','FontSize',fs);
 axis([-L L -dhl high+dhl])
 if(plotrealimag)
  plot(x,real(u),'g',x,imag(u),'r','LineWidth',lw)
  title('IC: u(x,t) (black), Real (green) and Imag (red) of u','FontSize',fs);
  axis([-L L -high-dhl high+dhl])
 end

 text(L/1.5,high+dhl/2,['T=' num2str(kk*dt) '/' num2str(maxtime)],'FontSize',fs);
 xlabel('x','FontSize',fs+2)
 ylabel('u(x,t)','FontSize',fs+2)
 set(gca,'FontSize',fs);
 allu = u;

 fprintf('\nt=%5.2f ',t0);
 fprintf('\nPress any key to contiue...');
 pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 beta = (1/3.0)*(2+2^(1/3)+2^(-1/3));
 ooN=1.0/N;
 rho = (pi/L)^2/2.0;
 symbol1 = exp(-i*rho*(kvec.*kvec)*beta*dt);
 symbol2 = exp(-i*rho*(kvec.*kvec)*(1-2*beta)*dt);
 idt1=-i*beta*dt/2.0;
 idt2=-i*(1-2*beta)*dt/2.0;

 for kk=1:maxstep

   [temp1] = NLS_spectral(N,ooN,u,    symbol1,idt1,alpha);
   [temp1] = NLS_spectral(N,ooN,temp1,symbol2,idt2,alpha);
   [u] =     NLS_spectral(N,ooN,temp1,symbol1,idt1,alpha);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   noise=0;
   if(noise>0)
    u  = u+noise*(rand(N,1)-0.5);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/dispscreen)==kk/dispscreen)
    fprintf('%5.2f ',t0+kk*dt);
    strlength=strlength+1;
    if strlength>=10
     strlength=0;
     fprintf('\n');
    end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/stopsave)==kk/stopsave)
    allu = [allu,u];
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
     figure(1); clf;
     if(plotsquare)
      plot(x,abs(u).^2,'k','LineWidth',lw)
     else
      plot(x,abs(u),'k','LineWidth',lw)
     end
     title('u(x,t) (black)','FontSize',fs);
     axis([-L L -dhl high+dhl])
     hold on
     if(plotrealimag)
      plot(x,real(u),'g',x,imag(u),'r','LineWidth',lw)
      title('u(x,t) (black), Real (green) and Imag (red) of u','FontSize',fs);
      axis([-L L -high-dhl high+dhl])
     end

     text(L/1.5,high+dhl/2,['T=' num2str(kk*dt) '/' num2str(maxtime)],'FontSize',fs);
     xlabel('x','FontSize',fs+2)
     ylabel('u(x,t)','FontSize',fs+2)
     set(gca,'FontSize',fs);
     pause(0.01);
   end

 end

 ran=1;
 fprintf('\n');