% L filter, discontinuous conduction mode
% 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;
phi=wt(1:npoint/2);
deg=wt*180/pi;
min=sin(wt);
mx0=abs(min);
s=sign(min);

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

for k=1:nf
    
    mout=1-(k-0.5)/nf*(1-2/pi);
    Mout(k)=mout;
    
    pattern=-cos(phi)-mout*phi;
    
    c1=2+sqrt(1-mout^2)+mout*(asin(mout)-pi);
    c2=sqrt(1-mout^2)+mout*(asin(mout));
    
    gamma=asin(mout);
        
    w1=c1+pattern;
    w1=w1.*(phi<gamma);
    w1=w1.*(w1>0);
    w2=c2+pattern;
    w2=w2.*(phi>gamma);
    w2=w2.*(w2>0);
    w=w1+w2;
    
    jl=w;
    jl(npoint/2+1:npoint)=w;
    state=jl>0;
    mx=mx0.*state+mout*(~state);
    jin=jl.*s;
    pin=jin.*min;
    px=jl.*mx;
    
    Jout(k)=mean(jl);
    Pout(k)=mout*Jout(k);
    Pin(k)=mean(pin);
    Px(k)=mean(px);
    Jinrms(k)=sqrt(mean(jin.*jin));
    S(k)=Jinrms(k)/sqrt(2);
    PF(k)=Pin(k)/S(k);
    a(k)=180*mean(state);
    g(k)=180/pi*gamma;
    b(k)=180*((g(k)+a(k))/180-floor((g(k)+a(k))/180));
        
    figure(1)
    
    subplot(3,2,1)
    plot(deg,jl)
    xlabel('w_0t [deg]')
    ylabel('j_L')
    axis([0 360 -0.5 0.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,3)
    plot(deg,mx)
    xlabel('w_0t [deg]')
    ylabel('m_X')
    axis([0 360 -1.1 1.1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,5)
    plot(deg,px)
    xlabel('w_0t [deg]')
    ylabel('p_X')
    axis([0 360 -0.1 0.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    
    subplot(3,2,2)
    plot(deg,jin)
    xlabel('w_0t [deg]')
    ylabel('j_{IN}')
    axis([0 360 -0.5 0.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,4)
    plot(deg,min)
    xlabel('w_0t [deg]')
    ylabel('m_{IN}')
    axis([0 360 -1.1 1.1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,6)
    plot(deg,pin)
    xlabel('w_0t [deg]')
    ylabel('p_{IN}')
    axis([0 360 -0.1 0.5])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

figure(2)
plot(Mout,Pin,'g',Mout,Px,'b',Mout,Pout,'r')
xlabel('M_{OUT}')
ylabel('P_{IN} [green], P_X [blue], P_{OUT} [red]')
axis([0.6 1 0 0.25])

figure(3)
plot(Mout,a,'r',Mout,b,'b',Mout,g,'g')
xlabel('M_{OUT}')
ylabel('alpha [red], beta [blue], gamma [green], all in [deg]')
axis([0.6 1 0 180])
set(gca,'ytick',[0 30 60 90 120 150 180])

figure(4)
plot(Mout,Pin,'b',Mout,S,'r')
xlabel('M_{OUT}')
ylabel('P_{IN} [blue], S [red]')
axis([0.6 1 0 0.25])

figure(5)
plot(Mout,PF)
xlabel('M_{OUT}')
ylabel('PF')
axis([0.6 1 0 1])

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

figure(7)
plot(Jout,Mout)
xlabel('J_{OUT}')
ylabel('M_{OUT}')
axis([0 0.25 0.6 1])