% fdtdWave1D.m - Matthieu Aussal - Copyright 2018 - Licence CC BY-NC-SA 4.0
% CMAP - Ecole Polytechnique, Route de Saclay, 91128 Palaiseau CEDEX France
% For any informations, please contact  : matthieu.aussal@polytechnique.edu

clear all
close all
clc

% Physical parameters
L  = 1;        % String size
x0 = 0.3;      % Initial position
tf = 1;        % Final time
c  = 1;        % Sound celerity

% Initial condition for magnitude and speed
ui = @(x) exp(-(x-x0).^2/1e-2); % gaussian
vi = @(x) zeros(size(x));       % null

% Numerical discretization
dx = 0.01;          % Space step
dt = 0.5*dx/c;      % Time step (CFL)

% Regular grid 
x = (0:dx:L(1))';
N = length(x);

% 1D Laplacian operator
e  = ones(N,1);
Mx = spdiags(e.*[1 -2 1]/(dx^2),[-1,0,1],N,N);
% figure
% spy(Mx)

% % Add Dirichlet condition to operator
% Mx(1,1:2)         = [1 0];
% Mx(end,end-1:end) = [0 1];

% Add Neumann condition to operator
Mx(1,1)     = -1/(dx^2);
Mx(end,end) = -1/(dx^2);

% Initial magnitude and speed (u0 and v0)
u0 = ui(x);
v0 = vi(x);

% Force initial condition
u0(1)   = 0;
u0(end) = 0;

% First step scheme (Forward Euler) 
u1 = u0 + dt * v0;

% Graphical representation
figure(10)
Hfig = plot(x,u1);
axis([0 L -1 1])
grid on

% Video output
vidObj = VideoWriter('1Dneumann.avi');
open(vidObj);

% Order two sheme in times
t = dt;
while t < tf
    % Leap-frog sheme
    u2 = 2*u1 - u0  + (c*dt)^2 * (Mx * u1);
    t  = t + dt;
   
    % Graphical representation
    set(Hfig,'YData',u2)
    title(['Wave propagation t = ',num2str(t,'%6.3f')])
    
    % Write each frame to the file.
    currFrame = getframe(gcf);
    writeVideo(vidObj,currFrame);
    pause(dt)
    
    % Incrementation
    u0 = u1;
    u1 = u2;
end

 % Close the file.
 close(vidObj);

% Error
ref = exp(-(x-(L-x0)).^2/1e-2);   
norm(u2-ref)./norm(ref)


