// common emitter amplifier, sinusoidal output voltage
// all voltages normalized to Vcc
// all currents normalized to Vcc/R

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

MceQ=1/2; // you can adjust this, move the quiescent operating point from the optimal position
JcQ=1-MceQ;

np=200; // number of frames, reduce if the simulation is too slow

Mm=min([MceQ, 1-MceQ]);

scf(1);

for k=1:np
  
    m(k)=k/np; // modulation index
    mce=MceQ+m(k)*Mm*sin(wt);
    jc=1-mce;
    mr=1-mce;
    pcc=1*jc; // 1* is for educational purposes, a symbol for Vcc*
    pd=mce.*jc;
    pr=mr.*jc;
    prac=(jc-JcQ).*(mr-(1-MceQ));
    Pcc(k)=mean(pcc);
    Pd(k)=mean(pd);
    Pr(k)=mean(pr);
    Prac(k)=mean(prac);
    efficiency(k)=Prac(k)/Pcc(k)*100;
    
    drawlater()
    
    clf(1);
    
    subplot(421)
    plot2d(deg,mce,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','vce/Vcc')
      
    subplot(423)
    plot2d(deg,jc,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','ic/(Vcc/R)')
    
    subplot(425)
    plot2d(deg,mr,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','vr/Vcc')
    
    subplot(422)
    plot2d(deg,pcc,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','pcc/(Vcc^2/R)')
    
    subplot(424)
    plot2d(deg,pd,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','pd/(Vcc^2/R)')
    
    subplot(426)
    plot2d(deg,pr,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','pr/(Vcc^2/R)')
    
    subplot(427)
    plot2d(mce,jc,style=0,rect=[0 0 1 1])
    xtitle('','vce/Vcc','ic/(Vcc/R)')
    
    subplot(428)
    plot2d(deg,prac,rect=[0 0 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','prac/(Vcc^2/R)')

    drawnow()
    
end

scf(2);
clf(2);
plot2d(m,Pcc,rect=[0 0 1 1.1*JcQ],style=color("black"))
plot2d(m,Pr,rect=[0 0 1 1.1*JcQ],style=color("blue"))
plot2d(m,Prac,rect=[0 0 1 1.1*JcQ],style=color("green"))
plot2d(m,Pd,rect=[0 0 1 1.1*JcQ],style=color("red"))
xtitle('','Vm/Vmmax','Pcc [black], Pr [blue], Pd [green], Prac [red], all /(Vcc^2/R)')

scf(3);
clf(3);
plot2d(m,efficiency)
xtitle('','Vm/Vmmax','efficiency [%]')