% lcl_dab_2igbt.m                                           	23-3-15

% Calculate and display (if plott = 1) voltages and currents in DABLCL.
% Between B1 and B2 is an LCL network of series L1 R1, shunt C1,  
%   a series L2 R2, shunt L3 to trfmr. tr = trfmr turns ratio, n1/n2.
% The magnetising reactance is XL3 = k4*xd*n. 
% The VB1 lead, and the modulation angles are thi*pi, m1*pi and m2*pi.
% N.B. Loss calcs assume forward operation, and 2 igbt bridges.
% Legs 1 and 3,(Q1,2) and (Q5,6) lead legs 2 and 4 (Q3,4) and (Q7,8).
Plott = 1;                                  % 1 => yy plots are done
tics = 800;                                 % number of points, columns
rows = 50;     n = (1:2:2*rows - 1)';     	% rows = odd freqencies
f = 50000;                                  % for plot axis t labels
temp = 75;          	% Enter IGBT (STGWA19NC60) junction temperature

% Enter the parameter values :
% The reactances of the network legs equal xb at the fs fundamental.
%	XL1 = k1*xd*n, XL2 = k2*xd*n, XC1 = k3*xd/n, XL3 = k4*xd*n.
k1 = 1.0; k2 = 1.0; k3 = 1.0; k4 = 1e3;
m1 = 0.68;                                  % modulation index, 0..1
m2 = 0.66;                                  % modulation index, 0..1
m1 = 1.0;
m2 = m1;
phi = 0.5;                                  % v1 lead index, +ve => fwd
vdc1 = 400;                             	% DC supply to Bridge 1
vdc2 = 400;                               	% DC supply to Bridge 2
tr = 1;
xd = 32.1;                                	% 41.48uH, 244.3nF, 10 kW
r1 = 0.0; r2 = 0.0;                         % i^2r is calcd after I

% k1 = 1.007; k2 = 0.997; k3 = 1.021; k4 = 22.7;   
% tr = 1.085; xd = 31.83; vdc1 = 400; vdc2 = 400;       
% r1 = [0.0297; 0.0923; 0.218; ones(rows - 3,1)];
% r2 = [0.0713; 0.148; 0.293; ones(rows - 3,1)];                     

% k1 = 1.0; k2 = 1.0; k3 = 1.0; k4 = 10; r1 = 0.1; r2 = 0.3;    % TEST
% vdc1 = 420; vdc2 = 360; thi = 0.35; m1 = 0.7; m2 = 0.8;       % TEST
% tr = 1.3; xd = 23;                                            % TEST 

% Calculate NETWORK voltage phasors of each harmonic, v1 = ref = 0 deg
v1n = (vdc1*4/pi)./n.*sin(n*pi*m1/2);
v2n = (vdc2*tr*4/pi)./n.*sin(n*pi*m2/2).*exp(-1i*n*phi*pi);

% Set the scales:
th = 0:2*pi/tics:(tics-1)/tics*2*pi;        % fundamental phase, 1 x tics
thn = n*th;                                 % harmonic phase, n x tics

% The approx. NETWORK solutions for r1, r2 = 0 are :
den = (k1*k3 + k2*k3 - k1*k2*n.^2);
i1n = 1i*(k2*v1n.*n.^2 - k3*v1n + k3*v2n)./(xd*n.*den);
i2n = -1i*(k1*v2n.*n.^2 + k3*v1n - k3*v2n)./(xd*n.*den);

ic1n = i1n - i2n;                             	% current in C1 shunt
imn = 1i*v2n./(k4*xd*n);                    % trfmr magnetising current
i4n = i2n + imn;                            % into ideal trfmr prim
ib2n = tr*i4n;
vc1n = -1i*ic1n*k3*xd./n;                     	% VC1
vl1n = 1i*i1n*k1*xd.*n;                        	% vL1 (IL1 = IR1)
vl2n = 1i*i2n*k2*xd.*n;                       	% vL2 (IL2 = IR2)
clear k1 k2 k3 k4 xb den rows

% Voltage and current values follow. Separate the phasor column vector
% into magnitude(m) and phase(p) vectors. Form duplicate matrices.
% Sum components at each point in time/angle => time function

v1m = abs(v1n); v1md = v1m*ones(size(th));              % bridge 1
v1p = angle(v1n); v1pd = v1p*ones(size(th)); 
v1 = sum(v1md.*cos(v1pd + thn));    
v1rms = sqrt(sum(v1m.^2)/2); v1max = max(v1); v1min = min(v1);
clear v1n v1m v1p v1md v1pd

i1m = abs(i1n); i1md = i1m*ones(size(th));            	% I1 (= IL1) 
i1p = angle(i1n); i1pd = i1p*ones(size(th));          
i1 = sum(i1md.*cos(i1pd + thn));     
i1rms = sqrt(sum(i1m.^2)/2); i1max = max(i1); i1min = min(i1);
clear i1n i1p i1md i1pd 

vl1m = abs(vl1n); vl1md = vl1m*ones(size(th));        	% vL1
vl1p = angle(vl1n); vl1pd = vl1p*ones(size(th));     
vl1 = sum(vl1md.*cos(vl1pd + thn));     
vl1rms = sqrt(sum(vl1m.^2)/2); vl1max = max(vl1); vl1min = min(vl1);
clear vl1n vl1m vl1p vl1md vl1pd

vl2m = abs(vl2n); vl2md = vl2m*ones(size(th));         	% vL2
vl2p = angle(vl2n); vl2pd = vl2p*ones(size(th));     
vl2 = sum(vl2md.*cos(vl2pd + thn));     
vl2rms = sqrt(sum(vl2m.^2)/2); vl2max = max(vl2); vl2min = min(vl2);
clear vl2n vl2m vl2p vl2md vl2pd

i2m = abs(i2n); i2md = i2m*ones(size(th));            	% IL2 = I2
i2p = angle(i2n); i2pd = i2p*ones(size(th));            
i2 = sum(i2md.*cos(i2pd + thn));      
i2rms = sqrt(sum(i2m.^2)/2); il2max = max(i2); il2min = min(i2);                         
clear i2n i2p i2md i2pd

ib2m = abs(ib2n); ib2md = ib2m*ones(size(th));      	% Ib2 from trfmr
ib2p = angle(ib2n); ib2pd = ib2p*ones(size(th));        %    secondary           
ib2 = sum(ib2md.*cos(ib2pd + thn));                 
ib2rms = sqrt(sum(ib2m.^2)/2); ib2max = max(ib2); ib2min = min(ib2);
clear ib2n ib2p ib2md ib2pd

vc1m = abs(vc1n); vc1md = vc1m*ones(size(th));        	% vC1  
vc1p = angle(vc1n); vc1pd = vc1p*ones(size(th));   
vc1 = sum(vc1md.*cos(vc1pd + thn));       
vc1rms = sqrt(sum(vc1m.^2)/2); vc1max = max(vc1); vc1min = min(vc1);
clear vc1n vc1m vc1p vc1md vc1pd

ic1m = abs(ic1n); ic1md = ic1m*ones(size(th));        	% IC1  
ic1p = angle(ic1n); ic1pd = ic1p*ones(size(th));   
ic1 = sum(ic1md.*cos(ic1pd + thn));       
ic1rms = sqrt(sum(ic1m.^2)/2); ic1max = max(ic1); ic1min = min(ic1);
clear ic1n ic1m ic1p ic1md ic1pd

vb2m = abs(v2n)/tr; vb2md = vb2m*ones(size(th));     	% RHS of trfmr
vb2p = angle(v2n); vb2pd = vb2p*ones(size(th));      	%   bridge 2
vb2 = sum(vb2md.*cos(vb2pd + thn));                 	
vb2rms = sqrt(sum(vb2m.^2)/2); vb2max = max(vb2); vb2min = min(vb2);
clear v2n vb2m vb2p vb2md vb2pd

p = i1.*v1;                         % instantaneous pwr/phase for plot
p0 = sum(p,2)/tics;              	% average pwr in    (zero loss)

% figure (2);
% plot(th,p)
% ylabel('Instantaneous Power (W)','fontsize',12)
% xlabel('Power Cycle Angle (rad)','fontsize',12)

disp('The bridge values V1, VB2, I1 and IB2 are :')
ac(1,:) = [v1rms v1max v1min];
ac(2,:) = [vb2rms vb2max vb2min];
ac(3,:) = [i1rms i1max i1min];
ac(4,:) = [ib2rms ib2max ib2min];
ac(1:4,:)

disp('The vl1, il1, ic1, vc1, vl2 and il2 rms, max and min values are')
ac(5,:) = [vl1rms vl1max vl1min];
ac(6,:) = [i1rms i1max i1min];
ac(7,:) = [ic1rms ic1max ic1min];
ac(8,:) = [vc1rms vc1max vc1min];
ac(9,:) = [vl2rms vl2max vl2min];
ac(10,:) = [i2rms il2max il2min];
ac(5:10,:)

if Plott == 1                 	% requested plot o/p
ext = uint16(tics/2);                               % extention, 0.5 cycle
t = 0:1e6/f/tics:(tics*1.5-1)*1e6/f/tics;     % time in usecs
v1d = [v1 v1(1:ext)];
i1d = [i1 i1(1:ext)];
mid1 = [vc1 vc1(1:ext)];      % SELECT DISPLAY FOR MIDDLE AXES (Plot 2),
mid2 = [ic1 ic1(1:ext)];    	%   from vl1t, vl2t, vc1t, vc2t
vb2d = [vb2 vb2(1:ext)];
ib2d = [ib2 ib2(1:ext)];
pd = [p p(1:ext)];

figure(1);
subplot(2,1,1)                      	% V1, I1
[AX,H1,H2] = plotyy(t,v1d,t,i1d);
set(get(AX(1),'Ylabel'),'String','v_{1}(V)','FontSize',14)
set(get(AX(2),'Ylabel'),'String','i_{1}(A)','FontSize',14)
set(H1,'LineStyle','-','LineWidth',2.0); 
set(H2,'LineStyle','-','LineWidth',2.0);
set(AX(1),'XTickLabel',[]);             % clear the x axes
set(AX(2),'XTickLabel',[]);
set(AX(1),'Ylim',[-500 500]);
yends = get(AX(1),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(1),'YTick',yends(1):yinc:yends(2))
set(AX(1),'YTickLabel',get(AX(1),'YTickLabel'),'FontSize',14)
set(AX(2),'Ylim',[-20 20]);
yends = get(AX(2),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(2),'YTick',yends(1):yinc:yends(2))
set(AX(2),'YTickLabel',get(AX(2),'YTickLabel'),'FontSize',14)
hleg = legend('v_1','i_1');
set (hleg,'FontSize',14,'location','southeast')
grid on;
subplot(2,1,2)                              % H handles are for line
[AX,H1,H2] = plotyy(t,vb2d,t,ib2d);           % AX handles are for axis
set(get(AX(1),'Ylabel'),'String','v_{2}(V)','FontSize',14)
set(get(AX(2),'Ylabel'),'String','i_{2}(A)','FontSize',14)
set(H1,'LineStyle','-','LineWidth',2)  	
set(H2,'LineStyle','-','LineWidth',2)  
xlabel('Time (\mus)','FontSize',14)             % "\mu" => greek "micro"
set(AX(1),'Ylim',[-500 500]);
yends = get(AX(1),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(1),'YTick',yends(1):yinc:yends(2))
set(AX(1),'YTickLabel',get(AX(1),'YTickLabel'),'FontSize',14)
set(AX(2),'Ylim',[-20 20]);
yends = get(AX(2),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(2),'YTick',yends(1):yinc:yends(2))
set(AX(2),'YTickLabel',get(AX(2),'YTickLabel'),'FontSize',14)
hleg = legend('v_{2}','i_{2}');
set (hleg,'FontSize',14,'location','southeast')
grid on;
clear ext t i1d iB2d pd mid1 mid2 v1d vB2d Plott hleg AX H1 H2 yends yinc
end

% Calculate the switching indices. The switch sequences are 14, 13, 23, 
%  24, 14 ... and 58, 57, 67, 68, 58 ... Leg 1 leads leg 2 and leg 3
%  leads leg 4. If m1 > 0, 0 < a,b < tics. 
a = uint16(tics*(2 - m1)/4 + 1);                        % S2 => on
b = uint16(tics*(2 + m1)/4 + 1);                        % S4 => on
c = uint16(tics*(2 + m2 + 2*phi)/4 + 1);                % S8 => on
d = uint16(tics*(6 - m2 + 2*phi)/4 + 1);                % S6 => on
if d > tics
    d = d - tics;
    if c > tics
        c = c - tics;                               % fold indices into
    end                                             %   range 1..tics
end

% Calculate the transistor currents. Each has a tics/2 conduction span.
iq2 = zeros(1,tics);
iq2(a:a+tics/2-1) = -i1(a:a+tics/2-1);        	% 1 < st,sp <= tics
iq4 = zeros(1,tics);                            %  
iq4(b:tics) = i1(b:tics);                       % conducts to end of cycle
iq4(1:b-tics/2-1) = i1(1:b-tics/2-1);       	% 	then finishes in next
iin = iq2 + iq4;                                % current from VDC1 to B1
sw1rms = sqrt(sum((iq2').^2)/tics);             % leg 1 rms current
sw2rms = sqrt(sum((iq4').^2)/tics);             % leg 2 rms current
% figure(3);  plot(th, iq2)                  	% diagnostic
% figure(4);  plot(th, iq4)                  	% diagnostic
% figure(9);  plot(th, iin)                   	% diagnostic
% x=1:tics*3/2;
% y=[i1 i1(1:tics/2)];
% z=[iq2 iq2(1:tics/2)];
% figure(10);  plot(x/tics, -y)                   	
% figure(9);  plot(x/tics, z)                   	

iq6 = zeros(1,tics);
if d > (tics/2 + 1)                            	% for start > stop
    iq6(d:tics) = ib2(d:tics);                  % conducts st..end then
    iq6(1:d-tics/2-1) = ib2(1:d-tics/2-1);   	%   beginning of next
else
    iq6(d:d+tics/2-1) = ib2(d:d+tics/2-1);     	% conducts st..sp
end
iq8 = zeros(1,tics);
if c > (tics/2 + 1)                            	% for start > stop
    iq8(c:tics) = -ib2(c:tics);               	% conducts st..end then
    iq8(1:c-tics/2-1) = -ib2(1:c-tics/2-1);   	%   beginning of next
else
    iq8(c:c+tics/2-1) = -ib2(c:c+tics/2-1);   	% conducts st..sp
end
iout = -iq6 - iq8;                              % current from B2 to VDC2 
sw3rms = sqrt(sum((iq6').^2)/tics);             % leg 3 rms current
sw4rms = sqrt(sum((iq8').^2)/tics);             % leg 4 rms current
% figure(7);  plot(th, iq6)                  	% diagnostic
% figure(8);  plot(th, iq8)                  	% diagnostic
% figure(9);  plot(th, iout)                  	% diagnostic

iin_avg = sum(iin)/tics;
iin_ms = sum(iin.^2)/tics;                      % mean square
iin_rms = sqrt(iin_ms);
iin_rms_ac = sqrt(iin_ms - iin_avg^2);          % rms ripple current
iin_max = max(iin); iin_min = min(iin);

iout_avg = sum(iout)/tics;
iout_ms = sum(iout.^2)/tics;                  	% mean square
iout_rms = sqrt(iout_ms);
iout_rms_ac = sqrt(iout_ms - iout_avg^2);     	% rms ripple current
iout_max = max(iout); iout_min = min(iout);

% Calculate the IGBT turn-on currents and therefore the switching  
%   energy losses, in uJ. If Isw = -ve prior to turn-on then that 
%   Sw has soft switch on and hard switch off, and vica versa. 
i2s = iq2(a);
[eoff,eon] = STGWA19N60_swg_loss(temp,vdc1,abs(i2s));
if i2s < 0
    Eon2 = 0;                           	% ZVS on, E approx = 0
    Eoff2 = eoff;
else
    Eon2 = eon;
    Eoff2 = 0;                            	% ZCS off, E approx = 0
end
i4s = iq4(b);
[eoff,eon] = STGWA19N60_swg_loss(temp,vdc1,abs(i4s));
if i4s < 0
    Eon4 = 0;                           	% ZVS on, E approx = 0
    Eoff4 = eoff;
else
    Eon4 = eon;
    Eoff4 = 0;                            	% ZCS off, E approx = 0
end
i6s = iq6(d);
[eoff,eon] = STGWA19N60_swg_loss(temp,vdc2,abs(i6s));
if i6s < 0
    Eon6 = 0;                           	% ZVS on, E approx = 0
    Eoff6 = eoff;
else
    Eon6 = eon;
    Eoff6 = 0;                            	% ZCS off, E approx = 0
end
i8s = iq8(c);
[eoff,eon] = STGWA19N60_swg_loss(temp,vdc2,abs(i8s));
if i8s < 0
    Eon8 = 0;                           	% ZVS on, E approx = 0
    Eoff8 = eoff;
else
    Eon8 = eon;
    Eoff8 = 0;                            	% ZCS off, E approx = 0
end
Esw1 = 2*(Eon2 + Eoff2 + Eon4 + Eoff4); 	%
Esw2 = 2*(Eon6 + Eoff6 + Eon8 + Eoff8);     % 	
Psw1 = Esw1*f*1e-6;                      % B1 switching power in W
Psw2 = Esw2*f*1e-6;                      % B2 switching power in W
Psw = Psw1 + Psw2                       	% total switching power

% Calculate the IGBT conduction losses. If the IGBT has -ve I its body
%   diode conducts. Otherwise its transistor VCE characteristic applies.
%   Switches 2,4,6,8 start their conduction at a,b,d,c and span tics/2.

ic = iq2(a:a+tics/2-1);                     % conduction current
icd = -ic(ic < 0); ict = ic(ic >= 0);       % diode, transistor current
pd = icd.*STGWA19N60_vd(temp,icd);          % diode v
pt = ict.*STGWA19N60_vt(temp,ict);      	% transistor v
pc2 = sum([pd pt]);

ic = zeros(1,tics/2);                       % half cycle
ic(1:tics-b+1) = iq4(b:tics);            	% conducts st..end and then
ic(tics-b+2:tics/2) = iq4(1:b-tics/2-1);	%   next cycle beginning
icd = -ic(ic < 0); ict = ic(ic >= 0);       % diode, transistor current
pd = icd.*STGWA19N60_vd(temp,icd);          % diode v
pt = ict.*STGWA19N60_vt(temp,ict);      	% transistor v
pc4 = sum([pd pt]);

if d > (tics/2 + 1)                       	% for start > stop
    ic(1:tics-d+1) = iq6(d:tics);           % conducts st..end, then
    ic(tics-d+2:tics/2) = iq6(1:d-tics/2-1);%   next cycle beginning
else
    ic(1:tics/2) = iq6(d:d+tics/2-1);       % conducts st..sp
end
icd = -ic(ic < 0); ict = ic(ic >= 0);       % diode, transistor current
pd = icd.*STGWA19N60_vd(temp,icd);          % diode v
pt = ict.*STGWA19N60_vt(temp,ict);      	% transistor v
pc6 = sum([pd pt]);

if c > (tics/2 + 1)                       	% for start > stop
    ic(1:tics-c+1) = iq8(c:tics);         	% conducts st..end, then
    ic(tics-c+2:tics/2) = iq8(1:c-tics/2-1);%   next cycle beginning
else
    ic(1:tics/2) = iq8(c:c+tics/2-1);       % conducts st..sp
end
icd = -ic(ic < 0); ict = ic(ic >= 0);       % diode, transistor current
pd = icd.*STGWA19N60_vd(temp,icd);          % diode v
pt = ict.*STGWA19N60_vt(temp,ict);      	% transistor v
pc8 = sum([pd pt]);

Pcon1 = 2*(pc2 + pc4)/tics;                 % B1 conduction losses
Pcon2 = 2*(pc6 + pc8)/tics;                 % B2 conduction losses
Pcon = Pcon1 + Pcon2                        % total conduction losses

disp('Bridge legs 1..4 (down) have switch-on, maximum, minimum and rms')
disp('device currents of :')
sw(1,:) = [iq2(a) max(iq2) min(iq2) sw1rms];
sw(2,:) = [iq4(b) max(iq4) min(iq4) sw2rms];  
sw(3,:) = [iq6(d) max(iq6) min(iq6) sw3rms];
sw(4,:) = [iq8(c) max(iq8) min(iq8) sw4rms];
sw

disp('The avg, rms, ac_rms, max and min values (across) for the bridge')
disp(' input and output dc currents are :')
dc(1,:) = [iin_avg iin_rms iin_rms_ac iin_max iin_min];
dc(2,:) = [iout_avg iout_rms iout_rms_ac iout_max iout_min];
dc

% r1 and r2 are resistances of L1 and L2/transformer
Pr1 = sum(i1m.^2.*r1)/2; Pr2 = sum(i2m.^2.*r2)/2; Pr = Pr1 + Pr2 
% Magnetising losses can't be modelled as a shunt R.
Pmag = (m2^1.4)*67                         % 67 W at m2 = 1
% Pmag = 0;
p0
p1 = p0 + Pr2 + Pcon2 + Psw1 + Pmag/2    	% adjust p1 and p2 for
p2 = p0 - Pr1 - Pcon1 - Psw2 - Pmag/2     	%  other losses 
loss = p1 - p2                              % total losses

% Calculate the volt.sec products
vsl1 = zeros(1,tics); vsl2 = vsl1; vsm1 = vsl1; vsm2 = vsl1;	
vsl1(1) = vl1(1); vsl2(1) = vl2(1); 
vsm1(1) = vc1(1); 
for i = 2:tics
  	vsl1(i) = vsl1(i-1) + vl1(i);           % L1
  	vsl2(i) = vsl2(i-1) + vl2(i);           % L2
    vsm1(i) = vsm1(i-1) + vc1(i);           % trfmr prim and integ L2
end
fact = 1e6/f/tics;                      	% scale with period in us
vsl1 = (max(vsl1) - min(vsl1))*fact;
vsl2 = (max(vsl2) - min(vsl2))*fact;
vsm1 = (max(vsm1) - min(vsm1))*fact;
vsm2 = vdc2*m2*5e5/f;
apl1 = i1rms*vsl1; apl2 = i2rms*vsl2; 
apm1 = i2rms*vsm1; apm2 = i2rms*vsm2*tr; apm3 = ib2rms*vsm2; 
apdisc = apl1 + apl2 + apm2 + apm3;
apinteg = apl1 + apm1+ apm3;
disp('The V.us products of L1, L2, Primary+L2, and Secondary are :')
[vsl1 vsl2 vsm1 vsm2]
disp('The V.us.Irms Products for L1, L2, Primary+L2 and Secondary are :')
[apl1 apl2 apm1 apm3]
disp('The V.us.Irms Discrete and Integ L2 Products are :')
[apdisc apinteg]
p2
Eff = p2/p1

% clear n p e f g h iin_ms iout_ms ans freq m1 m2 thi tics % a b c d 

% Test:    With tics = 12000, rws = 562, freq = 50k, temp = 75,
%           ri = 0.1, r2 = 0.3, k1..3 = 1, k4 = 10, STGWA19N60
% V1,2 = 420,360, thi = 0.35, m1,2 = 0.7,0.8, tr = 1.3, xd = 23

%  => p0 = 5230.6, p1 = 5442.5, p2 = 5110.8, 
% i1rms = 17.4345, i2rms = 14.7131, ib2rms = 18.2323, 
% Psw = 42.1, Pcon = 188.1, Pr = 95.3, Pmag = 6.3, Ptot = 331.7 
% vsl1 = 3702.7   vsl2 = 3305.6   vsm1 = 5664.7   vsm2 = 2882.2   


%   PLECS tests with LCL_Loss_TestDAB.slx, using IGBTD_STGWA19N60
%   => All results very similar except P1,2 






