% L filter, continuous 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);
mx=abs(min);
s=sign(min);

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

pattern=1-cos(phi)-2/pi*phi;
ideal=-2/(3*pi)*sin(2*wt);

for k=1:nf
       
    Jout(k)=0.21+k*0.002;
    mout=2/pi;
    Mout(k)=mout;
    
    w=Jout(k)+pattern;
    jl=w;
    jl(npoint/2+1:npoint)=w;
    jin=jl.*s;
    pin=jin.*min;
    px=jl.*mx;
    jlideal=Jout(k)+ideal;
    jinideal=jlideal.*s;
    pinideal=jinideal.*min;
    pxideal=jlideal.*mx;
    
    Pin(k)=mean(pin);
    Pinideal(k)=mean(pinideal);
    Px(k)=mean(px);
    Pxideal(k)=mean(pxideal);
    Jinrms(k)=sqrt(mean(jin.*jin));
    Jinrmsideal(k)=sqrt(mean(jinideal.*jinideal));
    S(k)=Jinrms(k)/sqrt(2);
    Sideal(k)=Jinrmsideal(k)/sqrt(2);
    PF(k)=Pin(k)/S(k);
    PFideal(k)=Pinideal(k)/Sideal(k);
    
    figure(1)
    
    subplot(3,2,1)
    plot(deg,jlideal,'m',deg,jl,'b')
    xlabel('w_0t [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,'b')
    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,pxideal,'m',deg,px,'b')
    xlabel('w_0t [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,jinideal,'m',deg,jin,'b')
    xlabel('w_0t [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,'b')
    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,pinideal,'m',deg,pin,'b')
    xlabel('w_0t [deg]')
    ylabel('p_{IN}')
    axis([0 360 -0.1 0.6])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

% input power
figure(2)
plot(Jout,Pinideal,'m',Jout,Pin,'b')
xlabel('J_{OUT}')
ylabel('P_{IN}')
axis([0.2 0.4 0 0.3])

% Px
figure(3)
plot(Jout,Pxideal,'m',Jout,Px,'b')
xlabel('J_{OUT}')
ylabel('P_X')
axis([0.2 0.4 0 0.3])

% apparent power
figure(4)
plot(Jout,Sideal,'m',Jout,S,'b')
xlabel('J_{OUT}')
ylabel('S')
axis([0.2 0.4 0 0.3])

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