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

// initialization
clear
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=200; // number of frames, reduce if the simulation is too slow

scf(1);

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;
    
    drawlater();
    
    clf(1);

    subplot(4,3,1)
    plot2d(deg,ic1,rect=[0 -0.25 360 1.25],nax=[0 7 0 7])
    xtitle('','wt [deg]','ic1 [A]')

    subplot(4,3,2)
    plot2d(deg,vce1,rect=[0 0 360 20],nax=[0 7 0 5])
    xtitle('','wt [deg]','vce1 [V]')
    
    subplot(4,3,3)
    plot2d(vce1,ic1,rect=[0 -0.25 20 1.25],nax=[0 5 0 7],style=0)
    xtitle('','vce1 [V]','ic1 [A]')
    
    subplot(4,3,4)
    plot2d(deg,ic2,rect=[0 -0.25 360 1.25],nax=[0 7 0 7])
    xtitle('','wt [deg]','ic2 [A]')
    
    subplot(4,3,5)
    plot2d(deg,vce2,rect=[0 0 360 20],nax=[0 7 0 5])
    xtitle('','wt [deg]','vce2 [V]')
    
    subplot(4,3,6)
    plot2d(vce2,ic2,rect=[0 -0.25 20 1.25],nax=[0 5 0 7],style=0)
    xtitle('','vce2 [V]','ic2 [A]')
        
    subplot(4,3,7)
    plot2d(deg,vin,rect=[0 -1 360 1],nax=[0 7 0 5])
    xtitle('','wt [deg]','vin [V]')

    subplot(4,3,8)
    plot2d(deg,vout,rect=[0 -5 360 5],nax=[0 7 0 5])
    xtitle('','wt [deg]','vout [V]')
 
    subplot(4,3,9)
    plot2d(deg,pout,rect=[0 -3 360 12],nax=[0 7 0 6])
    xtitle('','wt [deg]','pout [W]')
        
    subplot(4,3,10)
    plot2d(deg,pd1,rect=[0 -3 360 12],nax=[0 7 0 6])
    xtitle('','wt [deg]','pd1 [W]')

    subplot(4,3,11)
    plot2d(deg,pd2,rect=[0 -3 360 12],nax=[0 7 0 6])
    xtitle('','wt [deg]','pd2 [W]')

    subplot(4,3,12)
    plot2d(deg,pcc,rect=[0 -3 360 12],nax=[0 7 0 6])
    xtitle('','wt [deg]','pcc [W]')
        
    drawnow()
       
end

scf(2);
clf(2);
plot2d(Vm,Pcc,rect=[0 0 5 7],style=color("green"))
plot2d(Vm,Pd1,rect=[0 0 5 7],style=color("blue"))
plot2d(Vm,Pd2,rect=[0 0 5 7],style=color("magenta"))
plot2d(Vm,Pout,rect=[0 0 5 7],style=color("red"))
plot2d(Vm,Pd1+Pd2,rect=[0 0 5 7],style=color("yellow"))
xtitle('','Vm [V]','Pcc [green], Pd1 [blue], Pd2 [magenta], Pout [red], Pd1+Pd2 [yellow], all in [W]')

scf(3);
clf(3);
plot2d(Vm,efficiency)
xtitle('','Vm [V]','efficiency [%]')

scf(4);
clf(4);
plot2d(vin,vout)
xtitle('','vin [V]','vout [V]')