// symmetrical class B amplifier with MOSFETSs
// Vdd=12V, m=1, n=1/2, R=2.5ohm, Vt=3V, B=0.5A/V^2

// 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);
m1=s>0;
m2=s<0;
splus=s.*m1;
sminus=-s.*m2;

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
    id1=Vm(k)/5*splus; // n/R*vout for vout>0
    id2=Vm(k)/5*sminus; // -n/R*vout for vout<0
    vin=(3+sqrt(4*id1)).*m1-(3+sqrt(4*id2)).*m2;
    idd=id1+id2;
    vds1=12-2*vout; // Vdd-vout/n
    vds2=12+2*vout; // Vdd+vout/n
    pd1=id1.*vds1;
    pd2=id2.*vds2;
    pdd=12*idd; // Vdd*idd
    pout=vout.*vout/2.5;
    Pdd(k)=mean(pdd);
    Pd1(k)=mean(pd1);
    Pd2(k)=mean(pd2);
    Pout(k)=mean(pout);
    efficiency(k)=Pout(k)/Pdd(k)*100;
    
    drawlater()
    
    clf(1);
    
    subplot(4,3,1)
    plot2d(deg,id1,rect=[0 -0.25 360 1],nax=[0 7 0 6])
    xtitle('','wt [deg]','id1 [A]')

    subplot(4,3,2)
    plot2d(deg,vds1,rect=[0 0 360 25],nax=[0 7 0 6])
    xtitle('','wt [deg]','vds1 [V]')
    
    subplot(4,3,3)
    plot2d(vds1,id1,rect=[0 -0.25 25 1],nax=[0 6 0 6],style=0)
    xtitle('','vds1 [V]','id1 [A]')
    
    subplot(4,3,4)
    plot2d(deg,id2,rect=[0 -0.25 360 1],nax=[0 7 0 6])
    xtitle('','wt [deg]','id2 [A]')
    
    subplot(4,3,5)
    plot2d(deg,vds2,rect=[0 0 360 25],nax=[0 7 0 6])
    xtitle('','wt [deg]','vds2 [V]')
    
    subplot(4,3,6)
    plot2d(vds2,id2,rect=[0 -0.25 25 1],nax=[0 6 0 6],style=0)
    xtitle('','vds2 [V]','id2 [A]')
        
    subplot(4,3,7)
    plot2d(deg,vin,rect=[0 -5 360 5],nax=[0 7 0 6])
    xtitle('','wt [deg]','vin [V]')

    subplot(4,3,8)
    plot2d(deg,vout,rect=[0 -5 360 5],nax=[0 7 0 6])
    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,pdd,rect=[0 -3 360 12],nax=[0 7 0 6])
    xtitle('','wt [deg]','pdd [W]')
        
    drawnow()
    
end

scf(2);
clf(2);
plot2d(Vm,Pdd,rect=[0 0 5 8],style=color("green"))
plot2d(Vm,Pd1,rect=[0 0 5 8],style=color("blue"))
plot2d(Vm,Pd2,rect=[0 0 5 8],style=color("magenta"))
plot2d(Vm,Pout,rect=[0 0 5 8],style=color("red"))
plot2d(Vm,Pd1+Pd2,rect=[0 0 5 8],style=color("yellow"))
xtitle('','Vm [V]','Pdd [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]')