//********************************************************************************
// Modélisation du pendule simple sans frottement
// Dominique Lefebvre - TangenteX.com
// septembre 2010
//********************************************************************************

// Définition des paramètres
l = 1.         // longueur en m
g = 9.81       // accélération de la pesanteur en m.s-2
m = 1.         // masse en kilogramme
Omega2 = g/l;  // pulsation du pendule

// Définition du système différentiel décrivant le modèle du pendule
function du = Pendule(t,u)
  du(1) = u(2);
  du(2) = - Omega2*sin(u(1));
endfunction

// Définition des conditions initiales
theta0 = input("Angle initial (en rd): ");
vtheta0 = input("Vitesse angulaire initiale (en rd/s): ");
E0 = (1/2)*m*(l^2)*vtheta0^2 + m*g*l*(1 - cos(theta0));  // Energie mécanique totale du pendule = Energei initiale

u0 = [theta0;vtheta0];
t0 = 0;

// Définition des paramètres de calcul
T = 10*%pi;
pas = 0.01;
t = t0:pas:T;

// intégration du système
[u] = ode(u0,t0,t,Pendule);  // u(1): angle(t)  u(2): vitesse angulaire(t)

// Affichage de la trajectoire theta(t)
clf;
subplot(2,1,1);
a=gca();
a.x_label.text = "Temps";
a.y_label.text = "Angle";
plot2d(t,u(1,:)',u(2,:)', style = 2, leg = " Trajectoire du pendule simple - Angle initial: " + string(theta0) + " rd - Em = " + string(E0) + " J");
xgrid(2);

// affichage du plan de phase
subplot(2,1,2);
b=gca();
b.x_label.text = "Angle";
b.y_label.text = "Vitesse angulaire";
plot2d(u(1,:)',u(2,:)', style = 2, leg = "Plan de phase");
xgrid(2);