% lcl_loss_2igbt.m                                             	6-6-15

% Called by lcl_loss_table_2igbt.m

% Calculate LCL DAB converter losses.
% 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 reactances of the network legs at the fs fundamental are:
%	XL1 = k1*xd*n, XL2 = k2*xd*n, XC1 = k3*xd/n, 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).
% The loss calculations assume that each bridge uses STGWA19N60 igbts. 

function [r] = lcl_loss_2igbt(phi,m1,m2)
tics = 800;                                 % number of points, columns
rws = 55;     n = (1:2:2*rws - 1)';     	% rows = odd freqencies

% Enter the parameter values for the 4 kW converter:
k1 = 1.007; k2 = 0.997; k3 = 1.021; k4 = 22.7;   
tr = 1.085; xd = 31.83; vdc1 = 400; vdc2 = 400;       
f = 5e4; temp = 75; 

% Use the following for resistance(f) as a result of skin effect:
r1 = ones(rws - 3,1);                     % on LHS of network
% append 6 values for f1..11*f1, = R1 transf + RC1 (say 0.01 const)
r1 = [0.0297; 0.0923; 0.218; r1];   	
r2 = ones(rws - 3,1);                     % on RHS (B2 side) of trfmr
% append 6 values for f1..11*f1, = R2 transf + RC2 (say 0.01 const)
r2 = [0.0713; 0.148; 0.293; r2];                      	
% r1 = 0; r2 = 0;

% 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):
den = (k1*k3 + k2*k3 - k1*k2*n.^2);             % denominator
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);

imn = 1i*v2n./(k4*xd*n);                    % trfmr magnetising current
i4n = i2n + imn;                            % into ideal trfmr prim
    % N.B. ib2 doesn't equal tr * i2
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));    
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);
clear i1n i1p i1md i1pd 

i2m = abs(i2n);                                     % IL2 = I2 into trfmr 
clear i2n                                           %  with L3 in //

i4m = abs(i4n); i4md = i4m*ones(size(th));       	% I5 into ideal
i4p = angle(i4n); i4pd = i4p*ones(size(th));        %   trfmr primary           
ib2 = tr*sum(i4md.*cos(i4pd + thn));                % IB2 is out of 
ib2rms = tr*sqrt(sum(i4m.^2)/2);                    %  trfmr secondary 
clear i4n i4m i4p i4md i4pd                         

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

% 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

% 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,vdc1,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
Es1 = 2*(Eon2 + Eoff2 + Eon4 + Eoff4); 	%
Es2 = 2*(Eon6 + Eoff6 + Eon8 + Eoff8);	% 	
Psw1 = Es1*f/1e6;                       % B1 switching power in W
Psw2 = Es2*f/1e6;                       % B2 switching power in W

% 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]);

pc2 = pc2/tics;  pc4 = pc4/tics;            % now in W
pc6 = pc6/tics;  pc8 = pc8/tics;            % now in W
Pcon1 = 2*(pc2 + pc4);                      % conduction losses, W
Pcon2 = 2*(pc6 + pc8);                      %  

Ptot1 = Psw1 + Pcon1; Ptot2 = Psw2 + Pcon2;

% 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
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

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

% Gather the output data into the return column vector:
r = [phi; m1; m2; Eon2; Eoff2; Eon4; Eoff4; pc2; pc4; Psw1; Pcon1; ...
    Ptot1; Eon6; Eoff6; Eon8; Eoff8; pc6; pc8; Psw2; Pcon2; ...
    Ptot2; Pmag; Pr; loss; p1; p2; Eff; i1rms; ib2rms]; 	% 29 rows









