% symmetrical class B amplifier with BJTs
% Vcc=10V, m=1, n=1/2, R=2.5ohm

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;
s=sin(wt);
q1=s>0;
q2=s<0;
splus=s.*q1;
sminus=-s.*q2;

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

for k=1:np
    
    % computation of voltages and currents
    Vm(k)=k/np*5; % amplitude of the output voltage
    vout=Vm(k)*s; % amplitude * pattern
    ic1=Vm(k)/5*splus; % n/R*vout for vout>0
    ic2=Vm(k)/5*sminus; % -n/R*vout for vout<0
    vin=0.025*((log(1+ic1*exp(30))).*q1-(log(1+ic2*exp(30))).*q2);
    icc=ic1+ic2;
    vce1=10-2*vout; % Vcc-vout/n
    vce2=10+2*vout; % Vcc+vout/n
    pd1=ic1.*vce1;
    pd2=ic2.*vce2;
    pcc=10*icc; % Vcc*icc
    pout=vout.*vout/2.5;
    Pcc(k)=mean(pcc);
    Pd1(k)=mean(pd1);
    Pd2(k)=mean(pd2);
    Pout(k)=mean(pout);
    efficiency(k)=Pout(k)/Pcc(k)*100;
    
    figure(1)
    
    subplot(4,3,1)
    plot(deg,ic1)
    xlabel('w_0t [deg]')
    ylabel('i_{C1} [A]')
    axis([0 360 0 1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,2)
    plot(deg,vce1)
    xlabel('w_0t [deg]')
    ylabel('v_{CE1} [V]')
    axis([0 360 0 20])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,3)
    plot(vce1,ic1,'.')
    xlabel('v_{CE1} [V]')
    ylabel('i_{C1} [A]')
    axis([-1 21 -0.1 1.1])
    set(gca,'xtick',[0 5 10 15 20])
        
    subplot(4,3,4)
    plot(deg,ic2)
    xlabel('w_0t [deg]')
    ylabel('i_{C2} [A]')
    axis([0 360 0 1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,5)
    plot(deg,vce2)
    xlabel('w_0t [deg]')
    ylabel('v_{CE2} [V]')
    axis([0 360 0 20])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,6)
    plot(vce2,ic2,'.')
    xlabel('v_{CE2} [V]')
    ylabel('i_{C2} [A]')
    axis([-1 21 -0.1 1.1])
    set(gca,'xtick',[0 5 10 15 20])
    
    subplot(4,3,7)
    plot(deg,vin)
    xlabel('w_0t [deg]')
    ylabel('v_{IN} [V]')
    axis([0 360 -1 1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,8)
    plot(deg,vout)
    xlabel('w_0t [deg]')
    ylabel('v_{OUT} [V]')
    axis([0 360 -5 5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,9)
    plot(deg,pout)
    xlabel('w_0t [deg]')
    ylabel('p_{OUT} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    subplot(4,3,10)
    plot(deg,pd1)
    xlabel('w_0t [deg]')
    ylabel('p_{D1} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,11)
    plot(deg,pd2)
    xlabel('w_0t [deg]')
    ylabel('p_{D2} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,12)
    plot(deg,pcc)
    xlabel('w_0t [deg]')
    ylabel('p_{CC} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

figure(2)
hold on
plot(Vm,Pcc,'g')
plot(Vm,Pd1,'b')
plot(Vm,Pd1,'m')
plot(Vm,Pout,'r')
plot(Vm,Pd1+Pd2,'y')
xlabel('V_m [V]')
ylabel('P_{CC} [green], P_{D1} [blue], P_{D2} [magenta], P_{OUT} [red], P_{D1}+P_{D2} [yellow], all in [W]')
axis([0 5 0 8])

figure(3)
plot(Vm,efficiency)
xlabel('V_m [V]')
ylabel('efficiency [%]')

figure(4)
plot(vin,vout)
xlabel('v_{IN}')
ylabel('v_{OUT}')