% Circuit RL. ELE4801. Jesus Morales. Hiver 2018.
close all; clear; clc
%------------ Donnes d'entré -------------------------------

R = 1;         L = 10e-3;                   % circuit RL
Vm = 1;        f = 60;            th = 30;  % source

dt = 2e-6;       % pas de integration
Tmax = 4/60;     % temps de simulation

%------------ Modèle de space d'etat -----------------------

A = -R/L;      B = 1/L;          

%------------ Discretisation (règle trapezoïdal) -----------

P = (1-dt/2*A)^-1 * (1+dt/2*A);
Q = (1-dt/2*A)^-1 * (dt/2*B);

%------------ Initialisation de variables ------------------
 
w = 2*pi*f;                 % fréquence angulaire de la source
t = 0:dt:Tmax;              % vecteur de temps
u = Vm*cos(w*t+th*pi/180);  % vecteur d'entrée 
N = length(t);              % nombre de samples
x = zeros(N,1);             % vecteur solution

%------------ Solution en régime permanent ------------------

% Solution en régime permanent avec phaseurs
% s = 1j*w;
% V_jw = Vm * ( cosd(th) + 1j*sind(th) );
% I_jw = V_jw/(R+s*L);

% Solution en régime permanent dans le domaine du temps
% it = abs(I_jw)*cos(w*t+angle(I_jw));  

% initialisation de la simulation en régime permanent 
% x(1) = real(I_jw);   

%------------ Integration numérique -------------------------

for k = 1:N-1
    x(k+1) = P * x(k) + Q * ( u(k) + u(k+1) ) ;
end

%------------ Graphique des résultats -----------------------

% figure (1); plot(t,u); grid; title('Entrée')
% xlabel('Temps (s)');  ylabel('Tension (V)')

figure (2);   
plot(t,x); grid; hold on;
% plot(t,it,'--'); 
% legend('Solution en régime transitoire','Solution en régime permanent')
title('Sortie');  xlabel('Temps (s)');  ylabel('Courant (A)')


%------------ Solution analytique ---------------------------

% beta = atan(w*L/R);
% 
% i_reg_permanent = Vm / ( sqrt(R^2+w^2*L^2) ) * ( cos(w*t + th*pi/180 - beta) );
% i_reg_transitoire = -Vm / ( sqrt(R^2+w^2*L^2) ) * ( cos(th*pi/180 - beta) ) * exp(-R*t/L);
% 
% i_analytique = i_reg_permanent + i_reg_transitoire;
% 
% figure; plot(t,i_reg_permanent); xlabel('Temps (s)');  ylabel('Courant (A)'); 
% title('Solution en régime permanent'); grid  
% 
% figure; plot(t,i_reg_transitoire); xlabel('Temps (s)');  ylabel('Courant (A)'); 
% title('Solution transitoire'); grid 
% 
% figure; plot(t,i_analytique); xlabel('Temps (s)');  ylabel('Courant (A)'); 
% title('solution totale'); grid  
%