
//////////////////////////////////////////////////////
//               Portrait de phase                  //
//////////////////////////////////////////////////////

///////////////////////////////////////////////////////
// La fonction portrait_phase definie dans ce programme
// a l'action suivante:
// - Elle trace le champ de vecteur defini par f
// (f etant une fonction a definir dans le programme
// par defaut fixer a celle du systeme de Volterra)
// - Quand on clique gauche sur un point du graphe
// affiche, elle trace la trajectoire partant de ce point.
// Cela continue tant qu'on clique gauche. Pour sortir 
// du programme, un LONG clic droit sur le graphique
// est necessaire.

///////////////////////////
//  Parametres

a=2

/////////////////////////////////////
// La fonction f a modifier sur demande.
// le temps est necessaire dans la syntaxe de 
// fchamp et de ode. Le mettre meme si comme pour
// Volterra le systeme est autonome. 

function Y =f(t,X)
// Mettre ici la definition de votre fonction
// Attention X et Y doivent etre des vecteur 
// 2 lignes et 1 colonnes.
// Pour l'instant, la fonction utilisee est 
// celle de Volterra
Y = [ X(1)*(1-X(2)); a*X(2)*(X(1)-1)]
endfunction

// Pour tracer les trajectoires vers les x negatifs.

function Y = fneg(t,X)
  Y = - f(t,X)
endfunction


/////////////////////////////////////////
// Une fonction utile pour ne pas tracer 
// de fleches trop grandes
function Y=fnorm(t,X)
  Y= f(t,X)
  if norm(Y)> .2 then Y= Y/norm(Y)
  else Y = Y/.2
  end
endfunction

/////////////////////////////////////////////////
//////////////////////////////////////////////////
// Une fonction qui trace le portrait 
// de phase du systeme y'=f(Y), pour la fonction 
// Y definie plus haut.
// xmin,ymin: coordonnees du point bas gauche du graphique
// xmax,ymax: coordonnees du point haut droit du graphique

function []=portrait_phase(xmin,xmax,ymin,ymax)

// Creation des vecteurs abscisse et ordonnees.
  dx= (xmax - xmin)/20 ; dy = (ymax- ymin)/20
  Xp=xmin:dx:xmax; Yp =ymin:dy:ymax
  rect1=[xmin,ymin,xmax,ymax] 
  clf()  //on efface les graphiques restants

// Trace du champ de vecteurs
  f0=scf(0)
  fchamp(fnorm,0,Xp,Yp)
   // f0.children.children($).colored="on"

// Trace des trajectoires.
  tmax=50; dt= 0.1
  T=0:dt:tmax

  while %t do 
    [c,x,y]=xclick()
    if (c==5)|(c==2)|(c==-1000) then break
    else if (x<xmax)&(x>xmin)&(y<ymax)&(y>ymin) then
           X0=[x;y]
           Xsol = ode(X0,0,T,f)
           plot2d(Xsol(1,:)',Xsol(2,:)',rect=rect1)
           f0.children.children(1).children.thickness=2
           f0.children.children(1).children.foreground=5          
              // pour augmenter la largeur des trajectoires.
          Xsol2 = ode(X0,0,T,fneg)
           plot2d(Xsol2(1,:)',Xsol2(2,:)',rect=rect1)
           f0.children.children(1).children.thickness=2
           f0.children.children(1).children.foreground=5  
         end
    end
  end
endfunction
