%% Applications of ODE: Inrush overcurrent in a linear inductor (I)
% © 2025, Antonio Sala, Universitat Politècnica de València. All rights reserved.
%
% Presentation in video: https://personales.upv.es/asala/YT/V/LinrushlinEN.html
% *Python code there as well.
% See also https://youtu.be/x24545H5HdU?si=8KGq9tnQAULMkciZ  from other channel.

L=1; R=1.5;

% Differential equation:
% \frac{di}{dt}=\frac{1}{L}\cdot (V-R\cdot i)
% Normalized form \dot x=Ax+Bu:
% \frac{di}{dt}=\frac{-R}{L}\cdot i+ \frac{1}{L}\cdot V

t_final=2.2;
syms t real
syms i(t)
omega=50*2*sym(pi); %50 Hz AC
ph0=0*pi/2; %initial phase (we'll play with it later on)
%ph0=-phaseI_ss;
V_peak=120;
V(t)=V_peak*sin(omega*t+ph0) %input voltage
vpa(V(t),4)

% Let us solve \frac{di}{dt}=\frac{1}{L}\cdot (V-R\cdot i):

sol=dsolve( diff(i)==(V(t)-R*i(t))/L , i(0)==0) %ODE symbolic solution
vpa(sol,4)
circuit_pole=-R/L %coefficient of the exponential
t_settling=fzero(@(t) exp(circuit_pole*t)-0.05,0) % 5 percent settling time
fplot(sol,[0 t_final],LineWIdth=1.5), grid on, yline(0), xline(t_settling)

%% No resistor particular case
% In this case \frac{di}{dt}=\frac{1}{L}V; hence i(t)=i(0)+\frac{1}{L}\cdot \int_0^tV(\tau)\,d\tau

syms tau real
i_noR=int(V(tau)/L,0,t)+0;
vpa(i_noR,4)
figure()
fplot(sol,[0 t_final/2], LineWidth=1.5), hold on
fplot(i_noR,[0 t_final/2]), hold off
grid on, yline(0)

%% Particular case 2: steady state sinusoidal regime

Z=R+L*omega*1j %impedance
amplI_ss=abs(1/Z)*V_peak; %Current amplitude
eval(amplI_ss)
phaseI_ss=angle(1/Z);
eval(phaseI_ss)*180/pi
I_ss=simplify(amplI_ss*sin(omega*t+phaseI_ss+ph0),120);
vpa(I_ss,4)
figure()
fplot(sol,[0 t_final],LineWidth=1.5), hold on
fplot(I_ss,[0 t_final]), hold off
yline([-2 -1 0 1 2]'*eval(amplI_ss),':'), grid on
legend("transient", "steady state", Location="best"),xline(t_settling)

%% Conclusions
% Initial response close to "no resistor" case, which has a non-vanishing DC component.
% \bullet If V(0)=0, we have twice initial current to that of steady state from impedance computations (2x inrush current); it takes \approx 3\frac{L}{R} seconds to vanish.
% \bullet If phase of voltage equal to that of conductance (inverse impedance) then we start "spot on" the steady-state sinusoid, inrush current avoided. It is close to 90º for a heavily inductive load.
% \bullet Inrush current in real transformers may be much higher (6x or more), due to magnetic saturation (in another video). Harmonics (non-sinusoidal waveforms) also appear in such a case. Fuse ratings, energization timing, serial resistance, etc. must be considered with care in, say, high-power transformer startup.