%linear Schroed. Eq. in 1D solved via Fourier transform.
%iu_t + u_xx = 0

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

all_figures_off=0;   %=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 = 20;  %temporal interval between two saves (and plots) of the solution

nsaves = floor(Tfin/savefreq);    %number of saves within the evolution time
Tfin = nsaves*savefreq;

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)

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


%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
    pause(0.4)
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
    u = fft(u);
    u = exp(-1i*k_vec.^2*savefreq).*u;
    u=ifft(u);
    t = t+savefreq;
    fprintf('time= %5.3e; final time = %5.2e \n',t,Tfin)
    
    t_all(m+1)=t;          %array for all the t values
    u_all(:,m+1) = u;      %array for all the solution vectors
    
    %conserved quantities
    energy(m+1) = dx*trapz(abs(u).^2);
    
    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))
        pause(0.4)
    end
end



if(save_data)
    save(file_out,'t_all','u_all','V','x','energy')
end

