% lcl_d_2igbt.m                                               24-9-15

% Calculates the i/p and o/p currents, the integrated volt.us 
%   products and the i/p and o/p powers for LCL DAB with L_leakage split.
% Between B1 and B2 is an LCL network of series L1 R1, shunt C1,  
%  series L2 R2, shunt L3, series L4 to trfmr. tr = trfmr ratio, n1/n2.
% XL1, XL2, XC1, XL3, XL4 proportions of Xd are k1, k2/2, k3, k4, k2/2.
% The magnetising inductance L3 = 2km*L2 = 2km*L4. Ideally (tuned) => 
%   k1 = 1, k2 = (2 + 4km)/(1 + 4km), k3 = 1, k4 = km*k2 
%   k4 = km(2 + 4km)/(1 + 4km), tr = (vdc1/vdc2)*(1 + 1/(2*km))
% The VB1 lead, and the modulation angles are phi*pi, m1*pi and m2*pi.
% Legs 1 and 3,(Q1,2) and (Q5,6) lead legs 2 and 4 (Q3,4) and (Q7,8).
% N.B. i3 is current in L2 (lumped LCL has i3 = iC1).

% Variant _d is for LCL DAB with L_leakage distributed about L_magnetising

function [p1,p2,i1rms,i3rms,ib2rms,ic1rms,vsl1,vsm1,vsm2] = ...
    lcl_d_2igbt(vdc1,vdc2,phi,m1,m2,xd,km,trr,f)

tics = 800;                                     % columns = time points
rws = 50;     n = (1:2:2*rws - 1)';              % rows = odd freqencies
temp = 75;          	% Enter IGBT (STGWA19NC60) junction temperature 

k1 = 1; k2 = (2 + 4*km)/(1 + 4*km); k3 = 1; k4 = km*k2; 
% N.B. TRANSFORMER RATIO used (tr) IS WEIGHTED:
tr = trr*(1 + 1/(2*km));                      	% tr = turns ratio n1/n2
% k1 = 1; k2 = 1.111; k3 = 1; k4 = 2.222;    	% For km = 2  tr = 1.25
% k1 = 1; k2 = 1.024; k3 = 1; k4 = 10.24;      	% For km = 10  tr = 1.05
% k1 = 1; k2 = 1; k3 = 1; k4 = 1e4;             % For km = 1e4  tr = 1

% 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

% Solution for r1, r2 = 0 (from lcl_netsolv_gen):
nsq = n.^2;
den = 1i*n*xd.*(k2^2*k3 + 2*k1*k2*k3 + 4*k1*k3*k4 ...
    + 4*k2*k3*k4 - n.^2*k1*k2*(4*k4 + k2));
i1n = (k3*v1n*(2*k2 + 4*k4) - nsq.*v1n*k2*(k2 + 4*k4) - 4*k3*k4*v2n)./den;
i2n = -(2*((2*k1 + k2 + 2*k4)*k3*v2n - 2*k3*k4*v1n  ...
    - (k1*k2 + 2*k1*k4)*nsq.*v2n))./den;
i3n = (2*(2*k1*k4*v2n.*nsq + k2*k3*v1n + 2*k3*k4*v1n - 2*k3*k4*v2n))./den; 
ic1n = i1n - i3n;                          	% current in C1 shunt
vc1n = -1i*ic1n*k3*xd./n;                 	% VC1 = vsm1
vl1n = 1i*i1n*k1*xd.*n;                  	% vL1
% vl2n = 1i*i3n*k2*xd.*n/2;              	% vL2
% vl3n = 1i*(i3n - i2n)*k4*xd.*n;          	% vL3

% 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

i1m = abs(i1n); i1md = i1m*ones(1,tics);          	% I1 (= IL1) 
i1p = angle(i1n); i1pd = i1p*ones(1,tics);          
i1 = sum(i1md.*cos(i1pd + thn));     
i1rms = sqrt(sum(i1m.^2)/2);

v1m = abs(v1n); v1md = v1m*ones(1,tics);        	% bridge 1
v1 = sum(v1md.*cos(thn));                           % REFERENCE, 0 deg

% i3m = abs(i3n); i3md = i3m*ones(size(th));    	% I3 = IL2 = IR2 
% i3p = angle(i3n); i3pd = i3p*ones(size(th));          
% i3 = sum(i3md.*cos(i3pd + thn));     
i3rms = sqrt(sum((abs(i3n)).^2)/2);

vl1m = abs(vl1n); vl1md = vl1m*ones(1,tics);     	% vL1
vl1p = angle(vl1n); vl1pd = vl1p*ones(1,tics);     
vl1 = sum(vl1md.*cos(vl1pd + thn));     

ic1m = abs(ic1n);                                   % IC1, needed for  
ic1rms = sqrt(sum(ic1m.^2)/2);                   	%   C1 costing

i2m = abs(i2n); i2md = i2m*ones(1,tics);	     	% into prim of trfmr 
i2p = angle(i2n); i2pd = i2p*ones(1,tics);          
ib2 = tr*sum(i2md.*cos(i2pd + thn));                % at bridge 2    
ib2rms = tr*sqrt(sum(i2m.^2)/2);                

vc1m = abs(vc1n); vc1md = vc1m*ones(1,tics);       	% vC1  
vc1p = angle(vc1n); vc1pd = vc1p*ones(1,tics);   
vc1 = sum(vc1md.*cos(vc1pd + thn));       

% v2m = abs(v2n); v2md = v2m*ones(size(th));    	% RHS of network
% v2p = angle(v2n); v2pd = v2p*ones(size(th));    	%   = vm2
% v2 = sum(v2md.*cos(v2pd + thn));                        

% p0 = sum(i1.*v1,2)/tics;                        % average pwr throughput 

% 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

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

% The following calcs are for the LOSSLESS system:
iin = iq2 + iq4;                                % current from VDC1 to B1
iin_avg = sum(iin)/tics;                       	% dc input current
p0 = vdc1*iin_avg;                            	% average pwr throughput 

% 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/1e6;                          % B1 switching power
Psw2 = Esw2*f/1e6;                          % B2 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.
% In all cases the transistor conducts for a half cycle, tics/2 points.

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

% Calculate the integrated volt.sec products: 
vsl1 = zeros(1,tics); vsm1 = vsl1; 	
vsl1(1) = vl1(1); vsm1(1) = vc1(1); 
for i = 2:tics
  	vsl1(i) = vsl1(i-1) + vl1(i);           % L1
    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;
vsm1 = (max(vsm1) - min(vsm1))*fact;
vsm2 = m2*vdc2*5e5/f;                    	% trfmr secondary
% apl1 = vsl1*i1rms;                     	% area product for L1
% aptr = vsm1*i3rms + vsm2*ib2rms;      	% area product for trfmr

% Pr1 = i1rms^2*r1;                           % resistive losses 
% Pr2 = i3rms^2*r2;                           % resistive losses
% Pmag = (m2^1.4)*67                         % 67 W at m2 = 1
% p1 = p0 + Pr2 + Pcon2 + Psw1 + Pmag/2;     	% adjust p1 and p2 for
% p2 = p0 - Pr1 - Pcon1 - Psw2 - Pmag/2;     	%  other losses 
p1 = p0 + Pcon2 + Psw1;             	% adjust p1 and p2 for swg
p2 = p0 - Pcon1 - Psw2;              	%  and cond losses

% LOSSLESS CALCULATIONS. Have iin_avg. Do same for iout:
iout = -iq6 - iq8;                        	% current to VDC1 from B2
iout_avg = sum(iout)/tics;                 	% dc input current
% Use p1 and p2, a better approx of iin and iout can be made:
iin_avg = iin_avg*p1/p0;                 	% correct current with ratio    
iout_avg = iout_avg*p2/p0;              	%   of actual to lossless p
% Calc ripple - the AC component remains the same, from AC link
iin_ms = sum(iin.^2)/tics;                	% mean square
iripl1 = sqrt(iin_ms - iin_avg^2);       	% rms input ripple current
iout_ms = sum(iout.^2)/tics;            	% mean square
iripl2 = sqrt(iout_ms - iout_avg^2);       	% rms output ripple current

if phi < 0  Eff = p1/p2;    else    Eff = p2/p1;    end
end

% Assume magnetising core loss is prop. to B^2.8 * f^1.46 (from 3C90).
% lcl_opt_2igbt_m(400,400,0.5,1,1,32,2,1.25,5e4,80) =>
% Pmag trfmr = 27 W, Pmag L1 = 13 W, Pr1 = 13 W, Pr2 = 27 W
% P1 = 4175.3, P2 = 3890.5, i1 = 11.3595, i3 = 11.3595, ib2 = 13.1168
% vsl1 = 3304.7, vsm1 = 4663.6


