% C filter, full-wave rectifier
% all voltages normalized by Vm
% all currents normalized by wCVm
% Vin=Vm cos(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=cos(wt);
mx=abs(min);
s=sign(min);
wjc=-sin(wt).*s;

nf=40; % number of frames, reduce if the simulation is too slow
Joutmax=0.25;

for k=1:nf
    
    jout=(k-0.5)/nf*Joutmax;
    Jout(k)=jout;
    
    beta=asin(jout);
    w1=sqrt(1-jout^2)-jout*(phi-beta);
    state=(phi>beta)&(w1>mx(1:npoint/2));
    state(npoint/2+1:npoint)=state;
    w1(npoint/2+1:npoint)=w1;
    
    mout=w1.*state+mx.*(~state);
    pout=jout.*mout;
    jc=wjc.*(~state)-jout*state;
    pc=mout.*jc;
    jd=jc+jout;
    px=jd.*mout;
    jin=jd.*s;
    pin=min.*jin;
    
    Pin(k)=mean(pin);
    Jinrms(k)=sqrt(mean(jin.*jin));
    S(k)=Jinrms(k)/sqrt(2);
    PF(k)=Pin(k)/S(k);
    Mout(k)=mean(mout);
    Pout(k)=jout*Mout(k);
        
    figure(1)
    
    subplot(4,2,1)
    plot(deg,mout)
    xlabel('wt [deg]')
    ylabel('m_{OUT}')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,3)
    plot(deg,jc)
    xlabel('wt [deg]')
    ylabel('j_C')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,5)
    plot(deg,jd)
    xlabel('wt [deg]')
    ylabel('j_D')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,7)
    plot(deg,jin)
    xlabel('wt [deg]')
    ylabel('j_{IN}')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    
    subplot(4,2,2)
    plot(deg,pout)
    xlabel('wt [deg]')
    ylabel('p_{OUT}')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,4)
    plot(deg,pc)
    xlabel('wt [deg]')
    ylabel('p_C')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,6)
    plot(deg,px,'r',deg,pin,'b')
    xlabel('wt [deg]')
    ylabel('p_X [red], p_{IN} [blue]')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
    subplot(4,2,8)
    plot(deg,min)
    xlabel('wt [deg]')
    ylabel('m_{IN}')
    axis([0 360 -6*Joutmax 6*Joutmax])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(1)
    
end

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

figure(3)
plot(Jout,Jinrms)
xlabel('J_{OUT}')
ylabel('J_{IN RMS}')
axis([0 Joutmax 0 2*Joutmax])

figure(4)
plot(Jout,Pin,'r',Jout,Pout,'b')
xlabel('J_{OUT}')
ylabel('P_{IN} [red], P_{OUT} [blue]')
axis([0 Joutmax 0 Joutmax])

figure(4)
plot(Jout,PF)
xlabel('J_{OUT}')
ylabel('PF')
axis([0 Joutmax 0 1])