%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    invRL.m   -  Calculates load phase currents for
%                 six-step inverter with R-L  if
%                 load current continuous.
%                 ( LSIM routine not required )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; figure(1);
ncyc=3; VB=300; R=3; L=0.001; freq=60; % ncyc - cycles of solution 
npts = 300; t = linspace(0, 1/freq, npts ); % Saves only last cycle 
x0 = [0 0 0]; k = length(t);       % Set ICs for first cycle
v13 = ones(1,npts/6)*VB/3/L;       % Form phase voltages
v23 = ones(1,npts/6)*2*VB/3/L;
van = [v23 v13 -v13 -v23 -v13 v13];
vbn = [-v13 v13 v23 v13 -v13 -v23];
vcn = [-v13 -v23 -v13 v13 v23 v13];
a = -eye(3)*R/L;    %  Build matrices for xdot = ax + bu
b = eye(3);
u = [van' vbn' vcn'];
c = eye(3); d = zeros(3);
for m=1:ncyc
clear x y;           % Implement solution for vector DEQ below
n = length(a); [P,D] = eig(a); PI = inv(P);
x = x0;
for i=2:length(t)
y = zeros(n,1); exp1 = zeros(n); exp2 = zeros(n);
for j=1:n
exp1(j,j) = exp( D(j,j) .*( t(i)-t(i-1) ) );
exp2(j,j) = ( exp( D(j,j) .*( t(i)-t(i-1) ) ) - 1 )/D(j,j);
end
y = P*exp1*PI*x0' + P*exp2*PI*b*( u(i,:)+u(i-1,:) )'/2;
x = [x;y']; x0 = y';
end
end
subplot(211), plot(0,VB,0,-VB,t,van*L,t,vbn*L,'--',t,vcn*L,'-.'); 
grid; title('LOAD PHASE VOLTAGES');
ylabel('VOLTAGE - V'); xlabel('TIME -s');
A=[1:length(t)]';
subplot(212),plot(t,x(A,1),t,x(A,2),'--',t,x(A,3),'-.');grid; 
title('R-L WYE LOAD FED BY 6-STEP INVERTER');
ylabel('LOAD PHASE CURRENTS - A'); xlabel('TIME - s');
text( 0.15, 0.91, ['R = ',num2str(R),'  L = ',num2str(L), ...
'  VB = ',num2str(VB)], 'sc');