% JDW linearadvection.m 4 Feb. 2016

%Solve the linear advection eqn in 1D on -1 <= x <= 1, with a periodic
%b/cond. Euler method, i.e. forward in time and upwind in space 

%Initial state is a sine wave.

%Advecting velocity is constant, set to unity.

U=1;

% Define x-axis. User chooses the number of nodes, N. I have assumed the index
% ranges 1 <= i <= N with x(1)=-1 and x(N)=1

N=100;
deltax=2.0/(N-1);
x=[-1.0:deltax:1];
%x(0)       % this was a check to make sure the lowest index is not zero by default
x(1);
x(N);

% Initialization. Note that the computed field is held at two times.

theta0=sin(pi*x);
theta_n=theta0;
theta_np1=theta_n;

%plot(x,theta0,'-o') % done to make sure the initial field looks good.
%xlabel('time')
%ylabel('\theta')

%In effect my choice of N dictated the value of delta-x. Instead of
%specifying delta-t I'll specify the Courant number instead -- and deduce
%the timestep

Courant=1.6;
deltat=Courant*deltax/U;

% How many time steps forward will we go?

%nmax=10/deltat

nmax=46;

% we start integrating forward in time
for n=1:nmax
    for i=1:N
       if (i > 1) 
           theta_np1(i)=theta_n(i)-Courant*(theta_n(i)-theta_n(i-1));       
       else
           % in this case i=1, and so i-1=0... which is a non-existent
           % gridpoint. However we use the value at N-1 instead, because
           % the periodic boundary condition applies.
           theta_np1(i)=theta_n(i)-Courant*(theta_n(i)-theta_n(N-1));
       end
    end
    plot(x,theta0,'k-',x,theta_np1);
    xlabel('time')
    ylabel('\theta')
        
    %for an animation (i.e. to see the wave advecting), use:
    refreshdata
    drawnow

    %to see a display that simultaneously plots the solution at ALL times, use:
    
    hold all

    
    %the last operation on each timestep is to feed the new field (_np1)
    %back into the old field
    
    theta_n=theta_np1;
end