// Simulation de l'évolution de l'énergie totale d'un circuit RLC
// Dominique LEFEBVRE - septembre 2010
// www.tangenteX.com

// Définition de la fonction à intégrer : d²u/dt² = -(2/tau)*du/dt - omega0²*u
// Il s'agit d'une dérivée d'ordre 2, on procède donc à un changement de variables
// f(1) = du/dt
// f(2) = -(R/L)*du/dt - (1/LC)*u
function du = RLC(t,u)
  du(1) =  u(2);                          // du/dt                    
  du(2) = -(2/tau)*u(2) - omega2*u(1);    // d²u/dt²
endfunction;

// Saisie des paramètres de la simulation
E = input('Tension de charge du condensateur (en V): ');
R = input('Résistance (en ohm - R > 0): ');
C = input('Capacité (en farad): ');
L = input('Inductance (en henry): ');

// Définition des conditions initiales de l'intégration
t0 = 0 ;                         // origine des temps = 0 s
u0 = [E;0];                      // la tension initiale est la tension de charge du condensateur

// Autres paramètres de la simulation
NbPeriode = 1;
t1 = NbPeriode*2*%pi*sqrt(L*C);         // date de la fin de l'expérience, fixée à 10 fois la valeur de 2*pi*sqrt(LC) (en seconde)
pas = 1000;                     // nombre de pas d'intégration
dt = (t1 - t0)/pas;             // pas de calcul d'intégration
mprintf('Période: %f\n',sqrt(L*C));

// Initialisation du vecteur de temps nécessaire à l'intégration
t=t0:dt:t1;

// calcul des paramètres du circuit RLC
omega2 = 1/(L*C);                // pulsation propre au carré
tau = 2*L/R                      // facteur d'amortissement
Q = (1/R)*sqrt(L/C)             // facteur de qualité          
mprintf('Pulsation propre: %f - Amortissement: %f - Facteur de qualité du circuit: %f\n',sqrt(omega2),tau, Q);

// affichage du régime de l'oscillateur
if (Q > 0.5) then mprintf('Regime de variation pseudo-periodique\n'); end;
if (Q < 0.5) then mprintf('Regime de variation aperiodique\n'); end;
if (Q == 0.5) then mprintf('Regime de variation critique\n'); end; 
     

// Résolution de l'EDO
// i est le vecteur contenant la valeur de l'intensité calculée pour chaque pas
u = ode(u0,t0,t,RLC);

// Calcul des énergies Ec et Eb (condensateur, bobine)

Ec = (C/2)*u(1,:).^2;
Eb = (L*C^2/2)*u(2,:).^2;
Et = Ec + Eb;

// Tracé de la solution
clf;                             // clear de la fenetre graphique
xtitle('Evolution des énergies dans un circuit RLC','Temps (s)','Energie (J)');
plot2d(t,Ec, style= 1);
plot2d(t,Eb, style= 2);
plot2d(t,Et, style= 5);