% ---------------------------------------------------------------
%     Comparison of flux limiters for 
%     constant-windspeed advection in a periodic domain
% ---------------------------------------------------------------

clear all;

% user specified parameters

nx = 50;          % number of grid points (not including the repeated point 
                   %    at the begining of the periodic domain).  
min_x = -1;        % coordinate of left edge of computational domain
max_x =  1;        % coordinate of right edge of computational domain
t_stop = 4;        % end time
spd=1;             % Flow speed (solution will propagate around the domain 
                   %    in 2 time units
cdx=0.5;           % Courant number (will be used to set time step)

square_wave = true;  % square wave if true, isolated peaks if false

% parameters for the plots

t_out_interval = 0.05;  % Interval at which to show animated plots
min_tracer = -0.5;      % minimum value for plot of dependent variable
max_tracer =  1.5;      % maximum value for plot of dependent variable

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

t_anim = 0.2;     % the wall-clock time in seconds between short pauses
                  % in the animation

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

% generate the discrete x axis

dx=(max_x-min_x)/nx;
xpt = linspace(min_x+dx,max_x,nx);

% time-related parameters
dt = cdx*dx/spd;        % recover time step from the Courant number
n_steps = t_stop/dt;    % number of time steps
n_outputs = t_stop/t_out_interval + 1;   % number of output times in animation
tspan = linspace(0,t_stop,n_outputs);    % output times for animation                   

% generate the initial field 

initial=zeros(nx,1);
for ix = 1:nx
    if square_wave
        initial(ix,1) =  ( 1 + square(pi*(xpt(ix)+0.5)) )/2;
    else
        initial(ix,1) = ( (1 + cos(pi*xpt(ix)))/2 )^4;    
    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]);

% set initial values

time=0; 
c_old = initial;
cflim_old = initial;

% -------------------------------------------------------------
% update flux-limiter solution for all output times
% -------------------------------------------------------------

for it=1:n_steps
    
    time=time+dt;
    
% -------------------------------------------------------------
% exact solution
% -------------------------------------------------------------
    exact=zeros(nx,1); 
    for ix = 1:nx
        x_shift= xpt(ix)-spd*time;
        if square_wave
            exact(ix,1) = ( 1 + square(pi*(x_shift+0.5)))/2 ; 
        else
            exact(ix,1) = ( (1 + cos(pi*x_shift))/2 )^4; 
        end
    end 
% -------------------------------------------------------------
% solve without limiting
% -------------------------------------------------------------
 
    c_new =Lax_Wendroff(c_old,cdx,ones(nx,1));
    
% -------------------------------------------------------------
% solve using flux-limiter method -- PUT CALL TO LIMITER HERE
% -------------------------------------------------------------

    Lim = ones(nx,1);
    cflim_new = Lax_Wendroff(cflim_old,cdx,Lim);

% -------------------------------------------------------------
% display animated plots 
% -------------------------------------------------------------
 if any(abs(time-tspan)<1e-5)==1  
    figure(2);
    plot(xpt, exact,'k:',xpt,c_new,'r',xpt,cflim_new,'b');
    nstr1 = num2str(time,3); title(['Tracer at time = ',nstr1]);
    axis([min_x_wndo max_x_wndo min_tracer max_tracer]);
    pause(t_anim); 
 end   
   
    
% -------------------------------------------------------------
% update arrays for next step
% -------------------------------------------------------------
    
    c_old=c_new; 
    cflim_old=cflim_new;
end



