%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  inrush.m - Calculates the inrush current to an unloaded
%             transformer by solution of the nonlinear DEQ.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
f=60; R1=0.6; w=2*pi*f; % Frequency, primary resistance
Vm=sqrt(2)*240; zeta=0*pi/180; % Voltage, phase shift
ncyc=12;              % No. source cycles
% No load voltage(rms) and magnetizing current(rms) arrays
Voc=[0 200 240 250 260 270 280 290 300 310 320 480];
Im=[0 0.52 0.65 0.7 0.75 0.85 1 1.5 3 6 12 120]*sqrt(2);
Lam=Voc/4.44/f;      % Flux linkage(max value) array   
np=length(Im);
% Vmax=300;          % Activate to check magnetization curve
% for i=1:length(Voc) 
%    if Voc(i)>Vmax; nn=i-1; break; 
%    else; end
% end
% plot(Im(1:nn),Voc(1:nn)); grid;
% ylabel('Voltage'); xlabel('Magnetizing current'); pause
% plot(Im(1:nn),Lam(1:nn)); grid;
% ylabel('Flux linkage'); xlabel('Magnetizing current'); pause
T=1/f; t0=0; x0=0; x=[x0]; t=[t0];
h=10e-06; incr=T/100; tsav=t0;
for i=1:fix(ncyc*T/h) % Solution for magnetizing current(x)
   if x0 == 0; x0=1e-08; else; end
   if abs(x0) > max(Lam)   % Lam-Im slope
      dLamdi=(Lam(np)-Lam(np-1))/(Im(np)-Im(np-1))/sqrt(2);
   else
      a=interp1(Im,Lam,abs(x0)); b=interp1(Im,Lam,1.01*abs(x0));
      dLamdi=(b-a)/(0.01*abs(x0));
   end
   xx=x0; tt=t0;
   k1=h*(-R1*xx+Vm*sin(w*tt-zeta))/dLamdi;  % Fourth-order
   xx=x0+k1/2;tt=t0+h/2;                    % Runge-Kutta
   k2=h*(-R1*xx+Vm*sin(w*tt-zeta))/dLamdi;  % integration
   xx=x0+k2/2; tt=t0+h/2;                   % with fixed
   k3=h*(-R1*xx+Vm*sin(w*tt-zeta))/dLamdi;  % increment h
   xx=x0+k3; tt=t0+h;
   k4=h*(-R1*xx+Vm*sin(w*tt-zeta))/dLamdi;
   x0=x0+(k1+2*k2+2*k3+k4)/6;
   t0=t0+h;
   if (t0-tsav) >= incr
      tsav=t0; t=[t t0]; x=[x x0];
   else; end
end
plot(t,x); grid; title('Inrush current');
ylabel('Primary current, A'); xlabel('Time, s');
text(0.7*max(t),0.9*max(x),['Voltage = ',num2str(Vm/sqrt(2))]);
text(0.7*max(t),0.8*max(x),['Phase angle = ',num2str(zeta*180/pi)]);
text(0.7*max(t),0.7*max(x),['No. cycles = ',num2str(ncyc)]);