%% Respuesta circuito RCR ante  tren de pulsos senoidales (señal seno rectificada)
% © 2023, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/trenpulRCR.html
% *Este código funcionó sin errores en Matlab R2022b (Linux)
% Objetivo: analizar y calcular la transformada de Laplace de un pulso (un semiperíodo) de forma senoidal, para su posible uso en problemas de respuesta temporal o/y rectificadores.
% Calcular la respuesta de un sistema con resistencias y condensadores ante un tren de pulsos senoidales.
% Tabla de Contenidos
% Modelado (repr. interna y FdT)
% Análisis de la entrada tipo pulso senoidal y transf. Laplace
%     Entrada en dominio del tiempo
%     Entrada "único pulso" en dominio de Laplace
%         Tren de pulsos
% Respuesta temporal ante PULSO senoidal único, CI Nulas
%     Sin usar , sólo con la salida ante senoide
%     Respuesta ante CI no nulas
% Respuesta "completa" ante tren de pulsos
% Régimen estacionario ante tren de pulsos senoidales
%
%
%% Modelado (repr. interna y FdT)
%
% En circuito RCR de la figura tiene de ecuaciones:
% \frac{dV_c}{dt}=\frac{1}{C}i_c,     i_c=i-i_L,     i=\frac{1}{R}(u-V_c),     i_L=\frac{1}{R_L}V_c
% Son 4 ecuaciones y 4 "incógnitas": V_c, i, i_c, i_L por lo que el modelo está completo.
% La representación interna con salida V_c será:
% \frac{dV_c}{dt}=\frac{1}{C}\left(\cdot\frac{1}{R}(u-V_c)-\frac{1}{R_L}V_c\right)
% o sea, ecuaciones de estado y salida:
% \frac{dV_c}{dt}=-\left(\frac{1}{RC}+\frac{1}{R_LC}\right)\cdot V_c+\frac{1}{RC}u
% y=V_c
% *El problema se "inspira" en un filtro capacitivo en un rectificador, pero no es el modelo del mismo, porque con diodos la intensidad i no puede ser negativa, o sea que i=\max\left(0,\frac{1}{R}(u-V_c)\right), no lineal.
% Demos valores a parámetros físicos constantes:

R=0.5;RL=20;C=2e-2;
A=-1/(R*C)-1/(RL*C); B=1/(R*C); C=1; D=0;

% En el dominio de Laplace, el modelo será, usando
% Y(s)=(C(sI-A)^{-1}B+D)\cdot U(s)+C(sI-A)^{-1}\cdot x(0)

syms s
FdT=simplifyFraction(C*1/(s-A)*B+D)
TCI=simplifyFraction(C*1/(s-A)) %lo que multiplica a V_c(0)
t_est_98=log(0.02)/A  %duración transitorio (tiempo establecimiento) criterio "98% disipado"

%% Análisis de la entrada tipo pulso senoidal y transf. Laplace
%% Entrada en dominio del tiempo
% La entrada "pulso senoidal" (un semiperíodo \pi/\omega) es
% u(t)=\begin{cases} 0 &t<0 \\ \sin(\omega t) &0\leq t\leq \frac{\pi}{\omega} \\ 0 &t>\frac{\pi}{\omega}\end{cases}
% Por ejemplo, trabajaremos con un segmento (semiperíodo) de onda de 50 Hz de corriente alterna.

omega=100*pi; %la onda de 50 Hz.
SemiPeriod=pi/omega
syms t real
u(t)= heaviside(t)*heaviside(SemiPeriod-t)*sin(omega*t);
fplot(u,[-0.01 0.03],LineWidth=2), grid on, title("Pulso senoidal (semiperíodo) de entrada"),xline(0,LineWidth=2)

%% Entrada "único pulso" en dominio de Laplace
% Matlab sabe su transformada de Laplace:

Pulso(s)=laplace(u)

%% Tren de pulsos
% Si el pulso se repite con un período T, entonces
% Tren_{pulsos}(t)=pulso(t)+pulso(t-T)+pulso(t-2T)+\dots
% Si U(s) es la transformada de Laplace de un único pulso, el tren completo tendrá:
% Tren_{pulsos}(s)=U(s)\cdot (1+e^{-Ts}+e^{-2Ts}+\dots)=U(s)\frac{1}{1-e^{-Ts}}.
% Hemos usado la fórmula de suma de una progresión geométrica 1/(1-r), con razón r=e^{-Ts}.

Tren(s)=Pulso(s)/(1-exp(-1/100*s))
fplot(ilaplace(Tren),[0 0.1]), grid on, title("Tren de pulsos senoidales")

% La respuesta del sistema ante ese tren de pulsos, con condiciones iniciales nulas, será:

coco=simplify(ilaplace(FdT*Tren),400)
fplot(coco,[0 0.1]) %Matlab se "raya" mucho...

%% Respuesta temporal ante PULSO senoidal único, CI Nulas
%% Sin usar U(s), sólo con la salida ante senoide U_1(s)

u1=sin(100*pi*t);
U1s=laplace(u1)
Y1=FdT*U1s
vpa(partfrac(Y1,FactorMode='real'),5)
Y1_t(t)=ilaplace(Y1);
vpa(Y1_t,5)
fplot([u1;Y1_t],[0 0.08]), grid on, legend("U1(t)","Salida (ante senoide completa)")

% Con y1, como la salida ante u2 (que es u1 retrasada) será y1 retrasada, la salida ante el pulso senoidal será:
% y(t)=y_1(t)+\begin{cases}0 &t<0.01 \\ y_1(t-0.01) &t\geq 0.01\end{cases}

Y_t_pulso(t)=Y1_t+subs(heaviside(t)*Y1_t,t,t-0.01);
vpa(Y_t_pulso(t),3)
fplot([u;Y_t_pulso],[0 0.06],LineWidth=2), grid on,legend("U(t)","Salida")

%% Respuesta ante CI no nulas
% Suponemos condensador inicialmente cargado

V0=0.95;

% Podemos usar lo anterior (forzada), porque por superposición, añadimos la respuesta libre y problema terminado:

resp_libre_Laplace=TCI*V0
resp_libre_t=ilaplace(resp_libre_Laplace)
Y_t_conCI=Y_t_pulso+resp_libre_t;
fplot([u;Y_t_conCI],[0 0.06],LineWidth=2), grid on,legend("U(t)","Salida")

% Al final, será esta exponencial por la condición inicial lo que hay que añadir:

cosaCI=ilaplace(TCI)

%
%% Respuesta "completa" ante tren de pulsos
% Se haría a trozos... Viendo en qué condición acaba "cada pulso" y repitiendo cálculos para el siguiente pulso en las CI en las que se quedara.
% Causalidad: como el pulso coincide con el "seno" en los primeros 0.01 segundos, las respuestas en ese intervalo serán idénticas. Trabajaremos con la respuesta Y1_t ante el seno y no con la respuesta ante el pulso completo [en el vídeo que explica este material se gasta Y_t_pulso, pero el resultado es idéntico y más sencillo con Y1_t].
% Por ejemplo, con CI = 0 (inicial de verdad en t=0):

Ytrozo1(t)=Y1_t(t); vpa(Ytrozo1(t),4) %Y1_t es la resp. ante sin(100·pi·t) y condiciones iniciales nulas.
V1=Ytrozo1(0.01); vpa(V1,4)

% Si la CI fuera la "final" del trozo 1:

Ytrozo2(t)=Y1_t(t)+V1*cosaCI; vpa(Ytrozo2(t),4)
V2=Ytrozo2(0.01); vpa(V2,4)

% Si la CI fuera la "final" del trozo 2:

Ytrozo3(t)=Y1_t(t)+V2*cosaCI;
V3=Ytrozo3(0.01); vpa(V3,4)
Ytrozo4(t)=Y1_t(t)+V3*cosaCI;
V4=Ytrozo4(0.01); vpa(V4,4)
Ytrozo5(t)=Y1_t(t)+V4*cosaCI;
V5=Ytrozo5(0.01); vpa(V5,4)
Ytrozo6(t)=Y1_t(t)+V5*cosaCI;
V6=Ytrozo6(0.01); vpa(V6,4)
Ytrozo7(t)=Y1_t(t)+V6*cosaCI;
V7=Ytrozo7(0.01); vpa(V7,4)
Ytrozo8(t)=Y1_t(t)+V7*cosaCI;
V8=Ytrozo8(0.01); vpa(V8,4) %realmente, un bucle "for" parecería una opción más inteligente: copiar-pegar es "herejía" en programación.

% Esto parece que "converge" a un régimen estacionario (no senoidal), que discutiremos abajo.

fplot(Ytrozo1(t),[0 0.01]), grid on, hold on
fplot(Ytrozo2(t-.01),[0.01 0.02]),
fplot(Ytrozo3(t-.02),[0.02 0.03]),
fplot(Ytrozo4(t-.03),[0.03 0.04]), 
fplot(Ytrozo5(t-.04),[0.04 0.05]), 
fplot(Ytrozo6(t-.05),[0.05 0.06]), 
fplot(Ytrozo7(t-.06),[0.06 0.07]), 
plot(0:0.01:0.07,[0 V1 V2 V3 V4 V5 V6 V7],'.k')
hold off

%
%
%% Régimen estacionario ante tren de pulsos senoidales
% Debe acabar en la misma tensión que empezó el ciclo.
% Si empezamos en V_0, entonces

syms V_0 real
Yperm(t)=Y1_t(t)+ilaplace(TCI*V_0); %válido el primer período.
vpa(Yperm(0),4)
vpa(Yperm(0.01),4)
sol=solve(Yperm(0.01)==Yperm(0))

% Pues esa es la condición inicial:

Yfinal(t)=subs(Yperm(t),V_0,sol);
vpa(Yfinal,4)
fplot([u; Yfinal],[0 0.01]), grid on, legend("Pulso de entrada","Régimen permanente NO senoidal")

% Eso se repetiría "ad infinitum"... Tendríamos un régimen estacionario periódico no-senoidal.

fplot(u,[0 0.01],Color="blue"), grid on
hold on
fplot(Yfinal,[0 0.01],LineWidth=2,Color="red"), grid on
fplot(u(t-0.01),[0.01 0.02],Color="blue"), grid on
fplot(u(t-0.02),[0.02 0.03],Color="blue"), grid on
fplot(u(t-0.03),[0.03 0.04],Color="blue"), grid on
fplot(Yfinal(t-0.01),[0.01 0.02],LineWidth=2,Color="red"), grid on
fplot(Yfinal(t-0.02),[0.02 0.03],LineWidth=2,Color="red"), grid on
fplot(Yfinal(t-0.03),[0.03 0.04],LineWidth=2,Color="red"), grid on
hold off, legend("Tren de pulsos entrada", "Régimen estacionario NO senoidal resultante")

% *Nota: esta NO sería la forma de onda de tensión de salida a la resistencia de "carga" R_L en un rectificador, porque la intensidad absorbida por la fuente no puede ser negativa si tiene diodos (y aquí sí que lo será cuando V_c>u). Sigue siendo, por tanto, un "ejercicio académico" por el momento.
%
%