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

clear all
close all

%initialization
npoint=720;
i=1:npoint;
t=(i-0.5)/npoint;
pattern(1:npoint/2)=4*t(1:npoint/2)-1;
pattern(npoint/2+1:npoint)=3-4*t(npoint/2+1:npoint);
%for j=1:npoint
%    if j<=npoint/2
%        pattern(j)=4*t(j)-1;
%    else
%        pattern(j)=3-4*t(j);
%    end
%end

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

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

nmi=20; % number of simulations, various modulation indices

for k=1:nmi+1
    m(k)=(k-1)/nmi; % modulation index
    mce=MceQ+m(k)*Mm*pattern;
    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(t,mce)
    xlabel('t/T')
    ylabel('v_{CE}/V_{CC}')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,3)
    plot(t,jc)
    xlabel('t/T')
    ylabel('i_C/(V_{CC}/R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,5)
    plot(t,mr)
    xlabel('t/T')
    ylabel('v_R/V_{CC}')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,2)
    plot(t,pcc)
    xlabel('t/T')
    ylabel('p_{CC}/(V_{CC}^2/R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,4)
    plot(t,pd)
    xlabel('t/T')
    ylabel('p_D/(V_{CC}^2/R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,6)
    plot(t,pr)
    xlabel('t/T')
    ylabel('p_R/(V_{CC}^2/R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    subplot(4,2,7)
    plot(mce,jc,'.')
    xlabel('v_{CE}/V_{CC}')
    ylabel('i_C/(V_{CC}/R)')
    axis([0 1 0 1])
    subplot(4,2,8)
    plot(t,prac)
    xlabel('t/T')
    ylabel('p_{R ac}/(V_{CC}^2/R)')
    axis([0 1 0 1])
    set(gca,'xtick',[0 0.25 0.50 0.75 1])
    pause(1)
end

figure(2)
hold on
plot(m,Pcc,'r')
plot(m,Pd,'b')
plot(m,Pr,'g')
plot(m,Prac,'y')
xlabel('V_m/V_{m max}')
ylabel('P_{CC} [red], P_D [blue], P_R [red], P_{R ac} [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 [%]')