%
%This code solves the NLS+Potential using a fourth order sympletic integrator.
%
%       i u_t + (1/2) u_xx - alpha*|u|^2 u = V(x) u
%
%with:
%
%  V(x) = V0*tanh(x)^2;
%  alpha = 1 is defocusing
%  alpha =-1 is focusing
%  
%Ex:  NLS_potential(n,dt,tf,L,vel,V0,alpha,dispframes,saveframes)
%
% npts=2^n, dt:time step, tf:t_final, domain:[-L,L]
% vel=velocity of soliton, V0=heith of tanh^2 potential
%
%Example:
% global allt allx allu v ran; ran=1
% n=7;dt=0.1;tf=100;L=10;alpha=-1;vel=0.2;V0=0.1;dispframes=50;saveframes=200;
% NLS_potential(n,dt,tf,L,vel,V0,alpha,dispframes,saveframes)
% wave_plotsurf(allt,allx,abs(allu).^2,L,1,length(allt),'no')
%
%DON'T FORGET: 
% global allt allx allu v ran; ran=1

function [] = NLS_potential(n,dt,maxtime,L,vel,V0,alpha,dispframes,saveframes)

 global allt allx allu v;
 global ran;

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

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

 maxstep=maxtime/dt;
 if(exist('allt'))
  if(isglobal(allt)~=1)
   fprintf('\n');
   error('Error: first do: global allt allx allu v 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=1;
 u=0; 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs and potential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 x0 = L/10;
 phi0=0;
 ip=[2:N,N]; im=[1,1:N-1];
 mu=1;

 vBS = V0*tanh(x).^2;
 uBS = sqrt(1-V0)*sech((x-x0)).*exp(1i*(vel*x+phi0));
 %uBS = sech(x-x0).*exp(1i*(vel*x+phi0));

 vDS = 0.5*V0^2*x.^2;
 uDS = tanh(x-x0);
 uTF = sqrt(max(mu-v,0));
 %smooth TF:
 for ii=1:15; uTF = (uTF(ip)+uTF(im))./2; end;uTF=mu*uTF./max(uTF);

 u = uBS;
 %u = uTF.*uDS;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('\n');
 fprintf('L \t\t= %f\n',L);
 fprintf('alpha \t\t= %f\n',alpha);
 fprintf('V0 \t\t= %f\n',V0);
 fprintf('vel \t\t= %f\n',vel);
 fprintf('\phi_0 \t\t= %f\n',phi0);

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

 text(L/1.75,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_potential_spectral(N,ooN,u,    v,symbol1,idt1,alpha);
   [temp1] = NLS_potential_spectral(N,ooN,temp1,v,symbol2,idt2,alpha);
   [u] =     NLS_potential_spectral(N,ooN,temp1,v,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',x,v,'b','LineWidth',lw)
     else
      plot(x,abs(u),'k',x,v,'b','LineWidth',lw)
     end
     title('u(x,t) (black), v(x) (blue)','FontSize',fs);
     axis([-L L -dhl high+dhl])
     hold on
     if(plotrealimag)
      plot(x,real(u),'g',x,imag(u),'r',x,v,'b','LineWidth',lw)
      title('u (black), Real (green) and Imag (red), v(x) (blue)','FontSize',fs);
      axis([-L L -high-dhl high+dhl])
     end

     text(L/1.75,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');
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pseudo spectral integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u] = NLS_potential_spectral(N,ooN,u,v,symbol,idt,alpha)

                mag = alpha*u.*conj(u)+v; %-0*0.1*1i;
                u = exp(idt*mag).*u;

                uh = fft(u)*ooN;
                udh = symbol.*uh;
                u = ifft(udh)*N;

                mag = alpha*u.*conj(u)+v; %-0*0.1*1i;
                u = exp(idt*mag).*u;

return