// common emitter amplifier, inductor applied to remove dc, sinusoidal output voltage
// all voltages normalized by Vcc
// all currents normalized by Icq
// resistor normalized to r=(Icq/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;
pattern=sin(wt);

nr=400; // number of points over r, reduce if the simulation is too slow
rstop=4; // the final value for r

scf(1);

for k=1:nr
  
    r(k)=rstop*k/nr; // modulation index
    rr=r(k);
    if rr < 1
        jc=1+pattern;
        mce=1+rr-rr*jc;
    else
        mce=1-pattern;
        jc=1+1/rr-1/rr*mce;
    end
    mout=1-mce;
    jout=jc-1;
    pcc=1*jc; // 1* is for educational purposes, a symbol for Vcc*
    pd=mce.*jc;
    pout=mout.*jout;
    pl=1*mout; // 1* is for educational purposes, a symbol for Icq*
    Pcc(k)=mean(pcc);
    Pd(k)=mean(pd);
    Pout(k)=mean(pout);
    Pl(k)=mean(pl);
    efficiency(k)=Pout(k)/Pcc(k)*100;

    drawlater()    

    clf(1);

    subplot(4,2,1)
    plot2d(deg,mce,rect=[0 0 360 2],nax=[0 13 0 11])
    xtitle('','wt [deg]','vce/Vcc')

    subplot(4,2,3)
    plot2d(deg,jc,rect=[0 0 360 2],nax=[0 13 0 11])
    xtitle('','wt [deg]','ic/Icq')

    subplot(4,2,5)
    plot2d(deg,mout,rect=[0 -1 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','vout/Vcc')
    
    subplot(4,2,2)
    plot2d(deg,pcc,rect=[0 0 360 2],nax=[0 13 0 11])
    xtitle('','wt [deg]','pcc/(Vcc*Icq)')

    subplot(4,2,4)
    plot2d(deg,pd,rect=[0 0 360 2],nax=[0 13 0 11])
    xtitle('','wt [deg]','pd/(Vcc*Icq)')

    subplot(4,2,6)
    plot2d(deg,pout,rect=[0 0 360 2],nax=[0 13 0 11])
    xtitle('','wt [deg]','pout/(Vcc*Icq)')

    subplot(4,2,7)
    plot2d(mce,jc,style=0,rect=[0 0 2 2])
    xtitle('','vce/Vcc','ic/Icq')

    subplot(4,2,8)
    plot2d(deg,pl,rect=[0 -1 360 1],nax=[0 13 0 11])
    xtitle('','wt [deg]','pl/(Vcc*Icq)')

    drawnow()

end

scf(2);
clf(2);
plot2d(r,Pcc,rect=[0 0 4 1],style=color("black"))
plot2d(r,Pd,rect=[0 0 4 1],style=color("blue"))
plot2d(r,Pout,rect=[0 0 4 1],style=color("green"))
plot2d(r,Pl,rect=[0 0 4 1],style=color("red"))
xtitle('','R*Icq/Vcc','Pcc [black], Pd [blue], Pout [green], Pl [red], all /(Vcc*Icq)')

scf(3);
clf(3);
plot2d(r,efficiency,rect=[0 0 4 50])
xtitle('','R*Icq/Vcc','efficiency [%]')