% traffic - Program to solve the generalized Burger  

% equation for the traffic at a stop light problem

clear;  help traffic; % Clear memory and print header

N = input('Enter the number of grid points - ');

L = 400;    % System size (meters)

h = L/N;    % Grid spacing for periodic boundary conditions

v_max = 25;    % Maximum car speed (m/s)

fprintf('Suggested timestep is %g\n',h/v_max);

tau = input('Enter time step (tau) - ');

coeff = tau/(2*h);   % Coefficient used by the schemes

coefflw = tau^2/(2*h^2);  % Coefficient used by Lax-Wendroff

%%%%% Initial condition is a square pulse		%%%%% 

%%%%% from x = -L/4 to x = 0               %%%%%

rho_max = 1;                   % Maximum density

Flow_max = 1/4*rho_max*v_max;  % Maximum Flow

rho = zeros(1,N);

for i=N/4:(N/2-1)

  rho(i) = rho_max;    % Max density in the square pulse

end

rho(N/2) = rho_max/2;  % Try running without this line

%%%%% Set loop and plot variables %%%%%

iplot = 1;

xplot = ((1:N)-1/2)*h - L/2;  % Record x scale for plot

rplot(:,1) = rho(:);          % Record the initial state

tplot(1) = 0;

nstep = ceil((L/4)/(v_max*tau)); % Stop after last car moves

%nstep = 5*nstep;  % Go further for last plots %******%

fprintf('Number of iterations = %g\n',nstep);

%%%%% MAIN LOOP %%%%%

ip(1:N) = (1:N)+1;  ip(N) = 1;   % ip = i+1 with periodic b.c.

im(1:N) = (1:N)-1;  im(1) = N;   % im = i-1 with periodic b.c.

for istep=1:nstep

  Flow = rho .* (v_max*(1 - rho/rho_max));

  %%% FTCS method %%%

  rho_new(1:N) = rho(1:N) - coeff*(Flow(ip)-Flow(im));

  %%% Lax method %%%

%  rho_new(1:N) = .5*(rho(ip)+rho(im)) ...

%                  - coeff*(Flow(ip)-Flow(im));

  %%% Lax-Wendroff method %%%

%  cp = v_max*(1 - (rho(ip)+rho(1:N))/rho_max);

%  cm = v_max*(1 - (rho(1:N)+rho(im))/rho_max);

%  rho_new(1:N) = rho(1:N) - coeff*(Flow(ip)-Flow(im)) ...

%   + coefflw*(cp.*(Flow(ip)-Flow(1:N)) ...

%             - cm.*(Flow(1:N)-Flow(im)));

  %%% ----------- %%%

  rho = rho_new;

  iplot = iplot+1;

  rplot(:,iplot) = rho(:);   % Record density for plot

  tplot(iplot) = tau*istep;  % Record time for plot

  plot(xplot,rho,'-',xplot,Flow/Flow_max,'--');

  title('Density (solid) and flow (dashed) versus position');

end

mesh(rplot,[-90 30])  % Wire-mesh plot of density

xlabel('Position'); ylabel('Time');

title('Density versus x and t');

pause;    % Pause between plots; strike any key to continue

% Use rot90 function to graph t vs x since

% contour(rplot) graphs x vs t.

levels = 0:.1:1;    % Contour levels in plot

cs = contour(rot90(rplot),levels,xplot,tplot); 

xlabel('Position');  ylabel('Time');              

clabel(cs,0:0.1:1);       % Put labels on contour levels            

title('Density contours');