% File  clc_loss_2igbt.m    CLC converter losses.           	6-6-15

% Called by clc_loss_table_2igbt.m

function [r] = clc_loss_2igbt(phi,m1,m2)
% For analysis purposes, between B1 and B2 is a CLC network of series 
%   L1 R1 C1, shunt L2, a series C2_ R2_ to trfmr to B2. The actual 
%   C2 and R2 are between trfmr RHS and B2. The network values are
%   C2_ = C2/tr^2, R2_ = R2*tr^2.  
% The reactance values at fs of L1, C1, L2 and C2' are k1*xb, k2*xb,
%   k3*xb and k4*xb respectively.  
% The VB1 lead, and the modulation angles are thi*pi, m1*pi and m2*pi.
% The loss calculations assume that each bridge uses igbts. 

tics = 800;                                 % number of points, columns
rows = 56;     n = (1:2:2*rows - 1)';     	% rows = odd freqencies
temp = 75; 

% Actual parameters for Conv 2,  approx. constant for f variation
%   (ideal 101.3 uH, 50 nF, 101.3 uH, 100 nF, tr = 1)
k1 = 0.993; k2 = 2.049; k3 = 0.997; k4 = 0.965;  % approx. f constant
tr = 1.01; xd = 31.83; f = 5e4; vdc1 = 400; vdc2 = 400; 
r1 = 0.133; r2 = 0.126;                     % r2 on B2 side of trfmr

% k1 = 1; k2 = 2; k3 = 1; k4 = 1;  
% tr = 1; xd = 31.83; f = 5e4; vdc1 = 400; vdc2 = 400; 
% 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

% The solutions for r1 = r2 = 0 are :
den = j*xd*(n.^2*(k1*k4 + k2*k3 + k3*k4) - k2*k4  - k1*k3*n.^4);
i1n = -((k3*v1n - k3*v2n).*n.^3 - k4*v1n.*n)./den;
i2n = ((k1*v2n - k3*v1n + k3*v2n).*n.^3 - k2*v2n.*n)./den;

% 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 => result.

v1m = abs(v1n); v1md = v1m*ones(size(th));       	% bridge 1, line
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 (= IC1, IL1) 
i1p = angle(i1n); i1pd = i1p*ones(size(th));          
i1 = sum(i1md.*cos(i1pd + thn));     
i1rms = sqrt(sum(i1m.^2)/2);
clear i1p i1md i1pd 

il2n = i1n - i2n;                                   % current in L2 shunt
il2m = abs(il2n);     

i2m = abs(i2n); i2md = i2m*ones(size(th));        	% start on trfmr LHS
i2p = angle(i2n); i2pd = i2p*ones(size(th));     	%   but then RHS 
ib2 = tr*sum(i2md.*cos(i2pd + thn));              	%   IB2 (= IC2)     
ib2rms = tr*sqrt(sum(i2m.^2)/2);
clear i2n i2p i2md i2pd

v2m = abs(v2n); v2md = v2m*ones(size(th));       	% bridge 2
v2p = angle(v2n); v2pd = v2p*ones(size(th)); 
vb2 = sum(v2md.*cos(v2pd + thn))/tr;    
vb2rms = sqrt(sum(v2m.^2)/2)/tr;
clear v2n v2m v2p v2md v2pd

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])/tics;

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])/tics;

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])/tics;

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])/tics;

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 transformer (second. referred)
Pr1 = sum(i1m.^2.*r1)/2; Pr2 = tr^2*sum(i2m.^2.*r2)/2; Pr = Pr1 + Pr2;
% Magnetising losses can't be modelled as a shunt R.
% Better approx is from Steinmetz: p = k B^2.8 f^1.46
fn = n.^1.46;                               % freq dependency
pl1 = 553e-5*sum(fn.*i1m.^2.8);           	% 15 W at m1,2 = 1
ptx = 346e-5*sum(fn.*il2m.^2.8);          	% 23 W at m1,2 = 1
Pmag = pl1 + ptx;                           % 38 W at 400 V
p1 = p0 + Pr2 + Pcon2 + Psw1 +  ptx/2;     	% adjust p1 and p2 for
p2 = p0 - Pr1 - Pcon1 - Psw2 - pl1 - ptx/2;	%  other losses 
loss = p1 - p2;                             % total losses

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

r = [m1; 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]; 	% 27 rows


% Test:     With tics = 12000, rws = 562, freq = 50k, temp = 75,
% r1,2 = fcn(freq), k2,3,4 ideal, STGWA19N60   
% [r] = clc_loss_2igbt(420,360,0.35,0.7,0.8,23,0.6,1.3) => 
%   0.7, 0.0, 64.4,	0.0, 670.9,	73.5, 17.1,	21.4, 76.8, 150.4,
%   0.0, 675.8,	0.0, 451.0,	112.7, 35.4, 40.4, 151.6, 264.2,
%   414.6, 136.6, 551.2, 5286.0, 5149.4, 18.4, 21.4

% Test: [sw,i1rms,ic2rms,p1,p2,Psw,Pcon,vsl1,vsm1,vsm2] = ...
%   clc_fcn_2igbt(420,360,0.35,0.8,0.8,13.7,0.6,1.3) =>









