//
// Esempio di cinematica (da Quarteroni, Saleri: Calcolo Scientifico,
// Pb. 2.3 pagina 42).
//
// Lunghezza delle aste (cm)
//
a1 = 10;
a2 = 13;
a3 = 8;
a4 = 10;
//
function y = f_cine(x)
    y = a1^2 + a2^2 - a3^2 + a4^2 + 2*a1*a4*cos(Beta)...
          - 2*a1*a2*cos(x) - 2*a2*a4*cos(x-Beta)
endfunction
//
function j = Dxf_cine(x)
    j = 2*a1*a2*sin(x) + 2*a2*a4*sin(x-Beta)
endfunction
//
Percorso = "PERCORSO newton1d.sci";  // <------ *** MODIFICARE ***
exec(Percorso + "newton1d.sci");
//
E_newt = 1d-5;
kmax = 20;
//
N = 100; // numero di incrementi di Beta
Beta_max = 2.636; // rad
dBeta = Beta_max/N; // incremento
x0_a = -0.24; // valore iniziale per soluzione a (per B = 0)
x0_b = 0.24; // valore iniziale per soluzione b (per B = 0)
xx_a = zeros(1,N+1); // riga delle soluzioni a
iinfo_a = zeros(1,N+1); // riga degli esiti per soluzione a
xx_b = zeros(1,N+1); // riga delle soluzioni b
iinfo_b = zeros(1,N+1); // riga degli esiti per soluzione b
for k = 1:N+1,
    Beta = (k-1)*dBeta;  // Beta = 0, dBeta,..., N*dBeta = Beta_max
    [xx_a(k),v,iinfo_a(k),XX_a] =... // cerca soluzione a
       newton1d(f_cine,x0_a,Dxf_cine,E_newt,kmax,"zitto");
    x0_a = xx_a(k); // punto iniziale prossima iterazione per sol. a
    [xx_b(k),v,iinfo_b(k),XX_b] =... // cerca soluzione b
       newton1d(f_cine,x0_b,Dxf_cine,E_newt,kmax,"zitto");
    x0_b = xx_b(k); // punto iniziale prossima iterazione per sol. b
end;
//
// Controllo che ciascuna applicazione di Newton sia andata a buon fine
//
if find(iinfo_a > 1) == [] then printf("\nsoluzione a: OK!");
   else printf("\nproblemi nella soluzione a!") end;
if find(iinfo_b > 1) == [] then printf("\nsoluzione b: OK!\n");
   else printf("\nproblemi nella soluzione b!\n") end;
//
// Grafico
//
scf(0); clf(0); show_window(0);
plot(dBeta*[0:N],xx_a,"-r",dBeta*[0:N],xx_b,"-b");
xlabel("beta");
ylabel("alfa");
legend("soluzione a","soluzione b")
xgrid();
//
// Aspetta segni di vita dall'utilizzatore
//
halt("\n\n *** Premere invio per avviare l''animazione");
//
// Inizio animazione
//
scf(77); show_window(77);
han = gcf();
han.pixmap = "on"; // Avvio del modo per l'animazione:
                   // prima genera gli oggetti da disegnare
                   // poi, con "show_pixmap", li mostra
for i = 1:N+1,
    attesa = 80;
    sleep(attesa); // si gingilla per attesa ms allo scopo
                   // di rallentare l'animazione
    clf();
    Beta = (i-1)*dBeta;
    subplot(121);
    xtitle("Animazione del quadrilatero articolato (a)");
    Aa = [a2*cos(xx_a(i)); a2*sin(xx_a(i))];
    Ba = [a1 + a4*cos(Beta);  a4*sin(Beta)];
    plot2d(%nan,%nan,frameflag=3,rect=[-0.1,-5,21,15]);
    xpoly([0,Aa(1),Ba(1),a1],[0,Aa(2),Ba(2),0],"lines");
    quadrilatero_a = gce();
    quadrilatero_a.thickness = 3;
    quadrilatero_a.foreground=5;
    giunti_mobili_ax = [Aa(1),Ba(1)];
    giunti_mobili_ay = [Aa(2),Ba(2)];
    plot2d(0,0,-3); plot2d(a1,0,-3); // giunti fissi
    plot2d(giunti_mobili_ax,giunti_mobili_ay,-9);
    subplot(122);
    xtitle("Animazione del quadrilatero articolato (b)");
    Ab = [a2*cos(xx_b(i)); a2*sin(xx_b(i))];
    Bb = [a1 + a4*cos(Beta);  a4*sin(Beta)];
    plot2d(%nan,%nan,frameflag=3,rect=[-0.1,-5,21,15]);
    xpoly([0,Ab(1),Bb(1),a1],[0,Ab(2),Bb(2),0],"lines");
    quadrilatero_b = gce();
    quadrilatero_b.thickness = 3;
    quadrilatero_b.foreground=2;
    giunti_mobili_bx = [Ab(1),Bb(1)];
    giunti_mobili_by = [Ab(2),Bb(2)];
    plot2d(0,0,-3); plot2d(a1,0,-3); // giunti fissi
    plot2d(giunti_mobili_bx,giunti_mobili_by,-9);
    show_pixmap();  // mostra il disegno
    if i == 1 then sleep(500); end;
end;
f.pixmap = "off";
//
