clear all
close all
datestr(now)

%-----------Parameter -----------------------------------   
%
L=2*pi; % x, y in [-L/2, L/2]   %initialize domain length
N=256; % Anzahl der Punkte in [-L/2,L/2] % number of spatial discretization points    

T=2100;   % Anzahl der Zeitschritten t in [0,T]  %number of time steps
h=0.01;  %dt
tt=10;   % Zeitintervall für Bilder-Output %time interval for image output 
c=0.3;   %advection speed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=(L/N)*(-N/2:N/2-1)'; %initialize spatial grid
dx=L/N; 
%
k=((2*pi/L)*[0:N/2 -N/2+1:-1]'); %k array corresponding to the fft transformed array
u0=exp(-2*pi*(x-L/4).^2); %initial condition
%plot initializiation
figure(1)
plot(x,u0,'r','Linewidth', 2)
xlabel('x','FontName','Arial','FontSize',16)
ylabel('u(x,t)','FontName','Times New Roman','FontSize',16)
axis([-L/2 L/2 0 1])
title(['Advection Equation, t=' num2str(0)],'Color','k')
hold on
plot(x,exp(-2*pi*(x-L/4).^2),'b-.','Linewidth', 2)
%legend('numerics','exact','Location','eastoutside','Orientation','vertical')
hold off
pause(0.1)
%
uhat=fft(u0); %fourier transform of initial condition
%
rhs = @(c,k,u) -i*c*k.*u; % RHS function 
%

for n=1:T
    k1 =rhs(c,k,uhat);
    k2 =rhs(c,k,uhat+0.5*h*k1);
    k3 =rhs(c,k,uhat+0.5*h*k2);
    k4 =rhs(c,k,uhat+h*k3);
    uhat = uhat + (h/6)*(k1 + 2*k2 + 2*k3 + k4); %Runge Kutta 4 time stepping
    %
    if mod(n,tt)==0
        u=real(ifft(uhat)); %backtransformation
        %animation
        figure(1)
        plot(x,u,'r','Linewidth', 2)
        xlabel('x','FontName','Arial','FontSize',16)
        ylabel('u(x,t)','FontName','Times New Roman','FontSize',16)
        axis([-L/2 L/2 0 1])
        title(['Advection Equation, t=' num2str(n)],'Color','k')
        hold on
        plot(x,exp(-2*pi*(x-L/4).^2),'b-.','Linewidth', 2)
        %legend('numerics','exact','Location','eastoutside','Orientation','vertical')
        hold off
        pause(0.1)
    end;
end;
