//
//  programme 9
//
write(%io(2),'prog 9: modele proies-predateurs de Dunbar ([Mur], page 315)')
//
// parametre du calcul 
//
Mu=0.6;val=(1+Mu)/2;Beta=0.4;dt=0.3;
N=100;Longueur=50;dx=Longueur/(N-1);
LD=dx.*(0:N-1);
//
U=[];V=[];c=[];A=[];
//
// conditions initiales
//
for I=1:N/2,
  U(I)=0;V(I)=1;
  U(I+N/2)=1-Mu;V(I+N/2)=Mu;
end
//
// matrice du laplpacien discretise
//
for I=2:N-1,
  A(I,I)=-2;A(I-1,I)=1;A(I+1,I)=1;
end
A(1,1)=-1;A(2,1)=1;
A(N,N)=-1;
A(N-1,N)=1;
//
AZ=eye(N,N)-A*(dt/(dx^ 2));
B=inv(AZ);
//
for J=1:300		// debut boucle en temps
//
// etape Euler
//
  for I=1:N
    U(I)=U(I)+dt*Beta*U(I)*(V(I)-Mu);
    V(I)=V(I)+dt*V(I)*(1-V(I)-U(I));
  end
//
// etape diffusion implicite
//
U=B*U;
//
// calcul de la vitesse du front
//
for I=2:N-1,
  Vpr(I)=(V(I+1)-V(I))/dx;
  c(I)=V(I)*(1-V(I)-U(I))/(Vpr(I)-0.0000001);
end
  c(1)=c(2);c(N)=c(N-1);
//
//   recherche du point X tq V(X)=val
//
Indice=N/2;
for I=1:N
 if V(I)>val then Indice=I;end
end;
//
// Translation du domaine de calcul
// (convergence supposee dans la partie rajoutee)
//
 Translation=-Indice+N/2;
 Translat=abs(Translation);
 if Translation <0 then
    for I=1:N-Translat
      V(I)=V(I+Translat);
      U(I)=U(I+Translat);
    end
    for I=N-Translat:N,
      V(I)=Mu;U(I)=1-Mu;
    end
  end
  if Translation >0  then
    for I=N:-1:1+Translation
      V(I)=V(I-Translat);
      U(I)=U(I-Translation);
    end
    for I=1:Translation
      V(I)=1;U(I)=0;
    end
  end
// 
// trace des resultats a l'instant t
//
  if or([int(J/20)==J/20,J==1]) then
   xbasc()
   xsetech([0,0,0.5,0.5])
   str=string(J*dt);
   xstring(0.5,0.5,['t=',str])
   plot2d(LD',U,1,'011',' ',[0,0,Longueur,1])
   xtitle('predateurs dans le repere mobile','z','u(z)')
   xsetech([0.5,0,0.5,0.5])
   xtitle('proies dans le repere mobile','z','v(z)')
   plot2d(LD',V,1,'011',' ',[0,0,Longueur,1])
   xsetech([0,0.5,1,0.5])
   xtitle('vitesse de deplacement','z','c(z)')
   plot2d(LD',c,1,'011',' ',[0,0,Longueur,1])
  end
//
// fin de la boucle en temps
//
end
