// common emitter amplifier, inductor applied to remove dc, influence of the
// inductor current ac ripple
// all voltages normalized by Vcc
// all currents normalized by Icq
// inductor normalized to l=(w L Icq/Vcc)
// optimal load resistor, R=Vcc/Icq

//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;
ps=sin(wt);
pc=cos(wt);

np=100; // number of points over l, reduce if the simulation is too slow
ndec=4; // number of decades to be covered
l0=10^(-ndec/2);
lstep=10^(ndec/np);

scf(1);

for k=1:np
  
    l0=l0*lstep;
    l(k)=l0; // actual normalized inductance
    p1=l0/sqrt(1+l0^2);
    Vm(k)=p1; // amplitude of the output voltage
    p2=1/sqrt(1+l0^2);
    mce=1-p1*ps;
    jl=1-p2*pc;
    jout=p1*ps;
    jc=jl+jout;
    mout=1-mce;
    pcc=1*jc; // 1* is for educational purposes, a symbol for Vcc*
    pd=mce.*jc;
    pout=mout.*jout;
    pl=jl.*mout;
    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(l,Pcc,rect=[min(l) -0.1 max(l) 1.1],style=color("black"))
plot2d(l,Pd,rect=[min(l) -0.1 max(l) 1.1],style=color("blue"))
plot2d(l,Pout,rect=[min(l) -0.1 max(l) 1.1],style=color("green"))
plot2d(l,Pl,rect=[min(l) -0.1 max(l) 1.1],style=color("red"))
xtitle('','w*L*Icq/Vcc','Pcc [black], Pd [blue], Pout [green], Pl [red], all /(Vcc*Icq)')

scf(3);
clf(3);
plot2d(l,efficiency,rect=[min(l) -25 max(l) 75])
xtitle('','w*L*Icq/Vcc','efficiency [%]')

scf(4);
clf(4);
plot2d(l,Vm,logflag='ll')
xtitle('','w*L*Icq/Vcc','Vm/Vcc')