% half-wave rectifier and L filter
% all voltages normalized by Vm
% all currents normalized by Vm/wL
% Vin=Vm sin(wt)

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;
min=sin(wt);

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

for k=1:nf
    
    mout=(k-0.5)/nf;
    Mout(k)=mout;
    
    o=sqrt(1-mout^2)+mout*asin(mout);
    
    jl=o-cos(wt)-mout*wt;
    d=(jl>0).*~((min<mout)&(wt<pi/2));
    jl=jl.*d;
    
    pin=min.*jl;
    pout=mout.*jl;
    ml=(min-mout).*(jl>0);
    pl=ml.*jl;
    
    Jout(k)=mean(jl);
    alpha(k)=360*mean(d);
    Pout(k)=mean(pout);
    Pin(k)=mean(pin);
    Jinrms(k)=sqrt(mean(jl.*jl));
    S(k)=Jinrms(k)/sqrt(2);
    PF(k)=Pin(k)/S(k);
        
    figure(1)
    
    subplot(3,2,1)
    plot(deg,jl)
    xlabel('w_0t [deg]')
    ylabel('j_L')
    axis([0 360 0 2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,3)
    plot(deg,min,'b',deg,mout,'r')
    xlabel('w_0t [deg]')
    ylabel('m_{IN} [blue], m_{OUT} [red]')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,5)
    plot(deg,pin)
    xlabel('w_0t [deg]')
    ylabel('p_{IN}')
    axis([0 360 -2 2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    
    subplot(3,2,2)
    plot(deg,ml)
    xlabel('w_0t [deg]')
    ylabel('m_L')
    axis([0 360 -1.5 1.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,4)
    plot(deg,pl)
    xlabel('w_0t [deg]')
    ylabel('p_L')
    axis([0 360 -2 2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,6)
    plot(deg,pout)
    xlabel('w_0t [deg]')
    ylabel('p_{OUT}')
    axis([0 360 -2 2])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

% power diagrams
figure(2)
plot(Mout,Pin,'b',Mout,Pout,'r')
xlabel('M_{OUT}')
ylabel('P_{IN} [blue], P_{OUT} [red]')
axis([0 1 0 0.15])

% input power, apparent power
figure(3)
plot(Mout,Pin,'b',Mout,S,'r')
xlabel('M_{OUT}')
ylabel('P_{IN} [blue], S [red]')
axis([0 1 0 1])

% power factor
figure(4)
plot(Mout,PF)
xlabel('M_{OUT}')
ylabel('PF')
axis([0 1 0 1])

% conduction angle
figure(5)
plot(Mout,alpha)
xlabel('M_{OUT}')
ylabel('conduction angle [deg]')
axis([0 1 0 360])

% regulation
figure(6)
plot(Mout,Jout)
xlabel('M_{OUT}')
ylabel('J_{OUT}')
axis([0 1 0 1])