function do_all=traffic
% Implementing individual traffic model from Mesterton-Gibbons, eqn
% 2.111 on page 80.
% Eli Tziperman, for APM 115.

close all; clear all

% initialize:
% length of time of simulation:
Tfinal=150;
% time step for numerical integration.  make sure results don't change
% it it is halved:
dt=0.01;
Nsteps=round(Tfinal/dt);
% number of cars:
Ncars=2;
% initialize variable with 
phi=zeros(Ncars,Nsteps);
tplot=zeros(Nsteps,1);
% delay due to driver response time
tau=2; % stanbdard value: 2.75  
ntau=round(tau/dt); 
% first car velocity dropped from v to v*(1-delta)
delta=0.1;  
% initial time step:
nstart=round(ntau)+1;
% velocity in feet per second (20Miles/hour):
v=88/3;
% density at which traffic stops: 65 cars/mile (5280ft/mile)
xmax=65/5280;

% integrate equation in time:
for n=nstart:Nsteps
  t=(n-nstart)*dt;
  tplot(n)=t;
  % first car us specified to slow down for one second:
  phi(1,n)=phi1(t,v,delta);
  % calculate the response of the other cars:
  for jj=2:Ncars
    phi(jj,n)=phi(jj,n-1)+dt*(88/3)*log(1+xmax*(phi(jj-1,n-ntau)-phi(jj,n-ntau)));
  end
end

% integrate taylor approx to delay eqn:
ode_options=odeset('RelTol',1.e-7,'abstol',1.e-8);
tspan=[tau:tau/10:Tfinal];
X0=[0 0]';
[time,X] = ode45(@(t,X) rhs(t,X,v,tau,xmax,delta),tspan,X0,ode_options);

% plot perturbation displacement for all cars, shift by 65 distance
% units per car to make plot more legible:
for jj=1:Ncars
  plot(tplot,phi(jj,:)+65*jj);
  hold on
end
plot(time,X(:,1)+65*2,'r');
xlabel('time')
ylabel('\phi_j(t)')
title('perturbation car displacement; red: Taylor approx')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=phi1(t,v,delta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p=v*delta*((t+1)*exp(-t)-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rr=rhs(t,X,v,tau,xmax,delta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi2=X(1);
xi=X(2);
rr=[xi; ...
    phi1(t-tau,v,delta)/tau^2-phi2/tau^2+xi*(1/tau-1/(v*xmax*tau^2))];
