// Initialisation

clear;

// données (D. Gruntz)

X0 = [0.7, 3.3, 5.6, 7.5, 6.4, 4.4, 0.3, -1.1];
Y0 = [4.0, 4.7, 4.0, 1.3, -1.1, -3.0, -2.5, 1.3];

// fonctions

// régression circulaire linéaire

function [C, r, R]=regression_circulaire(X, Y)
    // X, Y : points expérimentaux (vecteurs ligne)
    // C : coordonnées du centre du cercle (vecteur)
    // r : rayon du cercle (scalaire)
    A = [X ; Y];
    d = X.^2 + Y.^2;
    // xetoile = 2*xcentre
    // yetoile = 2*ycentre
    // zetoile = r^2 - (xcentre^2 + ycentre^2)
    [xyetoile, zetoile, sigma] = reglin(A, d');
    // A'\d' si l'on ne s'intéresse pas à l'écart
    C = 0.5*xyetoile;
    r = sqrt(zetoile + C(1)^2 + C(2)^2);
    R = sigma;
endfunction

// **********
// programme principal
// **********

// lecture des données

// Xdef, Ydef : vecteurs ligne
Xdef = X0;
Ydef = Y0;

// regression
[centre, rayon, ecart] = regression_circulaire(Xdef, Ydef)

// tracé

clf;
// points
plot2d(Xdef, Ydef, style = -1) //, frameflag=4)
isoview(-2, 8, -4, 6);
xstring(-1.5, 5.5, 'C('+string(centre(1))+' ; '+string(centre(2))+') ; r = '+string(rayon))
// modèle
plot(centre(1), centre(2), 'r+')
diametre = 2*rayon;
xarc(centre(1) - rayon, centre(2) + rayon,...
    diametre, diametre,...
    0, 360*64)
a = get('hdl'); // ellipse
a.foreground = 5; // couleur


print(%io(2), centre)
print(%io(2), rayon)
print(%io(2), ecart)