// Bille.Sce
clear all
function y = phi(M)
    x = M(1); y = M(2);
    y = 0.5*sqrt(x^4 + y^4);
endfunction

function g = gradphi(M)
    x = M(1); y = M(2);
    g = zeros(2,1);
    g(1) = x^3 / sqrt(x^4 + y^4);
    g(2) = y^3 / sqrt(x^4 + y^4);
endfunction

function zdot = Fbille(t, z)
    m = 1;
    // z=(M,V)=(x,y,vx,vy) vecteur d'etat
    M = z(1:2); 
    V = z(3:4);
    zdot = zeros(4,1)
    zdot(1:2,1) = V;
    zdot(3:4,1) = - gradphi(M) / m;
endfunction

function [td, z] = EulerExpl(t0, z0, F, T, N)
    // N : nb de points
    // t0 : instant initial
    // z0 : données initiale
    // T : temps final
    // F : fonction F(t, z)
    h = (T-t0) / N;
    d = length(z0)
    td =  t0 : h : T; // ! td(1)=t0, td(N+1)=T
    z = zeros(d, N+1);
    z(:,1) = z0;
    for n = 1:N
        z(:,n+1) = z(:,n) + h*F(td(n), z(:,n))
    end; // end for
endfunction

function [td, z] = EulerCauchy(t0, z0, F, T, N)
    // N : nb de points
    // t0 : instant initial
    // z0 : données initiale
    // T : temps final
    // F : fonction F(t, z)
    //
    // zhat = zn + h * F(tn, zn)
    // znew = zn + h * [F(tn, zn) + F(t_{n+1}, zhat)]
    h = (T-t0) / N;
    d = length(z0)
    td =  t0 : h : T; // ! td(1)=t0, td(N+1)=T
    z = zeros(d, N+1);
    z(:,1) = z0;
    for n = 1:N
      zhat     = z(:,n) + h*F(td(n), z(:,n))
      z(:,n+1) = z(:,n) + 0.5*h*( F(td(n), z(:,n)) + F(td(n+1), zhat)  );
    end; // end for
endfunction

// Initialisation
M0 = [2; 1]; V0 = [1;-1]; z0 = [M0; V0];
t0 = 0; T = 40;
N = 10000;
// Resolution du système EDO
//[td, z] = EulerExpl(t0, z0, Fbille, T, N)
[td, z] = EulerCauchy(t0, z0, Fbille, T, N)
//
xt = z(1,:);
yt = z(2,:);
vxt = z(3,:); // !!
vyt = z(4,:); // !!
//figure(1)
//subplot(2,1,1), plot(td, xt, '.-', );
//subplot(2,1,2), plot(td, yt, '.-', );
figure(2)
plot(xt, yt, '.-b'); xgrid;
//figure(3)
//comet(xt, yt);
figure(4)
zt = 0.5*sqrt(xt.^4 + yt.^4);
param3d(xt, yt, zt, 'o-');
mtlb_axis("equal"); 

// Calcul de l'energie totale
m = 1
Etot = 0.5*m*(vxt.^2+vyt.^2) ...
     + 0.5*sqrt(xt.^4 + yt.^4);
figure(5)
plot(td, Etot, 'o-k'); 
xgrid; xtitle('NRJ totale')