// ODE_11_o4_d6_Le_pendule_3D.sce
// Résolution d'une équation différentielle
// par la méthode de Runge-Kutta qui est d'ordre 4.
// Un pendule attaché à un point fixe oscille.
// Il n'est soumis qu'à la gravitation et à 
// la tension dans le fil.

funcprot(0); // pour ne pas indiquer une redéfinition des fonctions

global glnCount; // compte le nombre d'appel à la fonction
glnCount = 0;

function z=f(x,y)
//===============
// pour résoudre vectoriellement :
// v' = g + (g*z - v^2)/L^2 * (x,y,z)
// L = la longueur du fil.
// L'origine est le point d'attache du fil.
// La coordonnée  z  est donc négatif pour de petites oscillations.
// y(1) = position x
// y(2) = vitesse  x
// y(3) = position y
// y(4) = vitesse  y
// y(5) = position z
// y(6) = vitesse  z
// z(j) = y'(j)
g = 9.81;   // en m/s^2
L = 0.5;    // longueur du fil en [m]
Coef = (g * y(5) - y(2)^2 - y(4)^2 - y(6)^2) / L^2;  // < 0 en général.
z(1) = y(2);
z(3) = y(4);
z(5) = y(6);
z(2) =      Coef * y(1);
z(4) =      Coef * y(3);
z(6) = -g + Coef * y(5);
global glnCount;
glnCount = glnCount + 1; // compte le nombre d'appel à la fonction
endfunction

glnCount = 0; // compte le nombre d'appel à la fonction
// Utilisation de la résolution
// données :
t0   = 0;
tFin = 60;
L    = 0.5;  // longueur du fil en mètres 
px0 = -0.2;
py0 =  0.00;
pz0 = -sqrt(L^2 - px0^2 - py0^2);
y0 = [px0; 0.0;.. // pos_x., vitesse_x
      py0; 0.1;.. // pos_y., vitesse_y
      pz0; 0.00];  // pos_z., vitesse_z
N  = 1000;

// résolution numérique
   ax = linspace(t0, tFin, N+1);  // vecteur LIGNE
     
   // résolution numérique en appelant la fonction  ode  de SciLab
   ay = ode( y0,..     // condition initiale
             t0,..     // x_initial,
             ax,..     // série des x(j) pour lesquels  y(j)  sera calculé 
             1e-4,..   // tolérance d'erreur relative
             1e-5,..   // tolérance d'erreur absolue
             f..       // fonction définissant l'équation différentielle
             );
   // L'algorithme essaye d'avoir une erreur inférieure au
   // maximum entre l'erreur absolue et l'erreur relative.
 
// graphique de la solution numérique
//===================================
clf;
plot(ay(1,1:2), ay(3,1:2));
ee1 = gce();
ee1.children(1).thickness = 9;
aaa = gca();
aaa.data_bounds=[-0.3, -0.3 ; 0.3, 0.3];
aaa.isoview="on";  // axes x et y ont la même échelle
aaa.margins = [0.01, 0.01, 0.01, 0.01];  // pour agrandir l'affichage.
aaa.box = "on"; // pour tracer une boîte autour du graphique.
// lorsque la fenêtre graphique existe,
// tapez  aaa  dans la console Scilab pour plus d'options.

plot(ay(1,:), ay(3,:));  // trajectoire dans le plan X-Y

disp("glnCount = " + string(glnCount));

return;

// Animation du corps
for j=1:1:100
  drawlater();
  ee1.children(1).data = [ay(1,j), ay(3,j) ; ay(1,j)+0.01, ay(3,j)];
  drawnow();
  end;

