% class B, complementary pair, reactive load
% all voltages normalized to Vcc
% all currents normalized to the output current amplitude
% Vbe assumed to be 10% of the Vcc

clear all
close all

%initialization
npoint=720; % reduce if the simulation is too slow
i=1:npoint;
wt=2*pi*(i-0.5)/npoint;
deg=wt*180/pi;
mout=sin(wt);
mce1=1-mout;
mec2=mout+1;

np=180; % number of frames, reduce if the simulation is too slow

for k=1:2*np
    
    % output current phase shift
    phi(k)=(-np+k-0.5)/np*pi/2;
    % computation
    jout=sin(wt-phi(k));
    q1=jout>0;
    q2=jout<0;
    min=mout+0.10*sign(jout);
    jc1=jout.*q1;
    jc2=-jout.*q2;
    jcc=jc1+jc2;
    pd1=jc1.*mce1;
    pd2=jc2.*mec2;
    pcc=jcc;
    pout=mout.*jout;
    Pcc(k)=mean(pcc);
    Pd1(k)=mean(pd1);
    Pd2(k)=mean(pd2);
    Pout(k)=mean(pout);
    efficiency(k)=Pout(k)/Pcc(k)*100;
    pdmax(k)=max(pd1);
    
    figure(1)
    
    subplot(4,3,1)
    plot(deg,jc1)
    xlabel('wt [deg]')
    ylabel('j_C_1')
    axis([0 360 -0.1 1.1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,2)
    plot(deg,mce1)
    xlabel('wt [deg]')
    ylabel('m_{CE1}')
    axis([0 360 -0.2 2.2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,3)
    plot(mce1,jc1,'.')
    xlabel('m_{CE1}')
    ylabel('j_{C1}')
    axis([-0.2 2.2 -0.1 1.1])
    
    subplot(4,3,4)
    plot(deg,jc2)
    xlabel('wt [deg]')
    ylabel('j_{C2}')
    axis([0 360 -0.1 1.1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,5)
    plot(deg,mec2)
    xlabel('wt [deg]')
    ylabel('m_{EC2}')
    axis([0 360 -0.2 2.2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,6)
    plot(mec2,jc2,'.')
    xlabel('m_{EC2}')
    ylabel('j_{C2}')
    axis([-0.2 2.2 -0.1 1.1])
    
    subplot(4,3,7)
    plot(deg,mout,'b',deg,min,'r')
    xlabel('wt [deg]')
    ylabel('m_{OUT} [blue], m_{IN} [red]')
    axis([0 360 -1.2 1.2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,8)
    plot(deg,mout,'b',deg,jout,'r')
    xlabel('wt [deg]')
    ylabel('m_{OUT} [blue], j_{OUT} [red]')
    axis([0 360 -1.2 1.2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,9)
    plot(deg,pout)
    xlabel('wt [deg]')
    ylabel('p_{OUT}')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    subplot(4,3,10)
    plot(deg,pd1)
    xlabel('wt [deg]')
    ylabel('p_{D1}')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,11)
    plot(deg,pd2)
    xlabel('wt [deg]')
    ylabel('p_{D2}')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,12)
    plot(deg,pcc)
    xlabel('wt [deg]')
    ylabel('p_{CC}')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

phi=phi*180/pi;

figure(2)
hold on
plot(phi,Pcc,'g')
plot(phi,Pd1,'b')
plot(phi,Pd1,'m')
plot(phi,Pout,'r')
plot(phi,Pd1+Pd2,'y')
xlabel('phi [deg]')
ylabel('P_{CC} [green], P_{D1} [blue], P_{D2} [magenta], P_{OUT} [red], P_{D1}+P_{D2} [yellow]')
axis([-90 90 0 0.7])
set(gca,'xtick',[-90 -60 -30 0 30 60 90])

figure(3)
plot(phi,efficiency)
xlabel('phi [deg]')
ylabel('efficiency [%]')
axis([-90 90 0 100])
set(gca,'xtick',[-90 -60 -30 0 30 60 90])

figure(4)
plot(phi,pdmax)
xlabel('phi [deg]')
ylabel('p_{D max}')
axis([-90 90 0 1.5])
set(gca,'xtick',[-90 -60 -30 0 30 60 90])