%% Respuesta de circuito RC ante entrada CC y senoidal, por Laplace
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentaciones en vídeo: http://personales.upv.es/asala/YT/V/rcstepsin.html
% Este código funcionó sin errores en Matlab R2022b (Linux)
% Objetivo: obtener respuestas ante entradas escalón y senoidal de tensión de condensador e intensidad en un circuito RC serie, mediante el método de Laplace. Supondremos condiciones iniciales nulas.
%
% Tabla de Contenidos
% Modelado
% Análisis comportamiento Control Systems Toolbox
%     Comportamiento en corriente continua
%     Comportamiento en corriente alterna (frecuencia variable)
% Comportamiento por Laplace (symbolic toolbox)
%     En CC (entrada escalón)
%     En CA
%     Comparación con otros enfoques
%
%
%% Modelado
%
% El modelo es
% \frac{dV_c}{dt}=\frac{1}{C}\cdot \underbrace{\frac{(u-V_c)}{R}}_{i} (ecuación de estado)
% y vamos a calcular tanto tensión de condensador como intensidad como salidas de interés:
% y=\left(\begin{array}{c}V_c \\ \frac{1}{R}(u-V_c)\end{array}\right) (ecuación de salida)

R=20; C=2e-4;

% Representación matricial \dot x=Ax+Bu, y=C_ex+Du

A=-1/(R*C); B=1/(R*C);
Ce=[1;-1/R]; D=[0; 1/R]; %C ya está gastado, llamo Ce

%% Análisis comportamiento Control Systems Toolbox

circuito=ss(A,B,Ce,D);
circuito.InputName='u';
circuito.OutputName={'V_c','i'};
tf(circuito)

%% Comportamiento en corriente continua

step(circuito,0.05), grid on

%% Comportamiento en corriente alterna (frecuencia variable)

bode(circuito), grid on

%% Comportamiento por Laplace (symbolic toolbox)

syms s %variable de Laplace
MdT=simplifyFraction(Ce*inv(s*eye(1)-A)*B+D)

% El término de condiciones iniciales Ce\cdot (sI-A)^{-1}x(0) es nulo, según explicitado al enunciar el problema.
%% En CC (entrada escalón)

Y_s=MdT*1/s
Y_partfrac=partfrac(Y_s)
Y_t=ilaplace(Y_partfrac)
fplot(Y_t,[0 0.05],LineWidth=2), grid on, legend("V_c","i",Location="best")

%% En CA

syms t real
%syms omega real %genérico
%omega=100*sym(pi); %los 50Hz usuales
omega=200; %ejemplo con números más "redondos"
U_t=sin(omega*t);
U_s=laplace(U_t)
Y_s=MdT*U_s
Y_partfrac=partfrac(Y_s,s) %fracciones en la variable de Laplace, no omega si fuera genérica
vpa(Y_partfrac,5)
Y_t=vpa(ilaplace(Y_partfrac),4)
fplot([Y_t;U_t],[0 0.1]), grid on, legend("V_c","i","u",Location="best")

% En régimen estacionario senoidal, la tensión del condensador va retrasada en fase:

atan(-0.4878/0.6098)*180/pi %(grados)

% Y la intensidad va "adelantada":

atan(0.0243/0.01951)*180/pi %(grados)

%% Comparación con otros enfoques
% La amplitud de la intensidad coincide con las fórmulas de impedancia en circuitos eléctricos:

Z= R + -1/(C*omega)*1j
1/Z %intens=tension/impedancia
angle(ans)*180/pi %el adelanto o retraso en fase (grados)

% Y con las fórmulas genéricas de respuesta en frecuencia G(j\omega) para sistemas lineales (la impedancia es un caso particular de FdT cuando entrada es corriente y salida es tensión eléctrica):

vpa(subs(MdT(2), s, 1j*omega), 3)