Finding n points that are equidistant around the circumference of an ellipse

Question

Hello all. This is probably a simple problem for you guys, but my geometry is a bit rusty and I am hoping that you can help.

I am trying to arrange an arbitrary number of objects around the circumference of an ellipse. My first stab at the problem resulted in the use of a simple rotational matrix like this (note that I am using graphical Cartesian coordinates where y increases down from the top of the screen).

def rotate(x, y, theta)
  x_p = (x * Math.cos(theta)) - (y * Math.sin(theta))
  y_p = (y * Math.cos(theta)) + (x * Math.sin(theta))
  return [x_p, y_p]
end

That is obviously a circle, so I use the x coordinate and feed it into the equation of my ellipse to get the y(s).

This works and does what it is supposed to, but the problem is that I am only incrementing theta to find my points, so obviously they are not equidistant around the circumference of the ellipse. That is what I am after, and any help would be greatly appreciated.

You can find a description of the problem with pictures in this thread:

http://www.bigresource.com/Tracker/Track-flash-DO1WzX6KNq/

Answer 1

This answer assumes you are interested in finding $n$ points on an ellipse such that the arc lengths between successive points are equal.

As others mentioned, this problem involves an elliptic integral, which has no elementary expression. However, many scientific computation libraries are able to compute this function numerically, and you can take advantage of that. It looks like you're using Python; one possibility here would be the GNU Scientific Library (gsl) for which Python bindings are available.

Say your ellipse is parametrized by $x(t) = a \cos t$, $y(t) = b \sin t$ where $a$ is the semimajor and $b$ the semiminor axis. Let $e = \sqrt{1-b^2/a^2}$ be the eccentricity. The arclength from $(x(0), y(0)) = (a,0)$ counterclockwise to $(x(t), y(t))$ is given by $s(t)=a E(t,e)$, where $E(t,e) = \int_0^t \sqrt{1-e^2 \sin^2\theta}d\theta$ is the "incomplete elliptic integral of the second kind". gsl implements this function as gsl_sf_ellint_E. The total arc length of your ellipse is $S = s(2\pi)$, which you can compute, so you need to find the value $t_k$ such that $s(t_k) = kS/n$, for each $k=0,\dots,n$. That is, you must solve for $t$ in the equation $s(t) = kS/n$.

But you are able to compute $s(t)$ numerically, and you also know the derivative $s'(t) = a \sqrt{1 - e^2 \sin^2 t}$ thanks to the fundamental theorem of calculus. Using this, it is a standard procedure to solve the equation iteratively using something like Newton's method. You could implement it yourself with a bit of reading, or use one of gsl's more sophisticated implementations, e.g. gsl_root_fdfsolver_steffenson.

Answer 2

I'll attempt to answer your question by misinterpreting it.

As pointed out in the comments, computing elliptic integrals is not going to be easy. But what if you wanted to find $n$ points arranged around an ellipse which form the vertices of an equilateral polygon? Now the answer to the question is given by a real algebraic variety. It's possible that this question may be computationally more tractable.

There are $2n$ variables $(x_i,y_i)$, $i=1,\ldots,n$, and $2n-1$ equations: $$ \frac{x_i^2}{a^2}+\frac{y_i^2}{b^2}=1,\ i=1,\ldots, n,$$ $$(x_i-x_{i+1})^2+(y_i-y_{i+1})^2=(x_{i+1}-x_{i+2})^2+(y_{i+1}-y_{i+2})^2, i=1,\ldots, n-1,$$ indices taken $(\mod n)$.

Also, for geometric reasons, one expects $n-1$ components to this variety, each of which is a circle. For example, if $n=5$, one would expect two solutions which are oriented in different directions, and two solutions which are star shaped. As one moves around the circle, the solution should move around.

Thus, one expects this variety to be a complete intersection defined by quadratic equations. There are methods from algebraic geometry to find solutions to such equations. There are versions of Newton's recursion which may be effective for finding a numerical solution. The dihedral symmetry might further constrain the solutions. Maybe someone could point you to some references if this sort of solution would suffice for your application?

Answer 3

The answer, with full Java code is located on StackOverflow here

Answered by:

edited Dec 11 '13 at 4:14 John Paul

answered Dec 11 '13 at 3:48 Dave

package com.math;

public class CalculatePoints {

public static void main(String[] args) {
    // TODO Auto-generated method stub

    /*
     * 
    dp(t) = sqrt( (r1*sin(t))^2 + (r2*cos(t))^2)
    circ = sum(dp(t), t=0..2*Pi step 0.0001)

    n = 20

    nextPoint = 0
    run = 0.0
    for t=0..2*Pi step 0.0001
        if n*run/circ >= nextPoint then
            set point (r1*cos(t), r2*sin(t))
            nextPoint = nextPoint + 1
        next
        run = run + dp(t)
    next
 */


    double r1 = 20.0;
    double r2 = 10.0;

    double theta = 0.0;
    double twoPi = Math.PI*2.0;
    double deltaTheta = 0.0001;
    double numIntegrals = Math.round(twoPi/deltaTheta);
    double circ=0.0;
    double dpt=0.0;

    /* integrate over the elipse to get the circumference */
    for( int i=0; i < numIntegrals; i++ ) {
        theta += i*deltaTheta;
        dpt = computeDpt( r1, r2, theta);
        circ += dpt;
    }
    System.out.println( "circumference = " + circ );

    int n=20;
    int nextPoint = 0;
    double run = 0.0;
    theta = 0.0;

    for( int i=0; i < numIntegrals; i++ ) {
        theta += deltaTheta;
        double subIntegral = n*run/circ;
        if( (int) subIntegral >= nextPoint ) {
            double x = r1 * Math.cos(theta);
            double y = r2 * Math.sin(theta);
            System.out.println( "x=" + Math.round(x) + ", y=" + Math.round(y));
            nextPoint++;
        }
        run += computeDpt(r1, r2, theta);
    }
}

static double computeDpt( double r1, double r2, double theta ) {
    double dp=0.0;

    double dpt_sin = Math.pow(r1*Math.sin(theta), 2.0);
    double dpt_cos = Math.pow( r2*Math.cos(theta), 2.0);
    dp = Math.sqrt(dpt_sin + dpt_cos);

    return dp;
}

}

Answer 4

an easier method to compute the arclength can be found under 'circumference of an ellipse on wikipedia - a good algebraic approx to the ecliptic integrals' solution envolving functions no more complex than sqrt (). I have been looking at this, it does give the correct circumferance, but I am having issues with making my incremental angle increase from the knowledge of the arc length increment! (this is probably trivial in comparison, any ideas?)

Answer 5

If a practical solution is sought For the problem of getting the euclidean distances equal then start with the approximate solution from the question, perhaps fix one point (or don't) and apply an iterative method to improve things. One possibility is a Newton's method for the sum of the squares of the differences of the squared distances.

$\sum(d_i-d_{i-1})^2$ where $d_i=(x_i-x_{i+1})^2+(y_i-y_{i+1})^2$

Perhaps Lagrange Multipliers for the minimization of $\sum d_i$ would yield a solution with points equally spaced and consecutive.