% 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;
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));
    
    beta=asin(mout);
        
    w1=c1+pattern;
    w1=w1.*(phi<beta);
    w1=w1.*(w1>0);
    w2=c2+pattern;
    w2=w2.*(phi>beta);
    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);
        
    figure(1)
    
    subplot(3,2,1)
    plot(deg,jl)
    xlabel('wt [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('wt [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('wt [deg]')
    ylabel('p_X')
    axis([0 360 -0.1 0.6])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    
    subplot(3,2,2)
    plot(deg,jin)
    xlabel('wt [deg]')
    ylabel('j_{IN}')
    axis([0 360 -0.7 0.7])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,4)
    plot(deg,min)
    xlabel('wt [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('wt [deg]')
    ylabel('p_{IN}')
    axis([0 360 -0.1 0.6])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

kk=k;
pattern=1-cos(phi)-2/pi*phi;
beta=180/pi*asin(2/pi);
gamma=180/pi*asin(2/pi);
alpha=180;

for k=kk+1:kk+nf
       
    Jout(k)=0.21+(k-kk)*0.002;
    mout=2/pi;
    Mout(k)=mout;
    
    w=Jout(k)+pattern;
    jl=w;
    jl(npoint/2+1:npoint)=w;
    mx=mx0;
    jin=jl.*s;
    pin=jin.*min;
    px=jl.*mx;
    
    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);
    
    figure(1)
    
    subplot(3,2,1)
    plot(deg,jl,'k')
    xlabel('wt [deg]')
    ylabel('j_L')
    axis([0 360 -0.7 0.7])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,3)
    plot(deg,mx,'k')
    xlabel('wt [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,'k')
    xlabel('wt [deg]')
    ylabel('p_X')
    axis([0 360 -0.1 0.6])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    
    subplot(3,2,2)
    plot(deg,jin,'k')
    xlabel('wt [deg]')
    ylabel('j_{IN}')
    axis([0 360 -0.7 0.7])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(3,2,4)
    plot(deg,min,'k')
    xlabel('wt [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,'k')
    xlabel('wt [deg]')
    ylabel('p_{IN}')
    axis([0 360 -0.1 0.6])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

% output plane
figure(2)
plot(Jout,Mout)
xlabel('J_{OUT}')
ylabel('M_{OUT}')
axis([0 0.4 0 1])

% power diagrams
figure(3)
plot(Jout,Pin,'g',Jout,Px,'b',Jout,Pout,'r')
xlabel('J_{OUT}')
ylabel('P_{IN} [green], P_X [blue], P_{OUT} [red]')
axis([0 0.4 0 0.3])

% input power, apparent power
figure(4)
plot(Jout,Pin,'b',Jout,S,'r')
xlabel('J_{OUT}')
ylabel('P_{IN} [blue], S [red]')
axis([0 0.4 0 0.3])

% power factor
figure(5)
plot(Jout,PF)
xlabel('J_{OUT}')
ylabel('PF')
axis([0 0.4 0 1])