//
// Studio delle configurazioni di equilibrio di un punto materiale
// pesante P, mobile su una guida liscia contenuta in un piano verticale
// e soggetto all'azione di una molla. La molla, per la quale si suppone
// valida la Legge di Hooke con lunghezza a riposo zero, ha l'estremo
// non collegato a P fisso in un punto A del piano.
//
// La massa m del punto e la costante elastica K della molla sono
// assegnate, come anche le coordinate del punto A e l'equazione della
// guida rispetto ad un fissato sistema di riferimento nel piano con
// asse x orizzontale verso destra e z verticale verso l'alto.
//
// Incognita del problema: ascissa delle configurazioni di equilibrio
// del punto P.
//
// Essendo il vincolo liscio e le forze attive conservative, le
// configurazioni di equilibrio sono punti stazionari della funzione
// energia potenziale
//
//                     1
//   V(x) = mg f(x) + --- K ( (x - x(A))^2 + (f(x) - z(A))^2 )
//                     2
//
// La funzione f che descrive la guida deve essere derivabile.
//
// **************************************************** Costanti:
m = 1; // kg
g = 9.8; // m/s^2
K = 10; // N/m
coord_A = [0,4]; // m, coordinate di A
// **************************************************** Fine costanti
//
protezione = funcprot();
funcprot(0);
//
// Funzione che descrive la guida e sua derivata
//
function z = f_guida(x)
    //z = (x.^2)/2;  // (1)
    z = 1 - cos(x)   // (2)
endfunction
//
function d1 = d1f_guida(x)
    //d1 = x    // (1)
    d1 = sin(x) // (2)
endfunction
//
// Energia potenziale del punto
//
function ep = V_punto(x)
    ep = m*g*f_guida(x) +...
         K*((x-coord_A(1)).^2 + (f_guida(x)-coord_A(2)).^2)/2
endfunction
//
// Derivata della funzione energia potenziale
//
function y = d1V_punto(x)
    y = m*g*d1f_guida(x) +...
         K*((x-coord_A(1)) + (f_guida(x)-coord_A(2)).*d1f_guida(x))
endfunction
//
// Ricerca intervalli che includono una configurazione di equilibrio.
// La ricerca avviene utilizzando il Teorema di esistenza degli zeri
// sulla funzione V'.
//
PuntiEq = []; // Colonna dei punti di equilibrio trovati (in m)
Intervalli = []; // Matrice le cui colonne sono gli intervalli (in
                 // teoria, ciascun intervallo include un punto di eq)
//
// *** Qui si sceglie in quale intervallo globale ricercare i punti
//     di equilibrio!
xx = linspace(-2*%pi,2*%pi,401); // I potenziali estremi di intervalli:
                         // si cercano tra questi zeri di V' oppure
                         // punti adiacenti nei quali V' assume valori
                         // di segno opposto.
// ***
Segnod1V = sign(d1V_punto(xx)); // sign(x) = 0 se x = 0 ecc...
if Segnod1V(1) == 0 then PuntiEq(1) = xx(1); end;
for i = 2:length(xx),
    if Segnod1V(i) == 0 then PuntiEq($+1) = xx(i);
       elseif Segnod1V(i-1) ~= 0 then
           if Segnod1V(i) ~= Segnod1V(i-1) then
              Intervalli(1:2,$+1) = [xx(i-1);xx(i)]; end;
    end;
end;
// Adesso ho qualche configurazione di eq (quelle trovate tra le xx)
// e intervalli che ne racchiudono altre
//
if PuntiEq ~= [] then
    for PE = PuntiEq,
            printf("\nx_eq = %3.2e\n",PE);
    end
end
//
// **************************** Ricerca configurazioni di eq
//
exec('..\Metodi\bisezione.sci');
tol = 2d-6; // Individua ciascun punto di eq con errore assoluto
            // inferiore ad 1d-6 m (un micron)
kmax = 55;
//
for Intervallo = Intervalli,
    a = Intervallo(1); b = Intervallo(2);
    [eq,v,info] = bisezione(d1V_punto,a,b,"assoluto",tol,kmax,"zitto");
    PuntiEq($+1) = eq;
    printf("\nx_eq = %3.2e",eq);
    printf("\nd1V(x_eq) = %3.2e",v);
    printf("\ninfo = %d\n",info);
end;
// **************************** Fine ricerca configurazioni di eq
//
printf("\nConfigurazioni di equilibrio trovate: %d", length(PuntiEq));
//
// Grafico energia potenziale e derivata, con configurazioni eq trovate
//
scf(0); clf(0);
subplot(221);
plot(xx,V_punto(xx),"b");
for eq = PuntiEq',
    plot(eq,V_punto(eq),"ro");
end;
xlabel("x");
ylabel("energia potenziale");
xgrid();
subplot(223);
plot(xx,d1V_punto(xx),"b");
for eq = PuntiEq',
    plot(eq,d1V_punto(eq),"ro");
end;
xlabel("x");
ylabel("derivata dell''energia potenziale");
xgrid();
//
// Grafico configurazioni di equilibrio
//
subplot(222);
plot2d(xx,f_guida(xx),frameflag=4);
for eq = PuntiEq',
    plot2d([eq,coord_A(1)],[f_guida(eq),coord_A(2)],3);
    DisegnoMolla = gce(); // il compound...
    LineaDisegnoMolla = DisegnoMolla.children; // ...la linea.
    LineaDisegnoMolla.thickness = 3;
    plot(eq,f_guida(eq),"ro");
end;
plot(coord_A(1),coord_A(2),"ko");
xlabel("x");
ylabel("z");
xtitle("configurazioni di equilibrio");
xgrid();
//
funcprot(protezione);
//
