function [W] = produit_vectoriel(U, V)
    // But : calcule W = U^V
    // Entrées : U, V, vecteurs ligne
    // Sortie : W, vecteur ligne
    W(1, 1) = U(2)*V(3) - U(3)*V(2);
    W(1, 2) = U(3)*V(1) - U(1)*V(3);
    W(1, 3) = U(1)*V(2) - U(2)*V(1);
endfunction

function [Xp, Yp, B, N, ux, uy] = regression_plan(X, Y, Z)
    // But : se placer dans le plan moyen
    // Moyen : à partir des données expérimentales X, Y et Z
    // déterminer l'équation du plan y = A(1)x + A(2)z + b
    // (soit A(1)x - y + A(2)z + b = 0)
    // par régression multilinéaire
    // puis projetter sur les projections des axes x et y sur le plan
    // Entrées : X, Y, Z : vecteurs colonne de réels
    // Sorties : Xp, Yp : vecteurs colonne de réels
    // N : vecteur normal, vecteur ligne de réels
    // ux, uy, BON du plan, vecteurs ligne
    
    // Régression plane
    // y = A(1)*x + A(2)*z + b
    [A, b, sigma] = reglin([X' ; Z'], Y');
    disp("Équation du plan moyen : ");
    disp("y = ("+string(A(1))+")*x + ("+string(A(2))+")*z + ("+string(b)+")");
    
    // Intersection du plan avec l'axe y
    B = [0, b, 0];
    
    // vecteur normal unitaire
    N = [A(1), -1, A(2)]/norm([A(1), -1, A(2)]);
    
    // vecteur unitaire de l'intersection du plan avec la plan (x, y)
    // soit le plan z = 0
    // droite d'équation A(1)x - y + b = 0, z = 0
    ux = [1, A(1), 0]/norm([1, A(1) 0]);
    
    // deuxième vecteur de la base du plan (produit vectoriel)
    uuy = produit_vectoriel(N, ux);
    uy = uuy/norm(uuy);
    // projection
    Xp = [X, Y, Z]*ux';
    Yp = [X, Y, Z]*uy';    
endfunction

function [xc, yc, r] = rc_kasa_coope(X, Y)
    // calcule les paramètres du cercle décrivant le mieux les points
    // selon la méthode linéaire de I. Kasa et Ian D. Coope
    // Entrée : X, Y, vecteurs colonne de réels
    // Sortie : xc, yc, réels, coordonnées du centre dans le plan
    // r : rayon du cercle
    
    // d = A(1)x + A(2)y + b avec d = x^2 + y^2
    d = X.^2 + Y.^2;
    [A, b, sigma] = reglin([X' ; Y'], d');
    
    // équation cartésienne sous la forme
    // (x^2 + y^2) - (r^2 - xc^2 - yc^2) - (2xc)x - (2yc)y
    xc = A(1)/2; yc = A(2)/2;
    r = sqrt(b + xc^2 + yc^2);
endfunction

function [S] = residus(A, X, Y)
    // But : calculer les mondres carrés totaux
    // si = sqrt((xi - xc)^2 + (yi - yc)^2)
    // S = 1/n*somme si
    // Entrées : A, paramètres de la loi, vecteur ligne de réels
    // A(1) = xc, A(2) = yc, A(3) = r
    // X, Y, vecteurs colonne de réels
    // Sortie : S, réel
    n = size(X, '*');
    s = (((X - A(1)).^2 + (Y - A(2)).^2).^0.5 - A(3)).^2;
    S = sqrt(sum(s)/n);
endfunction

function [xc, yc, r] = regression_circulaire(X, Y)
    // calcule les paramètres du cercle décrivant le mieux les points
    // Entrée : X, Y, vecteurs colonne de réels
    // Sortie : xc, yc, réels, coordonnées du centre dans le plan
    // r : rayon du cercle
    
    // première estimation par la méthode linéaire de Kasa et Coope
    [xc0, yc0, r0] = rc_kasa_coope(X, Y);
    //disp("Régression circulaire, Kasa et Coope : ")
    //disp("C ("+string(xc0)+" ; "+string(yc0)+") ; r = "+string(r0));
    
    // régression par les moindres carrés totaux
    Ainit = [xc0, yc0, r0];
    [sigma, Aopt] = leastsq(list(residus, X, Y), Ainit);
    xc = Aopt(1) ; yc = Aopt(2) ;
    r = Aopt(3);
    disp("Régression circulaire, moindrees carrés totaux : ")
    disp("C ("+string(xc0)+" ; "+string(yc0)+") ; r = "+string(r0));
endfunction

chemin_fcts = "mypath/";

chdir(chemin_fcts);

n = 24; // nombre de points
pas = 2*%pi/n;

xc = 0; yc = 0; zc = 0; // centre du cercle
R = 1; // rayon du cercle

dx = 5e-2; dy = 1e-4; dz = 2e-2; // décalage du centre

N = 3; // ordre de symétrie de la déformation
dR = 1e-3; // demie amplitude de la déformation radiale
dY = 5e-4; // " axiale 

dThetaX = 2e-3*%pi/180; // inclinaison selon l'axe x
dThetaZ = 1e-3*%pi/180; // inclinaison selon l'axe z
etR = 8e-4; // écart type sur R
etY = 1e-4; // écart type sur Z

phi = linspace(0, 2*%pi-pas, n)';

rayon = R*(1 + dR*cos(N*phi) + etR*rand(n, 1)); //
X = rayon.*cos(phi)*cos(dThetaZ) + dx;
Y = (etY*rand(n, 1) + dY*cos(N*phi))*cos(dThetaX)*cos(dThetaZ) + dy
Z = rayon.*sin(phi)*cos(dThetaX) + dz;

matrice = [X, Y, Z];

csvWrite(matrice, "cercle_deforme.csv")

// **********
// Constantes
// **********

echelle = 100; // amplification de l'écart au cercle
chemin_fcts = "mypath/";

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

// cercle idéal

xc0 = 0; yc0 = 0; zc0 = 0;
r0 = 1;

// cercle déformé

chdir(chemin_fcts);
matrice = csvRead("cercle_deforme.csv")
X1 = matrice(:, 1);
Y1 = matrice(:, 2);
Z1 = matrice(:, 3);

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

exec("fonctions_cercle_mmt.sci"); // chargement des fonctions

function [Xp, Yp, xc, yc, r, B, N, ux, uy] = cercle_projete(X, Y, Z)
    // But : détermine le plan moyen,
    // projette les points expérimentaux
    // et détermine les caractéristique du meilleur cercle
    // Entrées : X, Y, Z, vecteurs colonne de réels
    // Sorties : Xp, Yp, vecteurs colonne de réels
    // coordonnées des points projetés
    // xc, yc, r, réels
    // paramètres du cercle
    
    // Projection sur le plan moyen
    [Xp, Yp, B, N, ux, uy] = regression_plan(X, Y, Z);
    // Régression circulaire
    [xc, yc, r] = regression_circulaire(Xp, Yp);
endfunction

function [] = trace_cercle_projete(xc, yc, zc, r, B, N, ux, uy, ech)
    // But : tracer le cercle de référence
    // et sa projection sur le plan moyen
    // Entréées : xc, yc, zc, r, paramètres du cercle, réels
    // B : intersection du plan moyen avec l'axe y, réel
    // N : vecteur normal au plan moyen, vecteur ligne de réels
    // ux, uy : BON du plan moyen, vecteurs ligne de réels
    // ceh : facteur d'échelle, réel
    
    // cercle de base
    theta1 = linspace(0, 0.45*%pi, 9);
    theta = [theta1, theta1+%pi/2, theta1+%pi, theta1+3*%pi/2, 2*%pi];
    X0 = xc + r*cos(theta); Y0 = yc*ones(1, 37); Z0 = zc + r*sin(theta);
    param3d(X0, Y0, Z0);
    h1 = gce();
    h1.line_style = 3;
    param3d(xc, yc, zc);
    h2 = gce();
    h2.mark_style = 1;
    
    // cercle projeté
    BA = [X0 - B(1) ; Y0 - B(2) ; Z0 - B(3)] // vecteur BA, A pt exp
    B2prime = yc + ech*(B(2) - yc);
    Ux = ux*BA ; Uy = uy*BA; // projections sur le plan, vect ligne
    X1 = B(1) + ux(1)*Ux + uy(1)*Uy;
    Y1 = B2prime + ech*(ux(2)*Ux + uy(2)*Uy);
    Z1 = B(3) + ux(3)*Ux + uy(3)*Uy;
    param3d(X1, Y1, Z1)
    h3 = gce();
    h3.foreground = 2;
    h3.line_style = 4;
    param3d(xc, B2prime, zc)
    h4 = gce();
    h4.mark_foreground = 2;
    h4.mark_style = 1;
endfunction

function [] = trace_cercle_ax(X, Y, Z, xc, yc, zc, r, ech)
    // But : trace les points expérimentaux en amplifiant l'écart radial
    // d'un facteur ech
    // Entrées : X, Y, Z, vecteurs colonne de réels, points relevés
    // xc, yc, zc, r, réels, paramètres du cercle idéal
    // ech, réel, facteur d'échelle
    
    // points relevés
    ymoy = mean(Y); // altitude moyenne
    Yamplifie = yc + ech*(Y - yc);
    param3d([X ; X(1)], [Yamplifie ; Yamplifie(1)], [Z ; Z(1)], 90, 5);
    h2 = gce();
    h2.foreground = 2;
    h2.mark_style=1;
    h2.mark_foreground=2;
    // titre
    xtitle("Écart axial ×"+string(ech));
endfunction

function [] = trace_cercle_rad(X, Y, xc, yc, r, ech)
    // trace le cercle idéal projeté
    // et les points expérimentaux projetés en amplifiant l'écart radial
    // d'un facteur ech
    // Entrées : X, Y, vecteurs colonne de réels, points relevés projetés
    // xc, yc, r, réels, paramètres du cercle idéal projeté
    // ech, réel, facteur d'échelle
    
    // cercle idéal
    diam = 2*r;
    xarc(xc - r, yc + r, diam, diam, 0, 360*64)
    h1 = gce();
    h1.line_style = 3;
    
    // cercle réel
    R = ((X-xc).^2 + (Y-yc).^2).^0.5; // distance de chaque point au centre
    angle = atan(Y-yc, X-xc);
    R1 = ech*(R-r) + r; // distance avec écart amplifié
    XX = R1.*cos(angle) + xc;
    YY = R1.*sin(angle) + yc;
    plot(XX, YY, "+") // tracé des points
    plot([XX ; XX(1)], [YY ; YY(1)], "-") // relie les points
    a = get("current_axes");
    a.isoview = "on";
    xtitle("Écart radial ×"+string(ech), "x", "y");
endfunction

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

// écart axial, précision micromètre
ecart_axial = Y1-yc0;
eamax = round(1000*max(ecart_axial))/1000;
eamin = round(1000*min(ecart_axial))/1000;

// Déformée

// Projection sur le plan moyen
// et détermination du centre du cercle idéal
[Xproj1, Yproj1, xc1, yc1, r1, B, N, ux, uy] = cercle_projete(X1, Y1, Z1)

// Angle par rapport au vecteur y, arrondi au 1/1000
alpha = 0.001*round(acos(N*[0;-1;0])*180000/%pi);
disp("Défaut d''orientation : "+string(alpha));

// tracé en perspective des points par rapport à leur position d'origine
// écart axial (altitude) amplifié
f0 = scf(0); clf;
//subplot(1, 2, 1);
trace_cercle_projete(xc0, yc0, zc0, r0, B, N, ux, uy, echelle)
trace_cercle_ax(X1, Y1, Z1, xc0, yc0, zc0, r0, echelle);

// Tracé du cercle idéal et des points déplacés
f1 = scf(1); clf;
//subplot(1, 2, 2);
trace_cercle_rad(Xproj1, Yproj1, xc1, yc1, r0, echelle);

// écart radial, précision micromètre
disp("Déformée");
ecart_radial = ((Xproj1 - xc1).^2 + (Yproj1 - yc1).^2).^0.5 - r0;
ermax = round(1000*max(ecart_radial))/1000;
ermin = round(1000*min(ecart_radial))/1000;

// tracé de l'écart en fonction de l'angle

// Affichage résultats
disp("Écart axial : "+string(eamin)+" ≤ ea ≤ "+string(eamax));
disp("Écart radial : "+string(ermin)+" ≤ er ≤ "+string(ermax));