% traffic simulation

% simulate traffic circulating around a one-lane race track
function traffic()
  % low values of numcells seem to reveal that this code still
  % needs an entropy fix to handle rarefaction waves;
  % I suspect that second-order corrections would constitute
  % an entropy fix (or provide enough smearing to make one unnecessary).
  numcells=80;
  % initial density distribution of cars
  Phigh=9.8;
  Plow=1.2;
  Pave=0.5*(Phigh+Plow);
  Pamp=0.5*(Phigh-Plow);
  P=Pamp*cos((1:numcells)*(2*pi/numcells))+Pave;
  Pmax_initial = max(P);
  plot(P);
  nframe=1;
  dx=1;
  dt=1;
  dt_dx = dt/dx;
  %desired_CFL=1;
  q=-1; % type q to quit
  input_data=0;
  while ~(numel(input_data)==0 && isequal(input_data,-1))
    if isempty(input_data)
      input_data = nframe;
    end
    if (input_data==q)
      return;
    end
    plot(P);
    axis([0 numcells 0 10]);
    P=step_time(P,dt_dx);
    nframe=nframe+1;
    input_data=input(['plot frame ' num2str(nframe) '?']);
  end
end

function P_out=step_time(P,dt_dx)
  P_out=zeros(size(P));
  [flux_vals,CFL] = flux(P,dt_dx);
  % assume periodic domain
  P_out = P+circshift(flux_vals,[0,1])-flux_vals;
end

% flux_val(j) is flux across boundary (j-onehalf)
function [flux_val,CFL] = flux(P,dt_dx)
  % define boundary state value to be the
  % average of the values of the states to the
  % left and the right of the boundary.
  P_shifted = circshift(P,[0,-1]); % values to right
  P_boundary = 0.5*(P+P_shifted); % avg of left and right
  %% define boundary wave speeds
  wave_vel = wave_velocity(P_boundary);
  %wave_vel = wave_velocity(P_shifted);
  % compute Courant number
  max_wave_vel=max(wave_vel);
  min_wave_vel=min(wave_vel);
  max_spd=max(abs([max_wave_vel min_wave_vel]));
  %
  % The Courant-Friedrichs-Levy number should not exceed 1
  % (i.e., since our explicit method only looks as far as
  % the neighboring cell for information, the time step
  % we use must be short enough that a density waves should
  % not travel a distance exceeding a cell width in one 
  % time step.)
  %
  % Another approach would be to adjust the time step
  % so that the CFL number is always some fixed value
  % (e.g. 1).
  %
  CFL=max_spd*dt_dx;
  disp(['CFL=' num2str(CFL) ...
        ', min_wave_vel=' num2str(min_wave_vel) ...
        ', max_wave_vel=' num2str(max_wave_vel) ]);
  % calculate the fluxes
  vel=velocity(P);
  flux_rate = P.*vel;
  %
  % up-wind the fluxes
  %
  % by default take flux values from states below
  flux_val=flux_rate*dt_dx;
  % for high traffic densities the density wave travels _backwards_
  % and we need to get the information from upstream.
  % so where wave_vel is negative up-wind fluxes from above
  upwind_indices = find(wave_vel<0);
  flux_val_downshift=circshift(flux_val,[0,-1]);
  flux_val(upwind_indices)=flux_val_downshift(upwind_indices);
end

function output = wave_velocity(P)
  output=get_velocity(P,2);
end

function output = velocity(P)
  output=get_velocity(P,1);
end

% This function assumes that car speed is
% a simple linear decreasing function of density.
%
% type: 2 = density wave speed
%   (the derivative of flux with respect to density)
% type: 1 = car speed
function output = get_velocity(P,type)
  % use a linear velocity function.
  % when P is zero V is Vmax.
  % when P is Pmax, V is zero.
  Vmax=1;
  Pmax=10;
  output = Vmax*(1-P*(type/Pmax));
end
