% C filter, full-wave rectifier, comparison
% 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.2;

for k=1:nf
    
    jout=(k-0.5)/nf*Joutmax;
    Jout(k)=jout;
    
    beta=asin(jout);
    w1=sqrt(1-jout^2)-jout*(phi-beta);
    moutideal=1-jout*phi;
    moutideal(npoint/2+1:npoint)=moutideal;
    state=(phi>beta)&(w1>mx(1:npoint/2));
    w1(npoint/2+1:npoint)=w1;
    state(npoint/2+1:npoint)=state;
    mout=w1.*state+mx.*(~state);
    
    Mout(k)=mean(mout);
    Moutideal(k)=mean(moutideal);
    
    figure(1)
    
    plot(deg,mout,'b',deg,moutideal,'m')
    xlabel('wt [deg]')
    ylabel('m_{OUT} [blue], m_{OUT ideal} [magenta]')
    axis([0 360 0 1.1])
    set(gca,'xtick',[0 60 120 180 240 300 360])
        
    pause(0.1)
    
end

figure(2)
plot(Jout,Mout,'b',Jout,Moutideal,'m')
xlabel('Jout')
ylabel('M_{OUT} [blue], M_{OUT ideal} [magenta]')
axis([0 Joutmax 0 1.1])