% First order Godunov, solving the traffic flow equation:
%
%     p_t + (p(1-p))_x = 0
%
% with periodic boundary conditions.

% Number of grid cells, dx and dt:

N=200; dx=2*pi/N; dt=dx/4;

% Printing times:

tpr = [0 1 2 4 8 10];

npr=size(tpr,2);

% Number of time intervals in between printing times:

Ndt=round((tpr(2:npr)-tpr(1:npr-1))/dt);

% Auxiliary vectors for looking to the left, middle and right:

L = [N 1:N-1]; M = 1:N; R = [2:N 1]; % (Notice the periodicity here.)

% Colors for the output:

col=['k' 'b' 'r' 'g' 'y' 'c']; nc=length(col);

% Adjusted printing times, since the original ones are typically not multiples of dt.

for jpr=2:npr
  tpr(jpr)=tpr(jpr-1)+Ndt(jpr-1)*dt;
end

% Position at the middle of each cell:

x=dx/2:dx:2*pi;

% Initial data:

p = 0.25*(1.5+sin(x-pi));

% Plot initial data.

figure(1)

plot(x,p,col(1)); hold on

% Main loop over printing times.

for jpr=2:npr

% Loop over time steps.

  for j = 1:Ndt(jpr-1)
      
% Characteristic speed.

    c = 1-2*p;

% Fluxes Q: from the right, from the left, 
% or at x/t=0 within a rarefaction, yielding c=0 => p=1/2 => Q=1/4.

    Fl = p.*(1-p);

    Q = (c(L)>0).*(c(M)>0).*Fl(L) + ...
        (c(L)<=0).*(c(M)<=0).*Fl(M) + ...
        (c(L)>0).*(c(M)<=0).*(((c(L)+c(M))>0).*Fl(L)+((c(L)+c(M))<=0).*Fl(M)) +...
        (c(L)<=0).*(c(M)>0)*0.25;
    
% (The third line corresponds to a shock moving with either positive or
% negative speed.)

% Update of p using conservation.

    p = p + dt*(Q(M)-Q(R))/dx;

  end

% Plot at printing time jpr.

  plot(x,p,col(mod(jpr-1,nc)+1))

end

xlabel('x'); ylabel('\rho'); 
title('Snapshots of a numerical solution to the traffic flow equations')
hold off


