// **********
// Initialisation
// **********

clear;

// **********
// Données
// **********

X0 = [1, 2, 5, 7, 9, 6, 3, 8];
Y0 = [7, 6, 8, 7, 5, 7, 2, 4];

// **********
// Fonctions
// **********

function [a] = regression_quadratique(X, Y) 
    // méthode de la distance algébrique
    // X, Y : points expérimentaux, matrices colonnes réelles
    // a : coefficients de la formule quadratique (matrice colonne réelle)
    YY = (X.*X); 
    XX = [(X.*Y) ; (Y.*Y) ; X ; Y]; 
    [aa, b] = reglin(XX, YY);
    a = [-aa, -b];
endfunction

function [phi]=trouve_rotation(A)
    // A : coefficients de la formule quadratique (matrice colonne réelle)
    // phi : angle que fait un axe de l'ellipse avec x (radians)
    delta = 1 - 1/(1 + (A(3) - A(1))^2/A(2)^2);
    absphi = acos(sqrt((1 + sqrt(delta))/2));
    cap = cos(absphi); sap = sin(absphi);
    signephi = sign(A(2)*(cap*cap - sap*sap)/(A(1) - A(3)));
    phi = signephi*absphi;
endfunction

function [x,y]=trouve_centre(A)
    // A : coefficients de la formule quadratique (matrice colonne réelle)
    // x, y : coordonées du centre de l'ellipse (réels)
    delta = A(2)*A(2) - 4*A(1)*A(3);
    x = (2*A(3)*A(4) - A(2)*A(5))/delta;
    y = (2*A(1)*A(5) - A(2)*A(4))/delta;
endfunction

function [rx, ry]=trouve_rayons(a, phi, xc, yc)
    // a : coefficients de la formule quadratique (matrice colonne réelle)
    // phi : angle que fait un axe de l'ellipse avec x
    // xc, yc : coordonnées du centre de l'ellipse
    // rx, ry : rayons (grand et petit demi-grands axes) de l'ellipse
    A = [a(1), a(2)/2 ; a(2)/2, a(3)];
    Q = rotate([1, 0 ; 0, 1], phi); // matrice de rotation
    t = [xc ; yc]; // matrice de translation
    Abar = Q'*A*Q;
    b = [a(4) ; a(5)];
    bbar = (2*t'*A + b')*Q;
    c = a(6);
    cbar = t'*A*t + b'*t + c;
    rx = sqrt(abs(cbar/Abar(1,1)));
    ry = sqrt(abs(cbar/Abar(2,2)));
endfunction

function [d]=conique(a, x, y)
    d = a(1)*x'*x + a(2)*x'*y + a(3)*y'*y + a(4)*x + a(5)*y + a(6);
endfunction

function [] = trace_conique(a)
    epsilon = 1.5e-1;
    Xbase = 0:0.1:10;
    Ybase = 1:0.1:9;
    i = 1;
    for x = Xbase
        for y = Ybase
            if abs(conique(a, x, y)) <= epsilon then
                X(i) = x;
                Y(i) = y;
                i = i + 1;
            end
        end
    end
    Xmoy = mean(X); Ymoy = mean(Y); // détermination d'un point intérieur
    theta = atan(Y-Ymoy, X-Xmoy); // angle des vecteurs
    posMat = gsort([theta, X, Y], "lr"); // classement dans le sens trigo
    [m, k] = max(abs(diff(posMat(:, 1)))); // saut de point
    X = [posMat(k+1:$, 2) ; posMat(1:k, 2)]; // décalage circulaire
    Y = [posMat(k+1:$, 3) ; posMat(1:k, 3)];
    plot2d(X, Y, style = 5);
endfunction

// **********
// Programme principal
// **********

// lecture des données

Xdef = X0;
Ydef = Y0;

// Régression
Aopt = regression_quadratique(Xdef, Ydef);
aopt = [1, Aopt];

// affichage des paramètres
disp(aopt)

phi = trouve_rotation(aopt);
phideg = phi*180/%pi;
[xc, yc] = trouve_centre(aopt);
[a, b] = trouve_rayons(aopt, phi, xc, yc);
disp("phi = "+string(phi)+" rad = "+string(phideg)+"°.");
disp("C("+string(xc)+" ; "+string(yc)+").");
disp("a = "+string(a)+" ; b = "+string(b)+".");

// tracé
clf;

plot(Xdef, Ydef, "b+")
replot([0, 10 ; 1, 9]);
plot(xc, yc, "r+")
trace_conique(aopt)