function mit18086_traffic
%MIT18086_TRAFFIC
%    Solves a traffic model equation by a characteristic
%    particle method.

% 02/2007 by Benjamin Seibold
% Feel free to modify for teaching and learning.

n = 400;        % number of space gridpoints without boundaries
dt = 2e-2;                                          % time step
tf = 8e-0;                                         % final time
nc = 80;                                       % number of cars
x = linspace(-12,4,n)';
x0 = x;
h = x(2)-x(1);
u = f(x);
u0 = u;
cu = cumsum(u); cu = cu/cu(end);
c = interp1(cu,x,linspace(0,1,nc));
for tn = 1:ceil(tf/dt)
   x = x+dt*(1-2*u);
   c = c+dt*interp1(x,1-u,c);
   [mv,mi] = min(diff(x));
   while abs(mv)<dt*1.3
      x = [x(1:mi-1);(x(mi)+x(mi+1))/2;x(mi+2:end)];
      u = [u(1:mi-1);(u(mi)+u(mi+1))/2;u(mi+2:end)];
      [mv,mi] = min(diff(x));
   end
   [mv,mi] = max(diff(x));
   while abs(mv)>h*2
      x = [x(1:mi);(x(mi)+x(mi+1))/2;x(mi+1:end)];
      u = [u(1:mi);(u(mi)+u(mi+1))/2;u(mi+1:end)];
      [mv,mi] = max(diff(x));
   end
   clf
   plot(x0,u0,'b:',x,u,'r-')
   hold on
   plot(c,c*0,'s','LineWidth',1,'MarkerEdgeColor','k',...
                  'MarkerFaceColor','g','MarkerSize',3)
   text(0,0,'cars','VerticalAlignment','top')
   hold off
   axis([-4 4 -.1 1.1])
   title(sprintf('time t=%0.2f',tn*dt))
   drawnow
end

function y = f(x)
% initial condition function
y = .5./(1+4*(x+3).^2)+.8./(1+5*(x+1).^4)+.3./(1+2*(x-1).^2)+.1;