% ---------------------------------------------------------------
%     Comparison of 1st-order upwind and 2nd-order centered diffencesfor 
%       differential-difference approximations to the 
%       constant-windspeed advection eqn in a periodic domain
% ---------------------------------------------------------------

clear all;

% user specified parameters

nx = 30;          % number of grid points (not including the repeated point 
                  % at the begining of the periodic domain).  
min_x = 0;        % coordinate of left edge of computational domain
max_x = 1;        % coordinate of right edge of computational domain

% for isolated spike                  
% t_stop = 5;         % end time for isolated spike
% two_waves = false;  % false gives spike at $x=2/3$
                    
% for 2 waves (7.5 + 10 delta x waves)
t_stop = 12;       % end time for 2-waves
two_waves = true;  % true gives 7.5 + 10 delta x waves

% windowing parameters for the plot

% for isolated spike
% min_tracer = -0.5;   % minimum value for plot of dependent variable
% max_tracer =  1.0;   % maximum value for plot of dependent variable

%for 2 waves
min_tracer = -1.2;   % minimum value for plot of dependent variable
max_tracer =  1.2;   % maximum value for plot of dependent variable

min_x_wndo = 0;   % coordinate of left edge of plotted subdomain
max_x_wndo = 1;   % coordinate of right edge of plotted subdomain

name1 = 'exact';                 % labels for first curve on plot
name2 = '1st-order upstream';    % labels for second curve on plot
name3 = '2nd-order centered';    % labels for third curve on plot
t_anim = 0.5;     % the wall-clock time in seconds between short pauses
                  % in the animation

% ---------------------------------------------------------------
%                     Initialization
% ---------------------------------------------------------------


tspan = linspace(0,t_stop,t_stop+1);    % output times
nt=size(tspan,2);                       % number of output times
exact=zeros(nt,nx);                     % intialize exact solution array

% generate the discrete x axis

dx=(max_x-min_x)/nx;
xpt = linspace(min_x+dx,max_x,nx);
cdx = 1;    % This is speed/(delta x): the wind speed shifts 
            % solution dx in unit time
spd = cdx*dx;  % recover the wind speed

% generate the initial field and the exact solution

if two_waves
    for ix = 1:nx
        initial(ix,1) = 0.5*(sin(2*pi*xpt(ix)/(1/4))+ ...
            sin(2*pi*xpt(ix)/(1/3))) ;
        for it = 1:nt
            exact(it,ix)=0.5*(sin(2*pi*(xpt(ix)-spd*tspan(it))/(1/4))+ ...
                sin(2*pi*(xpt(ix)-spd*tspan(it))/(1/3))) ;
        end
    end
else
    initial = zeros(nx,1);
    initial(15)=1;
%   assumes spike translates integer multiple of dx between display times
    for it = 1:nt   
            exact(it,:)=circshift(initial,tspan(it));
    end
end

% plot initial condition

figure(1);
plot(xpt,initial,'k:');
title('Initial Condition');
axis([min_x_wndo max_x_wndo min_tracer max_tracer]);


% -------------------------------------------------------------
% solve differential-difference equations for all output times
% -------------------------------------------------------------

% 1st-order upstream


[plt_time, tracer1] = ode45(@(t,y)first_biased(t,y,cdx),...
    tspan,initial);

% 2nd-order centered

[plt_time, tracer2] = ode45(@(t,y)second_centered(t,y,cdx),...
    tspan,initial);


% -------------------------------------------------------------
% display animated plot 
% -------------------------------------------------------------

figure(2);
for it=1:nt
    plot(xpt,exact(it,:),'k:', xpt,tracer1(it,:),'r-',xpt,...
        tracer2(it,:),'b--');
    nstr1 = num2str(tspan(nt)); title(['Time = ',nstr1]);
    axis([min_x_wndo max_x_wndo min_tracer max_tracer]);
    pause(t_anim); 
end

% ---------------------------------------------------------------
%                  plot final result with line labels
% ---------------------------------------------------------------

plot(xpt,exact(nt,:),'k:', xpt,tracer1(nt,:),'r-',xpt,tracer2(nt,:),'b--');
nstr1 = num2str(tspan(nt)); title(['Time = ',nstr1]);
axis([min_x_wndo max_x_wndo min_tracer max_tracer]);
legend(name1,name2,name3,0);

