% symmetrical class B amplifier with MOSFETSs
% Vdd=12V, m=1, n=1/2, R=2.5ohm, Vt=3V, B=0.5A/V^2

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;
s=sin(wt);
m1=s>0;
m2=s<0;
splus=s.*m1;
sminus=-s.*m2;

np=20; % number of frames, reduce if the simulation is too slow

for k=1:np
    
    % computation of voltages and currents
    Vm(k)=k/np*5; % amplitude of the output voltage
    vout=Vm(k)*s; % amplitude * pattern
    id1=Vm(k)/5*splus; % n/R*vout for vout>0
    id2=Vm(k)/5*sminus; % -n/R*vout for vout<0
    vin=(3+sqrt(4*id1)).*m1-(3+sqrt(4*id2)).*m2;
    idd=id1+id2;
    vds1=12-2*vout; % Vdd-vout/n
    vds2=12+2*vout; % Vdd+vout/n
    pd1=id1.*vds1;
    pd2=id2.*vds2;
    pdd=12*idd; % Vdd*idd
    pout=vout.*vout/2.5;
    Pdd(k)=mean(pdd);
    Pd1(k)=mean(pd1);
    Pd2(k)=mean(pd2);
    Pout(k)=mean(pout);
    efficiency(k)=Pout(k)/Pdd(k)*100;
    
    figure(1)
    
    subplot(4,3,1)
    plot(deg,id1)
    xlabel('w_0t [deg]')
    ylabel('i_{D1} [A]')
    axis([0 360 0 1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,2)
    plot(deg,vds1)
    xlabel('w_0t [deg]')
    ylabel('v_{DS1} [V]')
    axis([0 360 0 25])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,3)
    plot(vds1,id1,'.')
    xlabel('v_{DS1} [V]')
    ylabel('i_{D1} [A]')
    axis([-2 27 -0.1 1.1])
    set(gca,'xtick',[0 5 10 15 20 25])
   
    subplot(4,3,4)
    plot(deg,id2)
    xlabel('w_0t [deg]')
    ylabel('i_{D2} [A]')
    axis([0 360 0 1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,5)
    plot(deg,vds2)
    xlabel('w_0t [deg]')
    ylabel('v_{DS2} [V]')
    axis([0 360 0 25])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,6)
    plot(vds2,id2,'.')
    xlabel('v_{DS2} [V]')
    ylabel('i_{D2} [A]')
    axis([-2 27 -0.1 1.1])
    set(gca,'xtick',[0 5 10 15 20 25])
    
    subplot(4,3,7)
    plot(deg,vin)
    xlabel('w_0t [deg]')
    ylabel('v_{IN} [V]')
    axis([0 360 -5 5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,8)
    plot(deg,vout)
    xlabel('w_0t [deg]')
    ylabel('v_{OUT} [V]')
    axis([0 360 -5 5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,9)
    plot(deg,pout)
    xlabel('w_0t [deg]')
    ylabel('p_{OUT} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    subplot(4,3,10)
    plot(deg,pd1)
    xlabel('w_0t [deg]')
    ylabel('p_{D1} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,11)
    plot(deg,pd2)
    xlabel('w_0t [deg]')
    ylabel('p_{D2} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,3,12)
    plot(deg,pdd)
    xlabel('w_0t [deg]')
    ylabel('p_{DD} [W]')
    axis([0 360 0 12])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

figure(2)
hold on
plot(Vm,Pdd,'g')
plot(Vm,Pd1,'b')
plot(Vm,Pd1,'m')
plot(Vm,Pout,'r')
plot(Vm,Pd1+Pd2,'y')
xlabel('V_m [V]')
ylabel('P_{DD} [green], P_{D1} [blue], P_{D2} [magenta], P_{OUT} [red], P_{D1}+P_{D2} [yellow], all in [W]')
axis([0 5 0 8])

figure(3)
plot(Vm,efficiency)
xlabel('V_m [V]')
ylabel('efficiency [%]')

figure(4)
plot(vin,vout)
xlabel('v_{IN}')
ylabel('v_{OUT}')