%!
/irnd 1.0 2147483647.0 div def
/myrand {rand irnd mul} bind def
/accel 0.66 def
/zweib 2.00 def
/L 450 def
/miny 100 def
/maxy 700 def
/nrcars 70 def
/nrcars1 nrcars 1 sub def
/nrcars2 nrcars1 1 sub def
/nrcars3 nrcars2 1 sub def
/vmax 5 def
/xoff 50 def
%
% the color map, so to speak
%
/rcol vmax 1 add array def
/gcol vmax 1 add array def
/bcol vmax 1 add array def
rcol 0 0.0 put  % v = 0 red 
gcol 0 0.0 put  % v = 0 green 
bcol 0 0.0 put  % v = 0 blue
rcol 1 1.0 put  % v = 1 red 
gcol 1 0.0 put  % v = 1 green 
bcol 1 0.0 put  % v = 1 blue
rcol 2 1.0 put  % v = 2 red 
gcol 2 1.0 put  % v = 2 green 
bcol 2 0.0 put  % v = 2 blue
rcol 3 1.0 put  % v = 3 red 
gcol 3 1.0 put  % v = 3 green 
bcol 3 0.0 put  % v = 3 blue
rcol 4 0.5 put  % v = 4 red 
gcol 4 1.0 put  % v = 4 green 
bcol 4 0.0 put  % v = 4 blue
rcol 5 0.0 put  % v = 5 red 
gcol 5 1.0 put  % v = 5 green 
bcol 5 0.0 put  % v = 5 blue
/Helvetica findfont
20 scalefont
setfont
newpath
xoff 10 sub maxy 5 add moveto
xoff L add maxy 5 add lineto
stroke
xoff L add 0.5 mul maxy 28 add moveto
(space) show
%
newpath
xoff 5 sub miny moveto
xoff 5 sub maxy 10 add lineto
stroke
xoff 28 sub miny maxy add 0.5 mul moveto
90 rotate
(time) show
-90 rotate
xoff maxy 50 add moveto
/rho 9 string def
/Symbol findfont
14 scalefont
setfont
nrcars L div rho cvs
(r = ) show
/Helvetica findfont
14 scalefont
setfont
rho show
/N 5 string def
nrcars N cvs
(, N = ) show
N show
(,  vmax = 5,  pbrake = 0.5, ruleset: smooth braking) show
pop 
pop
/vloc nrcars array def
/vnew nrcars array def
/xold nrcars array def
/xnew nrcars array def
/initSystem 
{
  0 1 nrcars1 
  {
    dup dup
    xold exch dup L nrcars div mul put      	% xold[i] =3*i
    vloc exch 5 put        	% vloc[i] = 0
    xnew exch 0 put		% xnew[i] = 0
  } for
    vloc nrcars1 0 put
    vloc nrcars2 0 put
    vloc nrcars3 0 put
} def 
%
%
%
/initSystem1 
{
  0 1 nrcars1 {
    /indx exch def
    xold indx indx 2 mul put      	% xold[i] = 2*i
    vloc indx 1 put        	% vloc[i] = 1
    xnew indx 0 put		% xnew[i] = 0
  } for
} def 
%
%
%
% mod for floats: 
%
/fmod
{
  /divisor exch def
  /dividend exch def
  divisor dividend le {/dividend dividend divisor sub def} if
  dividend 
} bind def
%
/updatePos 
{
  0 1 nrcars1 {
    dup
    dup
    xold
    exch
    get
    exch
    vloc
    exch
    get
    add
    L
    fmod
    xnew
    3 1 roll
    put
  } for
} bind def
%
% next procedure: prepareTimeStep; assures xold[i] = xnew[i]
%
/prepareNextStep 
{
  0 1 nrcars1 {
    dup
    xnew
    exch
    get
    xold
    3 1 roll
    put
  } for
} bind def
%
%
%
/vnewtovloc 
{
  0 1 nrcars1 {
    dup
    vnew
    exch
    get
    vloc 
    3 1 roll
    put
  } for
} bind def
%
% procedure drawLines: draws space-time lines in specified color
%
/drawLines 
{
  /y exch def   	% local variable y
  0 1 nrcars1 {
    /indx exch def
    /x2 xnew indx get xoff add def 	% local variable x2
    /x1 xold indx get xoff add def  	% local variable x1
    x1 x2 le {
      newpath
      x1 y 1 add moveto
      /v vloc indx get cvi def    % color...
      v 3.0 gt {/rgbind 4  def} if
      v 3.0 le {/rgbind v 1 idiv def} if
      rcol rgbind  get
      gcol rgbind  get
      bcol rgbind  get
      setrgbcolor
      % vloc indx get dup vmax div dup mul 0.4 mul 
      % exch 0 gt {1}{0}ifelse 1 exch sethsbcolor
      x2 y lineto
      stroke
    } if
  } for
} bind def
%
% next procedure is calcNewV: applies Stefan Krauss's rules to every car
%
/calcNewV 
{
  0 1 nrcars2 
  {
    /indx exch def
    /v vloc indx get def
    /thisx xold indx get def
    /nextx xold indx 1 add get def
    /nextv vloc indx 1 add get def
    /gap nextx thisx sub 1 sub def
    gap 0 lt {/gap gap L add def} if 	% boundary condition
  %  /v2 v 1 add v mul 0.5 mul def                   % v2 = v(v+1)
  %  /v2p nextv 1 sub nextv mul gap add 0.5 mul def  % v2p = v'(v'-1) + gap
  %  v2 v2p le
  %  {
  %    v2 v2p lt
  %    {
%	v vmax lt {v 1 add /v exch def} if
 %     } if
  %    v 0 gt { /v v myrand round cvi sub def } if
   % }
   % {
   %   v 0 gt {/v v 1 sub def} if
   % } ifelse
    /v v accel add def
    v vmax gt {/v vmax def}if   
    /vsafe 
        gap nextv sub
        v nextv add zweib  add
        div
        zweib mul
        nextv add
    def
    v vsafe gt {/v vsafe def}if
   /sigma myrand accel mul def
   /v v sigma  sub def
    v 0 lt {/v 0 def}if
  
    vnew indx v put
  } for		% main loop over all the cars
%
% last car is somewhat special
%
  /v vloc nrcars1 get def
  /thisx xold nrcars1 get def
  /nextx xold 0 get def
  /nextv vloc 0 get def
  /gap nextx thisx sub 1 sub def
  gap 0 lt {gap L add /gap exch def} if

 % /v2 v 1 add v mul 0.5 mul def                   % v2 = v(v+1)
 % /v2p nextv 1 sub nextv mul gap add 0.5 mul def  % v2p = v'(v'-1) + gap
 % v2 v2p le
 % {
 %   v2 v2p lt
 %   {
 %     v vmax lt {v 1 add /v exch def} if
 %   } if
 %   v 0 gt { /v v myrand round cvi sub def } if
 % }
 % {
 %   v 0 gt {/v v 1 sub def} if
 % } ifelse
 
 /v v accel add def
    v vmax gt {/v vmax def}if   
    /vsafe 
        gap nextv sub
        v nextv add zweib  add
        div
        zweib mul
        nextv add
    def
  v vsafe gt {/v vsafe def}if
 /sigma myrand accel mul def
 /v v sigma sub def
 %v  pstack
  v 0 lt {/v 0 def}if
  vnew nrcars1 v put
  vnewtovloc           % required for parallel update...
} bind def
%
% Initialization
%
usertime
srand
initSystem
1 1 1 
{			% dummy iteration, without drawing
  pop			% pop iteration index from stack
  calcNewV
  updatePos
  prepareNextStep
} for
%
% do the simulation...
%
maxy -1 miny 
{
  calcNewV
  updatePos
  drawLines		% assumes time on stack
  prepareNextStep
} for
showpage
%%EOF