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

clear;
chdir('monchemin/');

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

// paramètres de l'arc de cercle
centre=[1;1];
rayon = 4;
angle = %pi/3:0.1:2*%pi/3;

// paramètres du bruit
var = 0.05;

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

// génération des données
ligne_rayon_bruite = rayon + var*rand(angle,'norm');
foo = [cos(angle);sin(angle)].* ...
    [ligne_rayon_bruite;ligne_rayon_bruite];
arc_cercle_bruite(1,:) = centre(1) + foo(1,:);
arc_cercle_bruite(2,:) = centre(2) + foo(2,:);

clf
plot(arc_cercle_bruite(1,:), arc_cercle_bruite(2,:), '+')
write('arc_cercle_bruite.txt', arc_cercle_bruite')

// **********
// Initialisation
// **********
 
clear;
chdir('monchemin/')
 
// **********
// 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);
    C = 0.5*xyetoile;
    r = sqrt(zetoile + C(1)^2 + C(2)^2);
    // A'\d' si l'on ne s'intéresse pas à l'écart
    R = sigma;
endfunction
 
// **********
// programme principal
// **********
 
// lecture des données
 
// Xdef, Ydef : vecteurs ligne
donnees = read('arc_cercle_bruite.txt', -1, 2)';
Xdef = donnees(1,:);
Ydef = donnees(2,:);
 
// regression
[centre, rayon, ecart] = regression_circulaire(Xdef, Ydef)
 
// tracé
 
clf;
// points
plot2d(Xdef, Ydef, style = -1) //, frameflag=4)
isoview(-2, 4, 1, 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)