% analysis of RLC circuit in series configuration

clc
clear all
close all

% define the parameters
R1=100;
C1=1e-6;
L=1e-3;

% define state space representation for RLC
A1=[0 1/C1; -1/L -R1/L ];
B=[0; 1/L];
C=[1 0];
D=0;

sys1=ss(A1,B,C,D);

% define the transfer function from the ss repr.
[num1, den1]=ss2tf(A1,B,C,D);

% plot the step response for sys1 (R=R1=100 Ohm)
figure(1)
step(sys1)
hold on
grid on

% define the value of Rc for getting two multiple poles (i.e., equal, p1=p2)  
% i.e., Delta of the characteristic equation =0, zita=1
Rc=2*sqrt(L/C1);

R2=Rc/10; % for this value, the poles are complex, and zita<1 (=zita2, see below)
% define ss for RLC with R2 value
A2=[0 1/C1; -1/L -R2/L ];
sys2=ss(A2,B,C,D);

% define the damping factor, zita2, for sys2,  (zita2=0.1)
 [wn2,zeta2,vp2] = damp(sys2);
zita2=zeta(1);
omega_n=wn2(1); % the natural frequency
T_osc=2*pi/omega_n; %the oscillation period  when zita=0


% plot the step response for sys2 (R2=Rc/10)
figure(1)
hold on 
step(sys2)

% check the poles for both systems (sys1 and sys2)...
vp_1=pole(sys1);  % vp_1=eig(A1);
vp_2=pole(sys2);  % vp_2=eig(A2);

% and the static gain
g1=dcgain(sys1);
g2=dcgain(sys2);


R3=0; % poles are purely immaginary, i.e., zita=0
% define ss for RLC with R3=0;
A3=[0 1/C1; -1/L -R3/L ];
sys3=ss(A3,B,C,D);

% plot the step response for sys3 (R3=0)
figure(1)
hold on
step(sys3)
grid on
xlim([0 10*T_osc])
legend

% Bode diagrams for sys2 (zita2=0.1)
figure
bode(sys2)
hold on
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% evaluate the frequency behavior - RLC circuit as a filter %%%
%%% assuming e.g. sys2, the case with R2=Rc/10, zita=0.1 %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% response to a sin. signal (u1) with omega << omega_n
omega1=omega_n/10;
T1=2*pi/omega1;
time1=0:T1/50:5*T1;
u1=sin(omega1*time1);
y_1=lsim(sys2,u1,time1);
figure
lsim(sys2,u1,time1)
hold on
grid on
str_title1=['response to a sin. signal with {\omega} =', num2str(omega1),' rad/s'];
title(str_title1)


% response to a sin. signal (u2) with omega >> omega_n
omega2=10*omega_n;
T2=2*pi/omega2;
time2=0:T2/50:80*T2;
u2=sin(omega2*time2);
figure
lsim(sys2,u2,time2)
hold on
grid on
str_title2=['response to a sin. signal with {\omega} =', num2str(omega2),' rad/s'];
title(str_title2)

% response to a sin. signal (u3) with omega=omega_n
% in this case, the response has a peak amplitude approx. to 1/(2*zita*sqrt(1-2*zita^2))
time3=0:T_osc/50:10*T_osc;
u3=sin(omega_n*time3);
figure
lsim(sys2,u3,time3)
hold on
grid on
str_title3=['response to a sin. signal with {\omega} =', num2str(omega_n),' rad/s'];
title(str_title3)