% clc_d_2igbt.m                                          	24-9-15

% Function calculates the network i/p and o/p currents, the volt.us 
%   products and the i/p and o/p powers in a symmetrical CLC DAB.
% For alalysis, between B1 and B2 is a CLC network of series L1 C1, 
%   shunt L3, a series L2 C2 to trfmr, B2. The ideal (tuned) 
%   reactances of L1 and L2, L3, C1 and C2 are, as a proportion 
%   of Xd: k1/2, 1, 1 + (k1/2). tr = trfmr turns ratio, n1/n2. 
% The VB1 lead, and the modulation angles are phi*pi, m1*pi and m2*pi.
% IGBTs in both bridges.

function [p1,p2,i1rms,ic2rms,vsm1,vsm2] = ...
    clc_d_2igbt(vdc1,vdc2,phi,m1,m2,xd,k1,tr,f)

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

% 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, (clc_netsolv_sym.m) are:
nsq = n.^2; den2 = k1 + 2 - nsq*(k1 + 4);
den = j*xd*(k1 + 2 - k1*nsq).*den2;
i1n = -2*n.*(2*v1n + k1*v1n + nsq.*(2*v2n - 2*v1n - k1*v1n))./den;
i2n = 2*n.*(2*v2n + k1*v2n + nsq.*(2*v1n - 2*v2n - k1*v2n))./den;
% vl3n = -2*nsq.*(v1n + v2n)./den2;
vm1n = v1n + 1i*i1n*(1 + k1/2)*xd./n;       % after C1
vm2n = v2n - 1i*i2n*(1 + k1/2)*xd./n;       % before C2 in analysis cct

% 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(1,tics);        	% bridge 1
v1 = sum(v1md.*cos(thn));                           % REFERENCE, 0 deg
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 i1n i1m i1p i1md i1pd 

vm1m = abs(vm1n); vm1md = vm1m*ones(size(th));    	% Vm1
vm1p = angle(vm1n); vm1pd = vm1p*ones(size(th)); 
vm1 = sum(vm1md.*cos(vm1pd + thn));    
clear vm1n vm1m vm1p vm1md vm1pd

vm2m = abs(vm2n); vm2md = vm2m*ones(size(th));    	% Vm2
vm2p = angle(vm2n); vm2pd = vm2p*ones(size(th)); 
vm2 = sum(vm2md.*cos(vm2pd + thn));    
clear vm2n vm2m vm2p vm2md vm2pd

% vl3m = abs(vl3n); vl3md = vl3m*ones(size(th));    	% vL3
% vl3p = angle(vl3n); vl3pd = vl3p*ones(size(th));     
% vl3 = sum(vl3md.*cos(vl3pd + thn));     
% clear vl3n vl3m vl3p vl3md vl3pd

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

% 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) = ic2(d:tics);                  % conducts st..end then
    iq6(1:d-tics/2-1) = ic2(1:d-tics/2-1);   	%   beginning of next
else
    iq6(d:d+tics/2-1) = ic2(d:d+tics/2-1);     	% conducts st..sp
end
iq8 = zeros(1,tics);
if c > (tics/2 + 1)                            	% for start > stop
    iq8(c:tics) = -ic2(c:tics);               	% conducts st..end then
    iq8(1:c-tics/2-1) = -ic2(1:c-tics/2-1);   	%   beginning of next
else
    iq8(c:c+tics/2-1) = -ic2(c:c+tics/2-1);   	% conducts st..sp
end

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*1e-6;                         % B1 switching power
Psw2 = Esw2*f*1e-6;                         % 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 volt.sec products
vsm1 = zeros(1,tics); vsm2 = vsm1;        	% duplicates
vsm1(1) = vm1(1); vsm2(1) = vm2(1);
for i = 2:tics
    vsm1(i) = vsm1(i-1) + vm1(i);           % trfmr prim with integ L1
    vsm2(i) = vsm2(i-1) + vm2(i);           % trfmr sec' with integ L2
end
fact = 1e6/f/tics;                          % scale for V.us units
vsm1 = (max(vsm1) - min(vsm1))*fact;
vsm2 = (max(vsm2) - min(vsm2))*fact/tr;     % trfmr secondary

p1 = p0 + Pcon2 + Psw1;     	% adjust p1 and p2 for
p2 = p0 - Pcon1 - Psw2;     	%  other 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

% Test: tics = 800, rows = 56, 
%   [p1,p2,i1rms,ic2rms,vsm1,vsm2] = 
%   clc_2igbt_sym (420,360,-0.65,0.9,0.8,27,0.6,1.1,5e4,50) =>
% p1 = 4362.60, p2 = 4175.6, loss = 186.4
% i1rms = 12.6885, ic2rms = 15.3067, vsm1 = 4341.9, vsm2 = 3181.7
% (pr = 26.2, pcon = 136.3, psw = 23.9)
% 
% Speed Test: 1000 loops => 4 seconds (16 s before conduction-loss
%       vectorization)









