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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 dt=.01;
 t=0;
N=400;    % number of cars
L=15;     % length of a car in feet--density  p=1/l5
v1=100;   % max vel for nondense traffic in fps
v2=.4*v1; % max vel for dense traffic
gg=L+60*(0:1:N)/N;
R=5280./gg;
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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r=1./g;

if k==1
 mmax=24000
else
 mmax=12000
end
for m=1:mmax

RR=(u<vs2)/s2+(1/s2+(1/s1-1/s2)*(u-vs2)/(vs1-vs2)).*(vs2<=u).*(u<=vs1)+(u>vs1)/s1;

x=x+dt*u;
a=a-dt*a/e+dt*(v2-v1)*(1-L*r).*(r>=RR)/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
y=[x-(q+1)*p,x-q*p]/5280;
uu=[u,u]*3600/5280;
rr=[r,r]*5280;

   ff1=sprintf('speed mph vs x in miles\nminutes  = %.2f', t/60);
   subplot(131)
   plot(y,uu,'r ',y(N+1),uu(N+1),'k +')
   axis([0, p/5280, 0,max(u*3600/5280)+2])
   title(ff1)
   subplot(132)
   plot(y,rr,'r ',y(N+1),rr(N+1),'k +')
   axis([0,p/5280,R(1)/2-1,R(1)+1])
   ff2=sprintf(' density in cpm vs x in miles\n[e, k] = %d, %d', e, k);
   title(ff2)
%%new code----------------------------------------------------------------------
umin=min(u);
umax=max(u);
zmin=(u==umin);
zmax=(u==umax);
gmin=min(g(zmin));
gmax=max(g(zmax));
rmax=1/gmin;
rmin=1/gmax;
s=(rmin-rmax)/(umax-umin);
ss=(1/s2-1/s1)/(vs2-vs1);
Uo=(-1/s2+1/gmin-s*umin+ss*vs2)/(ss-s)-2;
Ro=1/s2+ss*(Uo-vs2);
%-----------------------------------------------------------------
u1=3600*([u,u(1)]-Uo)/5280;
r1=5280*([r,r(1)]-Ro);
D=(u1.^2+r1.^2).^.5;

cc1=u1./D;
ss1=r1./D;
Uo=Uo*3600/5280;
Ro=Ro*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;

[-s,   min(D)]
%% end new code-----------------------------------------------------------------------------------------------------------------------
subplot(133)
plot(V1,R,'b',V2,R,'b',U,R1,'r',[u,u(1)]*3600/5280,5280*[r,r(1)],'k',u(1)*3600/5280,5280*r(1),'r o')
axis([0 45,R(1)/2-1,R(1)+1])
title(' car density in cpm vs velocity in mph')

 pause(0)
 else
 end
end
pause
end
pause
end