% clc_dab_2igbt.m                                               6-6-15

% Calculate and display (if plott = 1) voltages and currents in CLC DAB.
% 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.  tr = trfmr turns ratio, n1/n2.
% The loss calculations assume that each bridge uses igbts.

Plott = 1;                                  % 1 => yy plots are done
tics = 800;                                 % number of points, columns
rows = 566;     n = (1:2:2*rows - 1)';     	% rows = odd freqencies
f = 5e4;                                    % for plot axis t labels
temp = 75;          	% Enter IGBT (STGWA19NC60) junction temperature

% Enter the parameter values :
% The reactances of the network legs equal xb at the fs fundamental.
%	XL1 = k1*n*xb, XC1 = (1 + k1)*xb/n, XL2 = xb*n, XC2_ = xb/n.
% k1 = 1;
k1 = 1.0; k2 = 1+k1; k3 = 1; k4 = 1;
m1 = 0.7;                                     % 0 < modulation index < 1
% m1 = 1.0;                                     % 0 < modulation index < 1
m2 = m1;
% m2 = 0.81;                                  % modulation index, 0..1
phi = -0.5;                                  % v1 lead index, -ve => fwd
vdc1 = 400;                               	% DC supply to Bridge 1
vdc2 = 400;                              	% DC supply to Bridge 2
tr = 1;                              	% tr = trfmr turns ratio n1/n2
xd = 31.83;                                  
% vdc1 = 765, vdc2 = 432, phi = -0.6, m1 = 0.7, m2 = 0.8,     % TEST
% k1 = 0.6, k2 = 1.7, k3 = 1.1, k4 = 1.2, tr = 2, xd = 55,    % TEST

% Actual parameters (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; r1 = 0.133; r2 = 0.126;
% 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;

% den = xd*(n.*(k2*r2_ + k4*r1) - n.^3.*(k1*r2_ + k3*r1 + k3*r2_)) ...
%     + 1i*(n.^2.*r1.*r2_ - xd^2*(k2*k4 - n.^2*(k1*k4 + k2*k3 + k3*k4) ...
%     + k1*k3*n.^4));
% i1n = -(xd*((k3*(v1n - v2n)).*n.^3 - k4*v1n.*n) - 1i*n.^2.*r2_.*v1n)./den;
% i2n = (xd*((k1*v2n - k3*(v1n - v2n)).*n.^3 - k2*v2n.*n) ...
%     - n.^2.*r1.*v2n*1i)./den; 

vc2n = -1i*i2n*xd./n/tr;                      	% actual C2, trfmr RHS
vl1n = 1i*i1n*k1*xd.*n;                        	% vL1 (I1 = IL1 = IC1)
il2n = i1n - i2n;                              	% current in L2 shunt
vl2n = 1i*il2n*xd.*n;                         	% vL2 
vc1n = -1i*i1n*(1 + k1)*xd./n;                  % vC1

% 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));    
v1rms = sqrt(sum(v1m.^2)/2); v1max = max(v1); v1min = min(v1);
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); i1max = max(i1); i1min = min(i1);
clear i1n i1p i1md i1pd 

vc1m = abs(vc1n); vc1md = vc1m*ones(size(th));          % vC1  
vc1p = angle(vc1n); vc1pd = vc1p*ones(size(th));   
vc1 = sum(vc1md.*cos(vc1pd + thn));       
vc1rms = sqrt(sum(vc1m.^2)/2); vc1max = max(vc1); vc1min = min(vc1);
clear vc1n vc1m vc1p vc1md vc1pd

vl1m = abs(vl1n); vl1md = vl1m*ones(size(th));          % vL1
vl1p = angle(vl1n); vl1pd = vl1p*ones(size(th));     
vl1 = sum(vl1md.*cos(vl1pd + thn));     
vl1rms = sqrt(sum(vl1m.^2)/2); vl1max = max(vl1); vl1min = min(vl1);
clear vl1n vl1m vl1p vl1md vl1pd

vl2m = abs(vl2n); vl2md = vl2m*ones(size(th));      	% vL2
vl2p = angle(vl2n); vl2pd = vl2p*ones(size(th));     
vl2 = sum(vl2md.*cos(vl2pd + thn));     
vl2rms = sqrt(sum(vl2m.^2)/2); vl2max = max(vl2); vl2min = min(vl2);
clear vl2n vl2m vl2p vl2md vl2pd

il2m = abs(il2n); il2md = il2m*ones(size(th));          % iL2  
il2p = angle(il2n); il2pd = il2p*ones(size(th));   
il2 = sum(il2md.*cos(il2pd + thn));       
il2rms = sqrt(sum(il2m.^2)/2); il2max = max(il2); il2min = min(il2);
clear il2n il2p il2md il2pd

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); ic2max = max(ib2); ic2min = min(ib2);
clear i2n i2p i2md i2pd

vc2m = abs(vc2n); vc2md = vc2m*ones(size(th));     	% vC2, trfmr RHS  
vc2p = angle(vc2n); vc2pd = vc2p*ones(size(th));   
vc2 = sum(vc2md.*cos(vc2pd + thn));       
vc2rms = sqrt(sum(vc2m.^2)/2); vc2max = max(vc2); vc2min = min(vc2);
clear vc2n vc2m vc2p vc2md vc2pd

v2r = v2n/tr;                                       % line voltage = VB2
v2rm = abs(v2r); v2rmd = v2rm*ones(size(th));     	% trfmr RHS =   	
v2rp = angle(v2r); v2rpd = v2rp*ones(size(th));     
vb2 = sum(v2rmd.*cos(v2rpd + thn));             	% line voltage
vb2rms = sqrt(sum(v2rm.^2)/2); vb2max = max(vb2); vb2min = min(vb2);                    	
clear v2n v2rn v2rm v2rp v2rmd v2rpd  

vm1 = v1 - vc1;

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

% figure (2);
% plot(th,p)
% ylabel('Instantaneous Power (W)','fontsize',12)
% xlabel('Power Cycle Angle (rad)','fontsize',12)

% disp('The average power input from VDC1, output to VDC2 are:');
% p = (p1 + p2)/2

disp('The bridge values V1, VB2, I1 (= IC1 = IL1), IB2 (= IC2) are :')
ac(1,:) = [v1rms v1max v1min];
ac(2,:) = [vb2rms vb2max vb2min];
ac(3,:) = [i1rms i1max i1min];
ac(4,:) = [ib2rms ic2max ic2min];
ac(1:4,:)

disp('The ic1, vc1, vl1, il2, vl2, ic2, vc2 rms, max and min values are')
ac(5,:) = [i1rms i1max i1min];
ac(6,:) = [vc1rms vc1max vc1min];
ac(7,:) = [vl1rms vl1max vl1min];
ac(8,:) = [il2rms il2max il2min];
ac(9,:) = [vl2rms vl2max vl2min];
ac(10,:) = [ib2rms ic2max ic2min];
ac(11,:) = [vc2rms vc2max vc2min];
ac(5:11,:)

if Plott == 1                 	% requested plot o/p
ext = uint16(tics/2);                               % extention, 0.5 cycle
t = 0:1e6/f/tics:(tics*1.5-1)*1e6/f/tics;     % time in usecs
v1d = [v1 v1(1:ext)];
i1d = [i1 i1(1:ext)];
v2d = [vb2 vb2(1:ext)];
i2d = [ib2 ib2(1:ext)];
pd = [p p(1:ext)];

figure(1);
subplot(2,1,1)                      	% V1, I1
[AX,H1,H2] = plotyy(t,v1d,t,i1d);
set(get(AX(1),'Ylabel'),'String','v_{1}(V)','FontSize',14)
set(get(AX(2),'Ylabel'),'String','i_{1}(A)','FontSize',14)
set(AX(1),'YTickMode','auto');
set(H1,'LineStyle','-','LineWidth',2.0); 
set(H2,'LineStyle','-','LineWidth',2.0);
set(AX(1),'XTickLabel',[]);             % clear the x axes
set(AX(2),'XTickLabel',[]);
set(AX(1),'Ylim',[-500 500]);
yends = get(AX(1),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(1),'YTick',yends(1):yinc:yends(2))
set(AX(1),'YTickLabel',get(AX(1),'YTickLabel'),'FontSize',14)
set(AX(2),'Ylim',[-20 20]);
yends = get(AX(2),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(2),'YTick',yends(1):yinc:yends(2))
set(AX(2),'YTickLabel',get(AX(2),'YTickLabel'),'FontSize',14)
hleg = legend('v_{1}','i_{1}');
set (hleg,'FontSize',12,'location','southeast')
grid on;
subplot(2,1,2)                              % H handles are for line
[AX,H1,H2] = plotyy(t,v2d,t,i2d);           % AX handles are for axis
set(get(AX(1),'Ylabel'),'String','v_{2}(V)','FontSize',14)
set(get(AX(2),'Ylabel'),'String','i_{2}(A)','FontSize',14)
set(H1,'LineStyle','-','LineWidth',2)  	
set(H2,'LineStyle','-','LineWidth',2)  
xlabel('Time (\mus)','FontSize',14)             % "\mu" => greek "micro"
set(AX(1),'Ylim',[-500 500]);
yends = get(AX(1),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(1),'YTick',yends(1):yinc:yends(2))
set(AX(1),'YTickLabel',get(AX(1),'YTickLabel'),'FontSize',14)
set(AX(2),'Ylim',[-20 20]);
yends = get(AX(2),'YLim');
yinc = (yends(2)-yends(1))/4;
set(AX(2),'YTick',yends(1):yinc:yends(2))
set(AX(2),'YTickLabel',get(AX(2),'YTickLabel'),'FontSize',14)
hleg = legend('v_{2}','i_{2}');
set (hleg,'FontSize',14,'location','southeast')
grid on;
clear ext t i1d i2d pd mid1 mid2 v1d v2d Plott hleg AX H1 H2 yends yinc
end

% 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
iin = iq2 + iq4;                                % current from VDC1 to B1
sw1rms = sqrt(sum((iq2').^2)/tics);             % leg 1 rms current
sw2rms = sqrt(sum((iq4').^2)/tics);             % leg 2 rms current
% figure(3);  plot(th, iq2)                  	% diagnostic
% figure(4);  plot(th, iq4)                  	% diagnostic
% figure(9);  plot(th, iin)                   	% diagnostic

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
iout = -iq6 - iq8;                              % current from B2 to VDC2 
sw3rms = sqrt(sum((iq6').^2)/tics);             % leg 3 rms current
sw4rms = sqrt(sum((iq8').^2)/tics);             % leg 4 rms current
% figure(7);  plot(th, iq6)                  	% diagnostic
% figure(8);  plot(th, iq8)                  	% diagnostic
% figure(9);  plot(th, iout)                  	% diagnostic

iin_avg = sum(iin)/tics;
iin_ms = sum(iin.^2)/tics;                      % mean square
iin_rms = sqrt(iin_ms);
iin_rms_ac = sqrt(iin_ms - iin_avg^2);          % rms ripple current
iin_max = max(iin); iin_min = min(iin);
p0 = vdc1*iin_avg;                              % average pwr throughput 

iout_avg = sum(iout)/tics;
iout_ms = sum(iout.^2)/tics;                  	% mean square
iout_rms = sqrt(iout_ms);
iout_rms_ac = sqrt(iout_ms - iout_avg^2);     	% rms ripple current
iout_max = max(iout); iout_min = min(iout);

% 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 in W
Psw2 = Esw2*f*1e-6;                    % B2 switching power in W
Psw = Psw1 + Psw2                       	% total 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])/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);                  	%  
Pcon = Pcon1 + Pcon2                      % total conduction losses

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

disp('Bridge legs 1..4 (down) have switch-on, maximum, minimum and rms')
disp('device currents of :')
sw(1,:) = [iq2(a) max(iq2) min(iq2) sw1rms];
sw(2,:) = [iq4(b) max(iq4) min(iq4) sw2rms];  
sw(3,:) = [iq6(d) max(iq6) min(iq6) sw3rms];
sw(4,:) = [iq8(c) max(iq8) min(iq8) sw4rms];
sw

disp('The avg, rms, ac_rms, max and min values (across) for the bridge')
disp(' input and output dc currents are :')
dc(1,:) = [iin_avg iin_rms iin_rms_ac iin_max iin_min];
dc(2,:) = [iout_avg iout_rms iout_rms_ac iout_max iout_min];
dc

% Calculate the volt.sec products
vsl1 = zeros(1,tics); vsm1 = vsl1; vsm2 = vsl1;             % duplicates
vsl1(1) = vl1(1); vsm1(1) = vm1(1); vsm2(1) = vl2(1);
for i = 2:tics
  	vsl1(i) = vsl1(i-1) + vl1(i);           % L1
    vsm1(i) = vsm1(i-1) + vm1(i);           % trfmr prim with integ L1
    vsm2(i) = vsm2(i-1) + vl2(i);           % trfmr primary (L2 = mut)
end
fact = 1e6/f/tics;                     	% scale with period in us
vsl1 = (max(vsl1) - min(vsl1))*fact;
vsm1 = (max(vsm1) - min(vsm1))*fact;
vsm2 = (max(vsm2) - min(vsm2))*fact;
apl1 = i1rms*vsl1; apm1 = i1rms*vsm1; 
apm2 = i1rms*vsm2; apm3 = ib2rms*vsm2/tr;
% indice 3 refers to the secondary voltage - derived from primary/tr
apdisc = apl1 + apm2 + apm3;
apinteg = apm1 + apm3;
disp('The V.us products of L1, Primary+L1, and Primary are :')
[vsl1 vsm1 vsm2]
disp('The V.us.Irms Products for L1, Primary+L1, Primary and Secondary are :')
[apl1 apm1 apm2 apm3]
disp('The V.us.Irms Discrete and Integ L1 Products are :')
[apdisc apinteg]
p2
eff=p2/p1

% clear n pwr p th a b c d e f g h iin_ms iout_ms

% Test:  With tics = 12000, rws = 562, freq = 50k, temp = 75,
% vdc1 = 765; vdc2 = 432; thi = -0.6; m1 = 0.7; m2 = 0.8;       % TEST
% k1 = 0.6; k2 = 1.7; k3 = 1.1; k4 = 1.2; tr = 2; xd = 55;      % TEST
% r1 = 0.2, r2 = 0.1, IGBTD_STGWA19N60  :

% p0 = 7216.1,  p1 = 7483.6, p2 = 7074.5
%  V1, VB2, I1, IB2 are : 639.96, 386.35, 12.4850, 22.5659 
% Psw = 110.1, Pcon = 195.6, Pr = 82.1, Pmag = 21.3, loss = 409.1
% sw =  13.1003   15.4117  -13.6547    8.8282
%       -8.6811   15.4117   -8.6811    8.8282
%       -13.9791   13.9698  -43.5027   15.9565
%       6.3967    6.3967  -43.5027   15.9565
% The V.us products = 3237.8    9663.4    6488.4
% Discrete and Integ AP = 194.6k, 193.9k




