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

clear;
chdir('monchenim/')

// 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];

// Paramètres de Newton-Raphson
precision = 1e-7; // condition d'arrêt
itermax = 60; // idem

Ainit = [0,0, 1] // cercle initial

// **********
// fonctions
// **********

function [ACx, ACy, AC] = caracteristiques(Xexp, Yexp, A)
    // calcule les caractéristiques du modèle
    // par rapport aux points expérimentaux
    // C(1) : xcentre (scalaire)
    // C(2) : ycentre (scalaire)
    // Xexp, Yexp : points mesurés (vecteur de nombres)
    // ACx : vecteur de composante x des vecteurs AC>
    // " y "
    // AC : longueurs des vecteurs AC>
    ACx = (A(1) - Xexp)';
    ACy = (A(2) - Yexp)';
    AC = ((ACx.^2 + ACy.^2).^(0.5));
endfunction

function [S] = residus(Xexp, Yexp, C)
    // Xexp, Yexp : coordonnées des points expérimentaux (vecteurs)
    // C : coordonnées du centre (vecteur)
    [ACx, ACy, AC] = caracteristiques(Xexp, Yexp, C);
    res = abs(C(3) - AC);
    S = sum(res);
endfunction

function []=affichage(X, Y, A, R, k)
    clf;
    subplot(2,1,1)
    plot2d(X, Y, style = -1) //, frameflag=4)
    isoview(-1.5, 7.5, -3.5, 5.5);
    xstring(-1.8, 5.3, 'C('+string(A(1))+' ; '+string(A(2))+') ; r = '+string(A(3)))
    // modèle
    plot(A(1), A(2), 'r+')
    diametre = 2*A(3);
    xarc(A(1) - A(3), A(2) + A(3),...
    diametre, diametre,...
    0, 360*64)
    a = get('hdl'); // ellipse
    a.foreground = 5; // couleur
    subplot(2,1,2)
    plot2d(R, rect=[1, 0, 8, 33])
    xstring(0.5, 0.1, 'α = '+string(k)+' ; R = '+string(R($)))
//    halt
endfunction

function [A, R]=regression_circulaire(Xexp, Yexp, Ainit, imax, epsilon)
    // Xexp, Yexp : points expérimentaux
    // A(1) : xcentre
    // A(2) : ycentre
    // A(3) : rayon
    // Ainit : valeur initiale de A
    // imax, epsilon : condition d'arrêt
    i = 1;
    k = 1;
    A = Ainit';
    dA = [1 1 1]';
    drapeau = %t; // pour la 1re itération
    R(1) = residus(Xexp, Yexp, A')
    affichage(Xexp, Yexp, Ainit, R, 1)
    nom = 'cercle_gruntz'+string(i)+'.png';
    xs2png(0,nom)
    while drapeau & (i <= imax)
        i = i+1;
        [ACx, ACy, AC] = caracteristiques(Xexp, Yexp, A);
        J = [ACx./AC, ACy./AC, -ones(ACx)]; // matrice des vecteurs unité
        f = -(AC - A(3)*(ones(Xexp))'); 
        dA = J\f;
        drapeau2 = %t // pour une 1re exécution
        while drapeau2 & (k>=0.1) // teste la divergence sur 1 étape
            A = A + k*dA; 
            R(i) = residus(Xexp, Yexp, A')
            drapeau2 = (R(i) >= R(i-1)) // vrai si diverge
            if drapeau2 then k = k*0.75; // atténue si diverge
            else k0 = k; // pour affichage de la valeur
                k = (1 + k)/2; // réduit l'atténuation si converge
            end
        end
        drapeau = (abs(R(i) - R(i-1)) > epsilon)
        disp('i = '+string(i)+' ; k = '+string(k0)+' ; R = '+string(R(i)))
        affichage(Xexp, Yexp, A', R, k0)
        nom = 'cercle_gruntz'+string(i)+'.png';
        xs2png(0,nom)
    end 
endfunction

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

// lecture des données

// Xdef, Ydef : vecteurs ligne

Xdef = X0;
Ydef = Y0;

fenetre = scf(0); // création de la fenêtre graphique
fenetre.figure_size = [400,800];

// regression
[Aopt, res] = regression_circulaire(Xdef, Ydef, Ainit, itermax, precision)
centre = Aopt(1:2);
rayon = Aopt(3);
R = res($);

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