Tracing curves with epicycles

Specific domains Graphics and plotting
1

Hello there,

There was a time where epicycles were used to build real machines to reproduce (and forecast) orbits of planets. Now tracing these curves allows to create completely useless but funny animations (for the maths, see View of Tracing closed curves with epicycles). I mean, these animations do not help you to understand Fourier Series, as considering the reconstruction with graphics on the complex plane is too complicated (showing things with graphs of real functions is simpler). However, seeing these twirling circles is really funny ! Here is a Scilab implementation with metropolitan France curve (script below, uncomment line 54 to have a closer view, focused on the center of fourth circle).

q = http_get("https://raw.githubusercontent.com/gregoiredavid/france-geojson/refs/heads/master/metropole.geojson");
pos = q.geometry.coordinates(1)(1)(1:$-1,:);
N= size(pos,1);

// approximate mercator projection
pos(:,2)=(pos(:,2)-mean(pos(:,2)))*1.5;
pos(:,1)=(pos(:,1)-mean(pos(:,1)));

// compute fft of complexified curve
u=fft(pos(:,1)+%i*pos(:,2))/N;

// sort Fourier coefficients w.r.t. decreasing modulus
[Rs,ks]=gsort(abs(u));
w=ks-1;
w(w>N/2)=w(w>N/2)-N;
tspan = linspace(0,2*%pi,1000);

set(gca(),"data_bounds",[-9 9 -9 9],"tight_limits","on",...
          "background",0,"margins",[0 0 0 0],"isoview","on");
set(gcf(),"background",0,"axes_size",[800 800]);
scale = 4;
thick = 2;

for NM=[10 50 100 200 300 400]

    sref = u(ks(1:NM)).'*exp(w(1:NM)*%i*tspan);
    xref = real(sref)';
    yref = imag(sref)';
    col = [addcolor(coolwarm(NM)) 0];
    R = [Rs(1:NM);0];
    circles = [zeros(NM+1,2) 2*R 2*R zeros(NM+1,1) 360*64*ones(NM+1,1)];

    delete(gca().children)
    l = 1;
    for t = tspan
        gcf().immediate_drawing = "off"
        s = [0; cumsum(u(ks(1:NM)).*exp(w(1:NM)*%i*t))];
        x = real(s); y = imag(s);
        circles(:,1:2)=[x-R y+R];
        if t == tspan(1)
            hb = plot(x,y,"-o","marksize",4,"thickness",thick);
            hc = xarcs(circles', col);
            hc.thickness = thick;
            hr = plot(xref(1),yref(1),"r","thickness",thick)
            ht = xstring(0,1,sprintf("N = %d",NM));
            ht.font_foreground=color("white");
            ht.font_size = 4;
        else
            hb.data = [x,y];
            hc.data = circles;
            hr.data = [xref(1:l) yref(1:l)];
        end
        gca().data_bounds=[x(4)+scale*[-1 1] y(4)+scale*[-1 1]];
        gcf().immediate_drawing = "on"
        l=l+1;
    end
end
3 Likes
2

Hi there,

In fact approximating curves with the above technique becomes very funny when applied to discontinuous curves, e.g.

z(t)=\operatorname{sign}\cos(t+\pi/4)+i \operatorname{sign}\sin(t+\pi/4),\;t\geq 0

which visits sequentially each value in the set \{1+i,-1+i,-1-i,1-i\}. The following video (and code below) focuses on the mass center of circles on the second run:

t=linspace(0,2*%pi,1500)(1:$-1)'+%pi/4;
pos=sign([cos(t) sin(t)])
N= size(pos,1);

// compute fft of complexified curve
u=fft(pos(:,1)+%i*pos(:,2))/N;

// sort Fourier coefficients w.r.t. decreasing modulus
[Rs,ks]=gsort(abs(u));
w=ks-1;
w(w>N/2)=w(w>N/2)-N;
tspan = linspace(0,2*%pi,N+1);

set(gca(),"data_bounds",[-2 2 -2 2],"tight_limits","on",...
          "background",0,"margins",[0 0 0 0],"isoview","on");
set(gcf(),"background",0,"axes_size",[600 600]);
scale = 1;
thick = 2;

for zoom = 0:1
    for NM=150
        sref = u(ks(1:NM)).'*exp(w(1:NM)*%i*tspan);
        xref = real(sref)';
        yref = imag(sref)';
        col = [addcolor(coolwarm(NM)) 0];
        R = [Rs(1:NM);0];
        circles = [zeros(NM+1,2) 2*R 2*R zeros(NM+1,1) 360*64*ones(NM+1,1)];
    
        delete(gca().children)
        l = 1;
        for t = tspan
            gcf().immediate_drawing = "off"
            s = [0; cumsum(u(ks(1:NM)).*exp(w(1:NM)*%i*t))];
            x = real(s); y = imag(s);
            circles(:,1:2)=[x-R y+R];
            if t == tspan(1)
                hb = plot(x,y,"-o","marksize",4,"thickness",thick);
                hc = xarcs(circles', col);
                hc.thickness = thick;
                hr = plot(xref(1),yref(1),"r","thickness",thick)
            else
                hb.data = [x,y];
                hc.data = circles;
                hr.data = [xref(1:l) yref(1:l)];
            end
            if zoom
                gca().data_bounds=[mean(x)+scale*[-1 1] mean(y)+scale*[-1 1]];
            end
            gcf().immediate_drawing = "on"
            l=l+1;
        end
    end
    sleep(500)
end