% L_d_2igbt.m                                             	1-8-15

% Calculates the i/p and o/p currents, the integrated volt.us products
%    and the i/p and o/p powers in Split L DAB.
% Between B1 and B2 is L1 R1 series, Lm shunt, transformer, L2 R2 series.
% Network solution is for L1 = L/2, Lm = km*L, L2 = L/2 tee conecting
%  to ideal transformer. Actual circuit has L2 = L/(2tr^2) on trfmr sec.
% The VB1 lead, and the modulation angles are phi*pi, m1*pi and m2*pi.
% 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).

function [p1,p2,i1rms,ib2rms,vsm1,vsm2] = ...
    L_d_2igbt(vdc1,vdc2,phi,m1,m2,xd,km,tr,f)

tics = 800;                          	% number of points, columns
n = (1:2:121)';                       	% rows = odd freqencies
temp = 75;          	% Enter IGBT (STGWA19NC60) junction temperature

% Calculate the voltage phasors of each harmonic
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, 62 x tics

den = n*xd*(1 + 4*km); 
i1n = -((v1n + 2*km*v1n - 2*km*v2n)*2*1i)./den;
ib2n = tr*((v2n - 2*km*v1n + 2*km*v2n)*2*1i)./den;   	% out trfmr sec
% vmn = (2*km*(v1n + v2n))/(4*km + 1);

% 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 v1m v1p v1md v1pd

i1m = abs(i1n); i1md = i1m*ones(size(th));           	% I1 (= IB1, 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 

% vmm = abs(vmn); vmmd = vmm*ones(size(th));            	% vm
% vmp = angle(vmn); vmpd = vmp*ones(size(th)); 
% vm = sum(vmmd.*cos(vmpd + thn));    
% clear vmm vmp vmmd vmpd

% vb2m = abs(v2n)/tr; vb2md = vb2m*ones(size(th));     	% bridge 2   	
% vb2p = angle(v2n); vb2pd = vb2p*ones(size(th));     
% vb2 = sum(vb2md.*cos(vb2pd + thn));     
% vb2rms = sqrt(sum(vb2m.^2)/2);                    	
% clear v2n vB2m vB2p vB2md vB2pd  

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); 
% clear ib2n ib2m ib2p ib2md ib2pd

% 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
Psw1 = 2*(Eon2 + Eoff2 + Eon4 + Eoff4)*f*1e-6;      % B1 switching power
Psw2 = 2*(Eon6 + Eoff6 + Eon8 + Eoff8)*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 magnetising volt.sec products: 
vsm1 = vdc1*m1*5e5/f;                    	% V.us product 
vsm2 = vdc2*m2*5e5/f;                       % V.us product 
% Pmag = (vs1^2.8 + vs2^2.8)*pm*f^1.46/55e13;	% trfmr mag loss    	

% Pr1 = i1rms^2*r1;                           % resistive losses 
% Pr2 = ib2rms^2*r2;                        	% resistive losses
% p1 = p0 + Pcon1 + Psw1 + Pr2;             	% adjust p1 and p2 for swg
% p2 = p0 - Pcon2 - Psw2 - Pr1;              	%  and cond losses
p1 = p0 + Pcon1 + Psw1;                     % adjust p1 and p2 for swg
p2 = p0 - Pcon2 - 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



