%2-phase traffic flow 
for n=1:4
e=2^(n-1);
for k=1:3

dt=.01; 
t=0;

N=400; 
L=15;       % length of a car in feet--density  p=1/l5
v1=100;     % max velocity in less dense traffic
v2=.4*v1;   % max vel for dense traffic

gg=L+60*(0:1:N)/N;
V1=3600*v1*(1-L./gg)/5280;
V2=3600*v2*(1-L./gg)/5280;
U=(0:v1/N:v1);
s2=4*L/3;
vs2=v2*(1-L/s2);
s1=s2;
vs1=v1*(1-L/s1);
l=.5*(s1+s2);
p=N*l;
y=l*(1:1:N);
x=y+sin(2*k*pi*y/p);
g=[diff(x),x(1)+p-x(N)];
 
Uo=.5*(vs2+vs1);
u=Uo*ones(size(x));
a=u-v1*(1-L./g);
step=60;

R1=(1/s2)*ones(size(U));

U=3600*U/5280;
R1=5280*R1;
th1=(1:1:N)*2*pi/N;
C=cos(th1);
S=sin(th1);

%update
if k==1
 mmax=24000;
else
 mmax=12000;
end

for m=1:mmax

  RR=(1/s2)*ones(size(u));
  G=1./RR;
  x=x+dt*u;
  a=a*(1-dt/e)+dt*(v2-v1)*(1-L./g).*(g<=G)/e;

  g=[diff(x),x(1)+p-x(N)];
  r=1./g;
  v=v1*(1-L*r);
  u=v+a;
  t=t+dt;

  j=m-step*floor(m/step);
  if j==0
    q=(x(1)-mod(x(1),p))/p; %%integer

  yy=1:1:N;
  uu=u*3600/5280;
  rr=r*5280;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
umin=min(u);
umax=max(u);
zmin=(u==umin);
zmax=(u==umax);
gmin=min(g(zmin));
gmax=max(g(zmax));
s=(gmax-gmin)/(umax-umin);
ss=(s2-s1)/(vs2-vs1);
Go=s2;
Uo=(s2-gmin+s*umin-ss*vs2)/(s-ss);
g1=[g,g(1)]-Go;
u1=[u,u(1)];
u1=(u1-Uo);
D=(u1.^2+g1.^2).^.5;

ss1=u1./D;
cc1=g1./D;
Uo=Uo*3600/5280;
th=asin(ss1);
th=(ss1>=0).*(th.*(cc1>=0)+(pi-th).*(cc1<0))+(ss1<0).*((pi-th).*(cc1<0)+(2*pi+th).*(cc1>=0));

phi=diff(th);
phi=phi.*(phi>-pi).*(phi<pi)+(phi+2*pi).*(phi<=-pi)+(phi-2*pi).*(phi>=pi);

w=-sum(phi)/2/pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   subplot(131)
   plot(yy,uu,'r ')
   axis([1,N , 0,floor(max(u*3600/5280))+2])
   ff1=sprintf('speed in mph vs car\nlabel: minutes = %.2f', t/60);
   title(ff1)

   subplot(132)
   plot(yy,g,'r')
   axis([1,N,15,min([100,floor(max(g))+2])])
   ff2=sprintf('car spacing in feet vs car\nlabel: [e,  k] = %d %d', e, k);
   title(ff2)

   subplot(133)
   plot(gg,V1,'b',gg,V2,'b',5280./R1,U,'r',[g,g(1)],[u,u(1)]*3600/5280,'k',g(1),u(1)*3600/5280,'r o')
   axis([L,min([100,max(g)+1]),0, min([100,max(3600*u/5280)+2])])
   title('velocity in mph vs car spacing in feet')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pause(0)
 else
 end
end
pause
end
pause
end