
clear all; clc;
ff =0;

%% Configure source wavelet
nu0 = 10;    %Hz

%% Evaluate wavelet and plot it
ts = linspace(0, 0.5, 1000);
ws = ricker_wave(ts, nu0);

figure;
plot(ts, ws, '-k', 'linewidth', 3);
xlabel('t'); ylabel('w(t)'); title('Ricker Wavelet');


%% Domain Parameters
x_limits = [ 0.0, 1.0 ];
nx = 201;
xx = linspace(x_limits(1), x_limits(2), nx);
dx = ( x_limits(2) - x_limits(1) ) / (nx - 1);

%% Source Parameter
x_s = 0.1;

%% Load the model
[C, C0] = basic_model( x_limits, nx, dx );

%% Build an example set of matrices
[M, A, K] = construct_matrices(C, nx, dx);

%% Set CFL safety constant
alpha = 1/6;

%% Define timestep parameters
T = 3;   % seconds
dt = alpha * dx / max(C);
nt = round(T/dt);

%% Generate the sources
for i = 0:nt-1
    t = i*dt;
    f = point_source( ricker_wave(t,nu0), x_s, nx, dx);
    sources(i+1,:) = f;
end

%% Generate wavefields
u = zeros(nt, nx);
tt = [1/nt:T/nt:T]; % time plot

figure;
for i = 1:nt-1
    ud = [0, u(i,1:end-1)];
    un = u(i,:);
    up = [u(i,2:end),0];
    
    if i == 1
        u_past = u(i,:);
    else
        u_past = u(i-1,:);
    end
    
    u(i+1,:) = M*dt^2 .* ( sources(i,:) + ( (up - 2*un + ud) .* K ) ) + 2*un - u_past;
    u(i+1,1) = A(1) * K(1) * dt * (up(1) - un(1)) + un(1);
    u(i+1,end) = A(end) * K(end) * dt * (-un(end) + ud(end)) + un(end);
    
    
    %% Plotting and saving the movie
    if i == 1
        subplot(1,2,1);
        hh = plot(xx,u(i+1,:));
        ylim([-1e-2 1e-2]);
        title_label = [ 'Wave Propagation at t = ',sprintf('%.2f',tt(i)), ' s' ];
        title(title_label, 'FontSize', 20);
        xlabel('Position (m)', 'FontSize', 18);
        ylabel('Amplitude', 'FontSize', 18);
        
        subplot(1,2,2);
        h = imagesc(xx,tt,u);
        colormap gray;
        colorbar;
        caxis([-8e-3,8e-3]);
        set(gca,'Ydir','normal');
        set(h,'EraseMode','none');
        title('Position and Time Record', 'FontSize', 20);
        xlabel('Position (m)', 'FontSize', 18);
        ylabel('Time (s)', 'FontSize', 18);
        
        
    else
        subplot(1,2,1);
        title_label = [ 'Wave Propagation at t = ',sprintf('%.2f',tt(i)), ' s' ];
        title(title_label, 'FontSize', 20);
        
        set(hh, 'XData', xx, 'YData', u(i+1,:))
        pause(0.0005);
        set(h,'CData',u)
        drawnow;
        
    end
    
    %% Saving pngs for the movie
%     if mod(i+19,20) == 0
%     ff=ff+1;
%     saveas(gcf, ['wave_propagation_',sprintf('%04d',ff),'.png'])
%     end
    
end


%% Plotting final image
figure;
set(gca,'Ydir','normal')
h = imagesc(xx,tt,u);
colormap gray;
colorbar;
caxis([-8e-3,8e-3]);
set(gca,'Ydir','normal');
title('Position and Time Record', 'FontSize', 20);
xlabel('Position (m)', 'FontSize', 18);
ylabel('Time (s)', 'FontSize', 18);



% saveas(gcf, ['wave_propagation_final.png'])

%% Plotting one trace at r = 0.15
figure;
trace_pos = 0.15;
plot(tt',u(:,round(  ( trace_pos-x_limits(1) )*nx ) ), 'linewidth', 3 );
xlabel('time (s)', 'Fontsize',18); ylabel('Amplitude', 'Fontsize',18); title('Trace at r = 0.15', 'Fontsize',20);

ylim([-1e-2 1e-2]);
pos = get(gcf, 'Position');
% saveas(gcf, ['trace_0.15.png']);