/////////////////////////////////////////////
// trafic routier: perturbation de l'équilibre  
// calcul de la valeur alpha0
/////////////////////////////////////////
alphaC=10;
alphaV=40;
//////////////////
/////////////////////////////////
N=50;
alphainf=20; Vinf=V*Fval(alphainf);  // 
//
V=30;alphainit=5;dt=0.2;Tfin=50;
///////////////////////////////////
function uval=Fval(alpha)
    uval=0;
    if (alpha>alphaC) then
      uval=1-exp(-(alpha-alphaC)/(alphaV-alphaC))
    end
endfunction
//
function uval=Fprimeval(alpha)
    uval=0;
    if (alpha>alphaC) then
      uval=1/(alphaV-alphaC)*exp(-(alpha-alphaC)/(alphaV-alphaC))
    end
endfunction
//
function uval=phi(alpha)
uval=Fprimeval(alpha)*(-alpha)+Fval(alpha)
endfunction
//
// calcul de alpha0
a=alphaC+0.01
b=20*alphaC
//
while((b-a)>0.001)
 c=(a+b)/2
  if(phi(c)>0) then
     b=c 
  else
     a=c     
   end  
end
//
disp('valeur de alpha0:');disp((a+b)/2)
Kmax=floor(Tfin/dt)
Pos=zeros(Kmax+1,N)
Pos(1,1:N/2)=linspace(0,alphainf*(N/2-1),N/2);
Pos(1,N/2+1:N)=linspace(alphainf*(N/2),0.8*alphainf*(N),N/2);
Pos(:,N)=(linspace(Pos(1,N),Pos(1,N)+Vinf*Tfin,Kmax+1))'
figure(1)
for k=1:Kmax
 drawlater      
 for n=N-1:-1:1
     Pos(k+1,n)=Pos(k,n)+dt*(V*Fval(Pos(k,n+1)-Pos(k,n)))
 end
 clf()
 plot2d(Pos(k+1,:),Tfin/2*ones(1,N),-1,rect=[0,0,N*alphainf+Tfin*Vinf,Tfin]) 
 drawnow
end 
///////////////
figure(2)
clf()
for k=1:Kmax
   plot2d(Pos(k+1,:),k*dt*ones(1,N),-1,rect=[0,0,N*alphainf+Tfin*Vinf,Tfin])       
end
//
