%linear Schroed. Eq. in 1D solved via Fourier transform.
%iu_t + u_xx -V(x)u= 0

%T. Dohnal
%May 2018, Halle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

all_figures_off=1;   %=1  if no plots are to be shown
save_data =0;        %=1  if the solution vector should be saved (at t = multiples of savefreq)

%data file into which the solution should be written out
%-------------------------------------------------------
file_out = 'DATA_FILES/tmp.mat';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SET-UP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 50000;  %number of spatila grid points (needs to be even)
L = 2000;   %size of the spatial domain
plot_bd = L/2;  %plot the solution between -plot_bd and plot_bd
dt = 0.01;  %temporal step
Tfin = 100;   %final time
savefreq = 10;  %temporal interval between two saves (and plots) of the solution

nsaves = floor(Tfin/savefreq);    %number of saves within the evolution time
Tfin = nsaves*savefreq;
savefreq_steps = ceil(savefreq/dt);  %number of temporal steps between 2 saves
dt = savefreq/savefreq_steps;
nsteps = savefreq_steps*nsaves;       %number of temporal steps within the evolution time


dx = L/N;
x = linspace(-L/2,L/2-dx,N)';   %spatial grid
k_vec = (2*pi/L)*[0:N/2 (-N/2+1):-1]';      %wavenumber array (for FOur. transform)

%potential V
%V = cos(x);
V=zeros(size(x));

%initial data
shift = 0;
phi = exp(-(x-shift).^2);
u=phi;

if(0)
    vel=1;
    u = u.*exp(1i*vel*x);
end


%define a filter for removing radiation that reaches the domain boundary
%(inaccurate way since it generates small reflections going back to the domain but it is often sufficient)
%(one should use a PML)
sl = 2;
rad_filt = (tanh(sl*(x(x<0)-9*(-L/2+shift)/10))+1)/2;
rad_filt = [rad_filt; (-tanh(sl*(x(x>=0)-9*(L/2+shift)/10))+1)/2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%plot the IC
%-------------------------------
if(all_figures_off==0)
    h1=figure(1);
    set(gcf,'Units','normalized'); set(h1,'Position',[.1 .1 .3 .3])
    plot(x,abs(u),'linewidth',2)
    hold on;
    set(gca,'fontsize',14)
    xlabel('x'); ylabel('|u|'); title('t = 0')
    xlim([-plot_bd plot_bd])
    hold off
    if(make_movie)
        Mov(1) = getframe(h1);   %save the figure into an array for a movie generation
    end
end

t=0;
t_all=t;           %array for all the t values
u_all(:,1) = u;    %array for all the solution vectors

%conserved quantities (for checking the code accuracy)
energy(1) = dx*trapz(abs(u).^2);
m=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TIME STEPPING
for m=1:nsaves
    
    %SPLIT-STEP METHOD (Strang splitting)
    %-----------------
    %the second derivative part  iu_t + u_xx =0
    u = fft(u);
    u = exp(-i*k_vec.^2*dt/2).*u;
    u=ifft(u);
    
    for j=1:savefreq_steps-1
        
        %the part: iu_t - V(x)  u  = 0
        u = u.*exp(-1i*V*dt);
        
        %the second derivative part   iu_t + u_xx =0
        u = fft(u);
        u = exp(-i*k_vec.^2*dt).*u;
        u=ifft(u);
        
        
        %radiation elimination
        if(mod(j,20)==0)
            u = u.*rad_filt;
        end
        
        t = t+dt;
    end
    %the part: iu_t - V(x) u  = 0
    u = u.*exp(-1i*V*dt);
    
    %the second derivative part  iu_t + u_xx =0
    u = fft(u);
    u = exp(-1i*k_vec.^2*dt/2).*u;
    u=ifft(u);
    t = t+dt;
    
    
    fprintf('time= %5.3e; final time = %5.2e \n',t,Tfin)
    
    t_all(m)=t;          %array for all the t values
    u_all(:,m) = u;      %array for all the solution vectors
  
    %conserved quantities
    energy(m) = dx*trapz(abs(u).^2);
    
  if(save_data)
      save(file_out,'t_all','u_all','V','x','energy')
    end
    
    if(all_figures_off==0)
        figure(1);
        plot(x,abs(u),'linewidth',2)
        hold on;
        plot(x,abs(phi),'--','linewidth',2)
        xlim([-plot_bd plot_bd])
        hold off
        title(sprintf('t = %g',t))
        if(make_movie)
            Mov(m) = getframe(h1);
        end
        pause(0.1)
    end
end
