// Résolution numérique d'une équation différentielle
// Navire au mouillage

clear;
//Donnée
Se=3;     //Surface dans l'eau
Sa= 10;   // Surface dans l'air
V=50 *1852/3600;     // Vitesse vent
Pl=2.25 * 9.81;    // Poid lineaire chaine
m=7000;   // masse Navire
h=4;      //Hauteur au dessus du fond

//Constante
K1 =0.6*Sa*V*V/m; // force vent
K2 = Pl/(2*h*m); //chaine
K3 = 500*Se/m; //res eau

// on réécrit l'équa diff sous la forme
// y'' = d - b y' - c y = f(y, t)
// on utilise une représentation d'état (voir espace d'état)
// c'est à dire on remplace l'équa diff du second ordre
// par un système de deux équations différentielles du premier ordre
// on pose u1=y et u2=y'

function [f]=Equa(t,u)
  f(1) = u(2);
  f(2) = K1 - K3*u(2)*abs(u(2)) - K2*u(1)*u(1);
endfunction

t0=0;
t=[0:.1:500];
u0 = [10 ; 0]; // y0=0 et y'0=0   declarer ici postion en x et vitesse a t=0

[U]=ode("rk",u0,t0,t,Equa);

// Affichage
clf();
yp = U(1,:); // extraction ligne 1 pos
yv = U(2,:); // extraction ligne 2 vit
subplot(221);
plot(yp,yv,"r");
xtitle( 'Phase', 'Position', 'Vitesse' ) ;
//scf();
subplot(223);
plot(t,yp);
xtitle( 'Position', 'Temps', 'Position' ) ;
//scf();
subplot(224);
plot(t,yv);
xtitle( 'Vitesse', 'Temps', 'Vitesse' ) ;

//recherche effort maximum
xmax=max(yp);
Tmax=Pl*xmax*xmax/(2*h);
k=Tmax/Pl;
a= 2*atan(tanh(xmax/(2*k)));
Fmax = Tmax / cos(a);
L=k*sinh(xmax/k);
//sortie
subplot(222);
stra=msprintf("longeur de chaine : %.1f m",L) ;
xstring (0,0,stra );
//subplot(313);
stra=msprintf("Effort maximum : %.0f daN",Fmax/10) ;
xstring  (0,.2,stra );

// FM4PN JP Lathuile 
// Droit Xynium codage
