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

clear all
close all

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

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

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;
    figure(1)
    subplot(4,2,1)
    plot(deg,mce)
    xlabel('wt [deg]')
    ylabel('v_{CE} / V_{CC}')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,3)
    plot(deg,jc)
    xlabel('wt [deg]')
    ylabel('i_C / (V_{CC} / R)')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,5)
    plot(deg,mr)
    xlabel('wt [deg]')
    ylabel('v_R / V_{CC}')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,2)
    plot(deg,pcc)
    xlabel('wt [deg]')
    ylabel('p_{CC} / (V_{CC}^2 / R)')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,4)
    plot(deg,pd)
    xlabel('wt [deg]')
    ylabel('p_D / (V_{CC}^2 / R)')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,6)
    plot(deg,pr)
    xlabel('wt [deg]')
    ylabel('p_R / (V_{CC}^2 / R)')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    subplot(4,2,7)
    plot(mce,jc,'.')
    xlabel('v_{CE} / V_{CC}')
    ylabel('i_C / (V_{CC} / R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.5 0.75 1])
    subplot(4,2,8)
    plot(deg,prac)
    xlabel('wt [deg]')
    ylabel('p_{Rac} / (V_{CC}^2/R)')
    axis([0 360 0 1])
    set(gca,'xtick',[0 30 60 90 120 150 180 210 240 270 300 330 360])
    pause(1)
end

figure(2)
hold on
plot(m,Pcc,'g')
plot(m,Pd,'b')
plot(m,Pr,'r')
plot(m,Prac,'y')
xlabel('V_m / V_{m max}')
ylabel('P_{CC} [green], P_D [blue], P_R [red], P_{Rac} [yellow], all / (V_{CC}^2/R)')
axis([0 1 0 1.1*JcQ])

figure(3)
plot(m,efficiency)
xlabel('V_m / V_{m max}')
ylabel('efficiency [%]')