// ODE_07_o2_d4_Trajectoire_avec_vent.sce
// Résolution d'une équation différentielle
// par la méthode de Runge qui est d'ordre 2.
// Un objet est lancé dans un plan  (x, 0 ,z) avec une
// vitesse initiale et subit une force de frottement
// proportionnelle au carré de sa vitesse relative au vent,
// dans la direction opposée à sa vitesse relative au vent.

function z=f(x,y)
//===============
// V' = g - (mu/m) * ||V_rel|| * V_rel
// V_rel = V - V_vent
// y(1) = position  x  de l'objet
// y(2) = vitesse   x  de l'objet
// y(3) = position  y  de l'objet
// y(4) = vitesse   y  de l'objet
// z(j) = y'(j)
// z(1)= = dérivée de la position en x
// z(2)= = dérivée de la vitesse en x
// z(3)= = dérivée de la position en y
// z(4)= = dérivée de la vitesse en y
g  = 9.81;  // [m/s^2]
m  = 5;     // masse de l'objet en [kg]
mu = 0.05;  // coefficient de frottement [N * s^2 / m^2]
Vvent_x = -15; // [m/s]
Vvent_y =  0; // [m/s]
Vrel_x  = y(2) - Vvent_x;
Vrel_y  = y(4) - Vvent_y;
Vrel    = sqrt(Vrel_x^2 + Vrel_y^2);
z(1) = y(2);
z(3) = y(4);
z(2) = -(mu/m) * Vrel * Vrel_x;
z(4) = -9.81 -(mu/m) * Vrel * Vrel_y;
endfunction

function y=Ordre2(x0,xFin,y0,nMax)
//================================
// pour la condition initiale  y(x0) = y0
// résoud l'équation différentielle y'(x)=f(x,y(x))
// pour  x  variant de  x0  à  xFin
// retourne un vecteur de  nMax+1  valeurs de  y.
y    = zeros(length(y0), nMax+1);
y(:,1) = y0;
dx    = (xFin - x0) / nMax;

for j=1:nMax do
  x = x0 + (j-1)*dx
  k1 = f(x, y(:, j));
  k2 = f(x + dx/2, y(:, j) + k1 * dx/2);
  y(:, j+1) = y(:, j) + dx*k2;
  end;
endfunction

// Utilisation de la résolution
// données :
a  = 0;
b  = 4.2;
y0 = [0; 10; 0; 25];
N  = 300;

// résolution numérique
y  = Ordre2(a,b,y0,N);

// graphique de la solution numérique
h = (b-a)/N;
x = a:h:b;

// Animation
//==========
clf;
a=gca();
a.isoview="on";  // axes x et y ont la même échelle
plot(y(1,1:2), y(3,1:2));
ee1 = gce();
ee1.children(1).thickness = 9;
plot([0, 25], [0,25], "w"); // pour définir les dimensions

//realtimeinit(2/N); // définit un intervalle de temps => temps total = ? [s]
// déjà trop lent, donc inutile de ralentir plus.
for j=1:N
//  realtime(j); // attend le temps définit dans realtimeinit(...)
  drawlater();
  ee1.children(1).data = [y(1,j), y(3,j) ; y(1,j)+0.5, y(3,j)];
  drawnow();
  if (y(3,j) < 0) then break; end;
  end;
  
// pour voir la trajectoire 
plot(y(1,:),y(3,:));
