An epicycloid is the curve a point on a circle makes as it rolls around another circle. A cyclogon is the shape a point on a regular polygon makes as it rolls across a plane. An epicyclogon is the curve traced by a point on one regular polygon as it rolls around another.

Write a program that draws an epicyclogon given r, r1, r2, n1, n2:

r = number of clockwise revolutions rolling polygon makes around stationary polygon (any real number as limited by float values) 
r1 = distance from center of stationary polygon to each of its vertices (positive real number)
r2 = distance from center of rolling polygon to each of its vertices (positive real number)
n1 = number of sides stationary polygon has (integer greater than 2)
n2 = number of sides rolling polygon has (integer greater than 2)


  • When r is negative the roller should go counterclockwise.
  • For r, one revolution occurs when the line connecting the centroids of the two shapes sweeps out a full 360 degrees. This notion is expanded to include all values of r. (So in a quarter revolution the line connecting the centroids sweeps out 90 degrees.)
  • These arguments should come from the command line or your program should prompt for them (e.g. with Python's input()).
  • r1 and r2 are relative to each other, not the dimensions of the image. So you can set one "unit" to be any number of actual pixels.

The point you must trace out is one of the vertices of the rolling shape. The shapes must start with this vertex touching a stationary vertex and two sides adjacent:

epicyclogon example

The exact starting vertices and the angle of the stationary polygon do not matter.


The output should go to an image that is at least 600x600 pixels (or some variable dimension than can be set to 600). It must show the entire epicyclogon curve specified by the parameters, well framed in the image.

The rolling and stationary polygons must also be drawn (with the roller in it's final state). The two shapes and the epicyclogon should be three noticeably different colors.

There must also be a simple way to not draw the polygons (a change of true to false in the code suffices).

Please show us at least 2 output images. It's ok to shrink them if necessary.


The shortest code that produces valid output images wins.


  • Minus 50 bytes if the output is an animated gif (or similar) of the curve being drawn.
  • Minus 150 bytes if you let n1 and n2 take the value 2 so the shapes become line segments of length 2 * r1 (or r2), "rolling" around each other. How you handle r when n1 and n2 are 2 is up to you since the centroids don't revolve around each other the way they do in other cases. (Not "rolling" at all does not count as handling it.)

Since I am quite eager to see this novel idea executed well (and it's not exactly a cakewalk), I am going to award 150 bounty rep to the winner. The contest will end the same day the bounty runs out.

The bounty will not be awarded to the winner if it is clear that they simply rewrote most of the code from another submission.

Library functions that already happen to do this (if there are any) are not allowed.

Note: This came from my leftover questions that anyone is free to post. But if no one else posts them there's a good chance I will in time. :P

  • \$\begingroup\$ I think anticlockwise should be positive instead. \$\endgroup\$ Commented Aug 24, 2014 at 5:21
  • 3
    \$\begingroup\$ @SohamChowdhury I think it hardly matters. \$\endgroup\$ Commented Aug 24, 2014 at 5:30
  • \$\begingroup\$ You're right, actually. Do you have any example images? I don't have the CDF player. \$\endgroup\$ Commented Aug 24, 2014 at 6:54
  • \$\begingroup\$ @githubphagocyte I see your point. Fixed. \$\endgroup\$ Commented Aug 24, 2014 at 16:21
  • \$\begingroup\$ @MartinBüttner Not too strict, it was just the first thing I though of. You can prompt for values in another way if necessary. \$\endgroup\$ Commented Aug 24, 2014 at 22:49

3 Answers 3


Java - 2,726 2,634 - 200 = 2434 characters

Improved from 3800 ish bytes

Thanks to all for your suggestions (especially pseudonym117), here is the new version.

I added a class P which is the point class and a class L which extends ArrayList

I also added some minor logic changes.

Here is the main class (not golfed):

import java.awt.*;
import java.awt.geom.*;

import javax.swing.*;
public class Polygons2 extends JPanel{
    public static void main(String[] args) throws InterruptedException{new Polygons2(args);}
    double q=Math.PI*2;
    int d=1;
    public Polygons2(String[] args) throws InterruptedException{
        double revolutions=Double.valueOf(args[0])*q;
        double stationaryRadius = Double.valueOf(args[1]);
        double rollingRadius = Double.valueOf(args[2]);
        int stationarySides = Integer.valueOf(args[3]);
        int rollingSides = Integer.valueOf(args[4]);    
        double dist = stationaryRadius+rollingRadius+70;
        P sp = new P(dist,dist);
        P rp = new P(sp.x,sp.y-rollingRadius-stationaryRadius);
        //get points for rolling polygon and stationary polygon
        int key=0;
        for(double stationaryAngle=-q/4;stationaryAngle<q-q/4;stationaryAngle+=q/stationarySides){
            P p=new P(Math.cos(stationaryAngle)*stationaryRadius+sp.x,Math.sin(stationaryAngle)*stationaryRadius+sp.y);
        for(double rollingAngle=q/4;rollingAngle<q+q/4;rollingAngle+=q/rollingSides){
            P p=new P(Math.cos(rollingAngle)*rollingRadius+rp.x,Math.sin(rollingAngle)*rollingRadius + rp.y);
        double g=(q/2)-((q/2-(q/rollingSides))/2) - ((q/2-(q/stationarySides))/2)-.05;
        for(P p:rollingPoints){p.r(getPoint(0), g);}
        //set up JFrame
        JFrame f = new JFrame();
        int[] pKeys= new int[]{stationaryPoints.get(0).k,rollingPoints.get(0).k};
        int index=1;
        P rc = rollingPoints.c();
        P sc =stationaryPoints.c();
        double currentRadian=Math.atan2(rc.y-sc.y,rc.x-sc.x);
        double totalRadian = 0;
            P rc2 = rollingPoints.c();
            P sc2 =stationaryPoints.c();
            double angle = Math.atan2(rc2.y-sc2.y,rc2.x-sc2.x);
            L clone=(L)path.clone();
            clone.add(new P(rollingPoints.get(1).x,rollingPoints.get(1).y));
            path = clone;
            for(P p:rollingPoints){
                int size = stationaryPoints.size();
                for(int i=0;i<size;i++){
                    P stationaryPointAtI = stationaryPoints.get(i);
                    P nextPoint=null;
                    if(p.b(stationaryPointAtI, nextPoint)==1&&containsKey(pKeys,p.k)==0){
                        //rolling point is between 2 stationary points
                    int size2=rollingPoints.size();
                    for(int h=0;h<size2;h++){
                        P nextPoint2=null;
                        if(stationaryPointAtI.b(rollingPoints.get(h), nextPoint2)==1&&containsKey(pKeys,stationaryPointAtI.k)==0){
                            //stationary point is between 2 rolling points
    volatile L path = new L();
    L rollingPoints = new L();
    L stationaryPoints = new L();
    P getPoint(int key){
        for(P p:rollingPoints){if(p.k==key){return p;}}
        for(P p:stationaryPoints){if(p.k==key){return p;}}
        return null;
    int containsKey(int[] keys,int key){
        for(int i:keys){if(key==i){return 1;}}
        return 0;
    public void paintComponent(Graphics g){
        Path2D.Double sPath = new Path2D.Double();
        sPath.moveTo(stationaryPoints.get(0).x, stationaryPoints.get(0).y);
        for(P p:stationaryPoints){
            sPath.lineTo(p.x, p.y);
        Path2D.Double rPath = new Path2D.Double();
        rPath.moveTo(rollingPoints.get(0).x, rollingPoints.get(0).y);
        for(P p:rollingPoints){
            rPath.lineTo(p.x, p.y);
        Graphics2D t = (Graphics2D)g;
        for(P p:path){g.fillOval((int)p.x-1, (int)p.y-1, 2, 2);}

And the golfed version:

import java.awt.*;import java.awt.geom.*;import javax.swing.*;import static java.lang.Math.*;class Polygons2Golfed extends JPanel{public static void main(String[]a)throws Exception{new Polygons2Golfed(a);}double q=PI*2;int d=1;public Polygons2Golfed(String[]a)throws Exception{double b,c,f;b=Double.valueOf(a[1]);c=Double.valueOf(a[2]);int d,e;d=Integer.valueOf(a[3]);e=Integer.valueOf(a[4]);f=b+c+100;P o=new P(f,f);P r=new P(o.x,o.y-c-b);int s=0;for(double u=-q/4;u<q-q/4;u+=q/d){P p=new P(cos(u)*b+o.x,sin(u)*b+o.y);p.k=s;s++;l.add(p);}for(double u=q/4;u<q+q/4;u+=q/e){P p=new P(cos(u)*c+r.x,sin(u)*c+r.y);p.k=s;s++;k.add(p);}double g=q/e/2+q/d/2-.05;for(P p:k){p.r(v(0),g);}JFrame j=new JFrame();j.add(this);j.setSize((int)f*2+60,(int)f*2+60);j.setVisible(true);m=new int[]{l.get(0).k,k.get(0).k};int ad=1;P rc=k.c();P sc=l.c();double ab,ac;ab=atan2(rc.y-sc.y,rc.x-sc.x);ac=0;while(abs(ac)<Double.valueOf(a[0])*q){P rc2=k.c();P sc2=l.c();double ah=atan2(rc2.y-sc2.y,rc2.x-sc2.x);if(ab-ah<2)ac+=(ah-ab);ab=ah;L ag=(L)n.clone();ag.add(new P(k.get(1).x,k.get(1).y));n=ag;for(P p:k){p.r(v(m[ad]),.01);int af=l.size();for(int i=0;i<af;i++){P aa=l.get(i);P w=null;if(i==af-1){w=l.get(0);}else{w=l.get(i+1);}if(p.b(aa, w)==1&&w(p.k)==0){if(ad==1)ad=0;else ad=1;m[ad]=p.k;}int ae=k.size();for(int h=0;h<ae;h++){P u=null;if(h==ae-1)u=k.get(0);else u=k.get(h+1);if(aa.b(k.get(h),u)==1&&w(aa.k)==0){if(ad==1)ad=0;else ad=1;m[ad]=aa.k;}}}}repaint();Thread.sleep(5);}}L n=new L();L k=new L();L l=new L();P v(int key){for(P p:k){if(p.k==key)return p;}for(P p:l){if(p.k==key)return p;}return null;}int[]m;int w(int key){for(int i:m){if(key==i)return 1;}return 0;}@Override public void paintComponent(Graphics g){Path2D.Double aq=new Path2D.Double();aq.moveTo(l.get(0).x,l.get(0).y);for(P p:l){aq.lineTo(p.x, p.y);}aq.closePath();Path2D.Double aw=new Path2D.Double();aw.moveTo(k.get(0).x, k.get(0).y);for(P p:k){aw.lineTo(p.x, p.y);}aw.closePath();g.setColor(Color.white);g.fillRect(0,0,getWidth(),getHeight());Graphics2D t=(Graphics2D)g;if(d==1){t.setColor(Color.black);t.draw(aq);t.setColor(Color.blue);t.draw(aw);}g.setColor(Color.green);for(P p:n){g.fillOval((int)p.x-1,(int)p.y-1,2,2);}}}

As well as classes P:

import java.awt.geom.*;class P{double x,y;public P(double a,double b){x=a;y=b;}int k;void r(P c,double g){double a,r;a=Math.atan2(y-c.y,x-c.x)+g;r=Math.sqrt((c.x-x)*(c.x-x)+(c.y-y)*(c.y-y));x=Math.cos(a)*r+c.x;y=Math.sin(a)*r+c.y;}public int b(P a,P b){if(Line2D.ptSegDist(a.x,a.y,b.x,b.y,x,y)<.5)return 1;return 0;}}

And L:

import java.util.*;public class L extends ArrayList<P>{public P c(){double x,y;x=0;y=0;for(P p:this){x+=p.x;y+=p.y;}return new P(x/size(),y/size());}}

Change int d to 0 or 1 to show polygons

arguments - 1 100 50 5 2

enter image description here

args - 1.5 100 100 7 3

enter image description here

args - 2 40 100 3 7

enter image description here

  • \$\begingroup\$ Is r really 50 in all your examples? That would mean the roller goes around 50 times. \$\endgroup\$ Commented Aug 25, 2014 at 1:58
  • \$\begingroup\$ @Calvin'sHobbies new example shows pi*3 \$\endgroup\$ Commented Aug 25, 2014 at 2:45
  • 1
    \$\begingroup\$ @StretchManiac That can't be right. 3π should take you a bit over 9 times around the stationary polygon. \$\endgroup\$ Commented Aug 25, 2014 at 20:13
  • 4
    \$\begingroup\$ It's funny how the class name is RotatingPolygonsGolfed in the "golfed" code while it's just RotatingPolygons in the normal code. ;) \$\endgroup\$ Commented Aug 25, 2014 at 20:37
  • 1
    \$\begingroup\$ you can save a good chunk of characters just by changing your imports to use * instead of specific classes... \$\endgroup\$ Commented Aug 29, 2014 at 13:34

Javascript, 1284 chars (-200 = 1084 chars)

Minified code is

function epi(B,r2,r1,n2,n1){K=Math;function C(t){return K.cos(t)}function S(t){return K.sin(t)}function A(y,x){return K.atan2(y,x)}P=K.PI;v=[[],[]];w=[[],[]];z=[];function Z(x,y,j){c=C(t=f*H+P/2);s=S(t);v[j][n]=c*x-s*y;w[j][n]=s*x+c*y;}function E(i){return{x:r1*S(t=p-i*q),y:r1*C(t)};}function D(x,y,X,Y,t){L=A(m.y,m.x);M=K.sqrt(m.x*m.x+m.y*m.y);N=K.sqrt(X*X+Y*Y);O=~~(t*(M>N?M:N)+1);for(i=J;i<=O;i++){J=1;z[n]=f*H+P+t*i/O;Z(x+M*C(T=L+t*i/O),y+M*S(T),0);Z(x+N*C(T=A(Y,X)+t*i/O),y+N*S(T),1);n++}}function F(x,y,n,r,L,s){I.strokeStyle=s;I.beginPath();for(i=0;i<n;i++)I[i?'lineTo':'moveTo'](x+r*C(t=L+(1-2*i)*P/n),y+r*S(t)*W);I.closePath();I.stroke()}p=P/n1;q=2*p;u=P/n2;H=2*u;s2=r2*S(u);g=f=l=n=J=h=0;R=300;while(l<=(B*2+1)*P/H){o=E(0);m=E(h);m.y-=o.y;m.x-=o.x;if(g<s2){D(g,-r2*C(u),-o.x,-o.y,q);h=(h+1)%n1;g+=2*r1*S(p)}else{m.x+=g-s2;D(s2,-r2*C(u),-o.x+g-s2,-o.y,H);g-=s2*2;f=(f+1)%n2;l++}}return function(e_,t,aa,W_){W=aa?-1:1;I=(e=e_).getContext('2d');I.fillStyle='black';I.fillRect(0,0,600,600);W_&1&&F(R,R,n2,r2,0,'white');T=A(w[1][0],v[1][0]);U=V=0;I.strokeStyle='teal';I.beginPath();I.moveTo(R+v[0][0],R+w[0][0]*W);while(U<t){_=A(w[1][V+1],v[1][V+1]);U+=_-T+(_+1<T?2*P:0);T=_;V++;I.lineTo(R+v[0][V],R+w[0][V]*W)}W_&2&&I.stroke();W_&4&&F(R+v[1][V],R+w[1][V]*W,n1,r1,z[V],'red')}}

Full code is

function epi( nr, r2, r1, n2, n1 ) {
function C( t )
    { return Math.cos( t ); }
function S( t )
    { return Math.sin( t ); }
function A( dy, dx )
    { return Math.atan2( dy, dx ); }

var iCCW, e, t_, xs = [[],[]], ys = [[],[]], ts = [], n = 0, iArc0 = 0;

function addpt( x, y, iBin ) {
    var c_ = C(t_ = iFrame*t2 + Math.PI/2 ),
        s_ = S(t_);

    xs[iBin][n] = c_*x-s_*y;
    ys[iBin][n] = s_*x+c_*y;

function poly1pt( iP )
    { return { x: r1*S(t_ = t1b2-iP*t1), y: r1*C(t_) }; }

function arc1( P_Arc_, xP_, yP_, xC_, yC_, t ) {
    var dx_, dy_, dxC, dyC;
    var t0 = A( dy_ = P_Arc_.y, dx_ = P_Arc_.x ),
        r_ = Math.sqrt( dx_*dx_ + dy_*dy_ ),
        t0C = A( dyC = yC_, dxC = xC_ ),
        rC = Math.sqrt( dxC*dxC + dyC*dyC ),
        nt = ~~(t*(r_>rC?r_:rC)+1);

    for( var i = iArc0; i <= nt; i++ ) {
        iArc0 = 1;
        ts[n] = iFrame*t2 + Math.PI + t*i/nt;
        addpt( xP_ + r_*C(t_ = t0+t*i/nt), yP_ + r_*S(t_), 0 );
        addpt( xP_ + rC*C(t_ = t0C+t*i/nt), yP_ + rC*S(t_), 1 );

function poly( x,y, n, r, t0, sColor ) {
    var Cx = e.getContext('2d');
    Cx.strokeStyle = sColor;
    for( var i = 0; i < n; i++ )
        Cx[i ? 'lineTo' : 'moveTo']( x + r*C(t_ = t0+(1-2*i)*Math.PI/n), y + r*S(t_)*iCCW );


var t1b2 = Math.PI/n1,
    t1 = 2*t1b2,
    t2b2 = Math.PI/n2,
    t2 = 2*t2b2,
    s1 = 2*r1*S(t1b2),
    s2 = 2*r2*S(t2b2),
    xPivot = 0,
    iPivot = 0,
    iFrame = 0,
    P_Pivot, P_Arc,
    nFrame = 0;

while( nFrame <= (nr*2+1)*Math.PI/t2 ) {
    P_Pivot = poly1pt( 0 );
    P_Arc = poly1pt( iPivot );
    if( xPivot < s2/2 ) {
        P_Arc.x -= P_Pivot.x;
        P_Arc.y -= P_Pivot.y;
        arc1( P_Arc, xPivot, -r2*C(t2b2), -P_Pivot.x, -P_Pivot.y, t1 );
        iPivot = (iPivot+1) %n1;
        xPivot += s1;
    } else {
        P_Arc.x -= (P_Pivot.x - (xPivot - s2/2));
        P_Arc.y -= P_Pivot.y;
        arc1( P_Arc, s2/2, -r2*C(t2b2), -P_Pivot.x + xPivot - s2/2, -P_Pivot.y, t2 );
        xPivot -= s2;
        iFrame = (iFrame+1) %n2;

function renderTo( eCanvas, t, isCCW, sWhat ) {
    iCCW = isCCW ? -1 : 1;
    var Cx = (e = eCanvas).getContext('2d');
    Cx.fillStyle = 'black';
    Cx.fillRect( 0,0, 600,600 );

    if( sWhat &1 )
        poly( 300,300, n2, r2, 0, 'white' );

    var tRef = A( ys[1][0], xs[1][0] ),
        tCum = 0,
        i0 = 0;

    Cx.strokeStyle = 'green';
    Cx.moveTo( 300+xs[0][0], 300+ys[0][0]*iCCW );
    while( tCum < t ) {
        t_ = A( ys[1][i0+1], xs[1][i0+1] );
        tCum += t_ - tRef + (t_ - tRef < -1 ? 2*Math.PI : 0);
        tRef = t_;
        Cx.lineTo( 300+xs[0][i0], 300+ys[0][i0]*iCCW );
    if( sWhat &2 )
    if( sWhat &4 )
        poly( 300+xs[1][i0], 300+ys[1][i0]*iCCW, n1, r1, ts[i0], 'red' );

return renderTo;

A fiddle to behold the routine in all its polygony glory (and to demonstrate animation) is found at


The script defines a function called epi that accepts the five listed parameters in the OP. epi returns a function with the signature (e,t,isCCW,flags) which accepts arguments:

  • e - a reference to a 600x600 HTML5 canvas element on which to render
  • t - the total angle (in radians) that the centroid of the second polygon should sweep around the centroid of the first. The argument passed in should not exceed 2 pi times the number of rotations passed to epi.
  • isCCW - boolean indicating whether the trace should proceed in a counterclockwise direction (as opposed to a clockwise one)
  • flags - a set of bit flags indicating which elements should be rendered
    • bit 1 - render polygon 1 if set
    • bit 2 - render trace if set
    • bit 3 - render polygon 2 if set

The function can be called any number of times with varying sets of arguments.

Some Notes:

  • The routine handles the degenerate cases where n1 = 2 and/or n2 = 2. When animating, certain combinations of lengths will cause sudden rapid advances in the trace. This is because the animation frames are indexed by the angle to the centroid of the second polygon, and d theta poly2 / d theta centroid becomes singular in cases where the centroid of 2-sided poly 2 is near a vertex of 2-sided poly 1. This does not affect the trace, however.

  • The parameter names in epi will seem confusing since throughout development, I referred to polygon 1 as "2", and polygon 2 as "1". When I realized the inconsistency between my convention and that of the OP, rather than swapping all of the indices in the code, I simply swapped the order of the arguments in epi.

  • The fiddle above imports jQuery, but this is to handle the UI. The epi function has no library dependencies.

  • The code handles CCW traces by simply inverting the Y axis. This is somewhat inelegant since polygon 2 starts in a Y-inverted position during CCW traces, but nobody said the routine had to be elegant. ;)

  • 1
    \$\begingroup\$ Beautiful! I've found the fullscreen link is easiest to work with: jsfiddle.net/7rv751jy/embedded/result \$\endgroup\$ Commented Aug 26, 2014 at 9:37
  • \$\begingroup\$ One tiny complaint is that the tracer vertex doe not start out on a stationary vertex. \$\endgroup\$ Commented Aug 26, 2014 at 9:45
  • \$\begingroup\$ Ha. I completely overlooked that in the specs. I say 'Ha' because the code was originally (inadvertently) to spec, but I changed the trace vertex because I figured the trace would look better if it started immediately. I've updated the code so that it's to spec, and updated the link to the fiddle to a fullscreen spec-compliant version. As a bonus, it knocks one character off the total count. \$\endgroup\$
    – COTO
    Commented Aug 26, 2014 at 13:22
  • \$\begingroup\$ How can I speed it up a bit? JS noob here. \$\endgroup\$ Commented Aug 29, 2014 at 0:01
  • \$\begingroup\$ @SohamChowdhury: Change the code nt = ~~(t*(r_>rC?r_:rC)+1) to nt = ~~(t*(r_>rC?r_:rC)/10+1) and it should speed things up a bit. \$\endgroup\$
    – COTO
    Commented Aug 29, 2014 at 7:33

MATLAB: 735 bytes - 200 bonus = 535

My program handles the n=2 case and draws a real-time animation. There are a few differences between the golfed and ungolfed versions:

The ungolfed version only has an option to save the animation to a file 'g.gif', by setting savegif = 1 in the code. It is off by default as it can be annoying for a few reasons:

  • Creation of an unwanted file
  • Possible lag
  • An error is generated if you have multiple monitors and the plot window is not on the right one... The gif saving had to be dropped in the golfed version as it took about 100 bytes, exceeding the size of the bonus.

The ungolfed version draws a circle on the tracer vertex. It also produces more frames and moves faster (though this can be adjusted in the golfed version by changing the numbers).


f(11,5,90,2,99,0) after program termination

golfed sample

epic(1.3,4,2,6,6,1) with gif output

ungolfed sample

Ungolfed code

%epicyclogon animation outputs to 'g.gif' if savegif=1 as well as animating in real time

function[] = epic(r,r1,r2,n1,n2,dispPoly)

savegif = 0;  %set to 1 to write .gif

cs = @(a) [cos(a);sin(a)];
vert = @(r, n, v) r * cs(2*pi*v/n);
polyPt = @(l, s, n, r) vert(r, n, floor(l/s)) + mod(l/s,1)*(vert(r, n, floor(l/s)+1) - vert(r, n, floor(l/s)));
polyPt2 = @(i, f, n, r) vert(r, n, i) + f*(vert(r, n, i+1) - vert(r, n, i));
rotm = @(a) [cos(a) -sin(a);sin(a) cos(a)];
arrpluspt = @(a, p) a + kron(p, ones(1,length(a)));
arg = @(p) atan2(p(2), p(1));

E = 1e-9;

dispPoly = dispPoly / dispPoly;

sgn = sign(-r);
r = abs(r);

s1 = 2*r1*sin(pi/n1);
s2 = 2*r2*sin(pi/n2);

%d1 = (r1*r1 - s1*s1*.25)^.5;
d2 = (r2*r2 - s2*s2*.25)^.5;

plotmax = r1+2*r2;

astep = .05; %determines amount of frames per rotation
delay = .01; % time per frame

l = 0;

lRem = 0;
lr = 0;

P1 = vert(r1, n1, 1:n1+1) * dispPoly; 
trace = [];

first = 1;
while 1

    if lr %exists while rotating about a corner of the stationary
        rotA = 2*pi/n1;
        rotA = 2*pi/n2;
    rotPt = polyPt(l, s1, n1, r1);
    lb = l + lRem;
    side1 = floor(l / s1 - E);
    side1up = side1 + lr;
    p2cen = polyPt2(side1, lb/s1 -side1 - .5 * s2/s1, n1, r1) + d2 * cs(2*pi*(side1+.5)/n1);
    if first
        p2cen0 = p2cen;
        r = r + arg(p2cen0)/(2*pi);

    for a = 0:astep:rotA    
        P2 = vert(r2, n2, 0:n2);
        P2 = rotm( pi +pi/n1 -pi/n2   +2*pi*side1/n1) * P2;
        P2 = arrpluspt(P2, p2cen);
        P2 = arrpluspt(P2, -rotPt);
        P2 = rotm(a) * P2;
        P2 = arrpluspt(P2, rotPt);
        trV = mod(floor(l/s2 + E) + lr, n2) + 1;

        cen = rotm(a) * (p2cen - rotPt) + rotPt;
        trace = [trace,P2(:,trV)]; 

        plot(P1(1,:), sgn*P1(2,:), P2(1,:)*dispPoly, sgn*P2(2,:)*dispPoly, trace(1,:),sgn*trace(2,:),P2(1,trV), sgn*P2(2,trV),'o');

        %plot(P1(1,:), P1(2,:), P2(1,:), P2(2,:), trace(1,:),trace(2,:),...
        %[0,p2cen0(1)],[0,p2cen0(2)],[0,cen(1)],[0,cen(2)], P2(1,trV), P2(2,trV),'o');

        axis square
       if savegif
           frame = getframe(1); % plot window must be on same monitor!
           img = frame2im(frame);
           [img1,img2] = rgb2ind(img,256);
       if first
           if savegif
               imwrite(img1,img2,'g','gif','DelayTime',2*delay); %control animation speed(but not really)
           first = 0;
           if savegif
               imwrite(img1,img2,'g','gif','WriteMode','append','DelayTime', 2*delay);

        adf = mod(arg(cen) - r*2*pi, 2*pi);
        if adf < astep & l/(n1*s1) + .5 > r


%cleanup for next iteration 
    jump = lRem + ~lr * s2; 
    lnex = l + jump; 

    if floor(lnex / s1 - E) > side1up 
        lnex = s1*(side1up+1);
        lRem = jump - (lnex - l);
        lr = 1;
        lRem = 0;
        lr = 0;
    l = lnex;

Golfed code

P=pi;T=2*P;F=@floor;C=@(a)[cos(a);sin(a)];g=@(i,f,n,r)r*C(T*i/n)*(1-f)+f*r*C(T*(i+1)/n);R=@(a)[C(a),C(a+P/2)];W=@(a,p)[a(1,:)+p(1);a(2,:)+p(2)];b=@(p)atan2(p(2),p(1));E=1e-9;d=d/d;S=1-2*(r>0);r=-r*S;x=2*h*sin(P/n);X=2*H*sin(P/N);M=h+2*H;l=0;z=0;L=0;A=h*C(T*(0:n)/n)*d;t=[];while 1
v=l/x;D=F(v-E);q=g(D,v-D,n,h);Z=D+L;c=g(D,v+z/x-D-.5*X/x,n,h)+H*cos(P/N)*C(T*D/n+P/n);r=r+~(l+L)*b(c)/T;for a=0:.1:T/(L*n+~L*N)


Save the function to a file with the same name, i.e. epic.m or f.m. Run it by calling the function from the Matlab console.

Usage: epic(r, r1, r2, n1, n2, dispPoly) where dispPoly is a Boolean variable (zero if false, a nonzero number if true) determining whether to draw the polygons.

Edit: Added bonus of 50 for animated image.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.