WaBis

walter.bislins.ch

Ellipse Connection between Focal Point Angle and Area

To calculate an Analemma I need an equation that gives me the angle from the Sun to the Earth at a specific time, taking the ellipsoidal orbit into account. According to Kepler's laws of planetary motion, a line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time.

So the function I need is one that inputs the n · 365.25th part of the whole swept area of a year for the day n and outputs the angle: φ = f( n · Ao / 365.25 ) where Ao is the area of the ellipse. As we will see, this is not a trivial function, but in effect can only be solved numerically. The inverse function A = f(φ) however has an analytic solution than can be used with a Newton-Solver to get the desired result.

I will derive all the equations and present a JavaScript to calculate φ = f(A) for any ellipse.

Distance to Ellipse Periphery

Ellipse Labels

I will assume that the coordinate system has its origin at the focal point of the ellipse, where the sun is located. The x axis points towards the nearest point of the ellipse, the periapsis.

The distance r from the focal point F1 to a point P of the ellipse as a function of the angle φ is:

(1)
with
Source

Wikipedia: Ellipsengleichung

where'
' =' 'distance from one focal point to the periphery of the ellipse at the angle
' =' 'angle from the X axis counter clockwise
' =' 'semi-major axis
' =' 'semi-minor axis
' =' 'linear excentricity

Ellipse Curve

A point on the ellipse is described by the following:

(2)
where'
' =' 'point coordinates; origin is at the right focal point of the ellipse
' =' 'radial distance from focal point F1 to P, see (1)
' =' 'angle measured from the x axis

Area Segment of Ellipse

To get the area of a segment from the x axis to an angle we have to define the area of an infinitesimal small triangle with height and top angle . Then we can add all such triangles from 0 to with an integral.

The infinitesimal small triangle area is:

(3)
where'
' =' 'infinitesimal small segment area at angle
' =' 'height of segment triangle
' =' ' = base lenght of segment triangle
' =' 'radial distance from focal point to P, see (1)
' =' 'angle measured from the x axis

The total area of the segment from the x axis to the angle is then the integral:

(4)
where'
' =' 'area of segment from 0 to
' =' 'radial distance from focal point F1 to P at angle , see (1)
' =' 'angle measured from the x axis

This integral can be solved analytically. I used Geogebra, see Media ellipse area - A(angle).zip Info. After restricting the angle to and simplifying the result from Geogebra we get the area as a function of the angle:

(5)
with
where'
' =' 'area of segment from 0 to
' =' 'full area of ellipse
' =' 'angle measured from the x axis
' =' 'semi-major axis
' =' 'semi-minor axis
' =' ' = linear excentricity

There is a numerical problem however. For φ = π the tan(φ/2) is undefined. So for φ near π we can not solve (5) numerically. We have to find at least a linear approximation around this point. From the first derivative of (5) we can get the slope at φ = π and the function value at this point is exactly half the ellipse area. So the linear approximation is:

(6)
where'
' =' 'a very small number, e.g. 10−13 · π for calculators with at least 14 to 15 digits precision

So for |φπ| < ε, where ε is a very small number, we have to use the approximation (6) instead of (5).

Here is what the area A as a function of the angle for the shown ellipse looks like:

N: Number of sub-areas of equal size Ao / N.

JavaScript

The following JavaScript code implements a class Ellipse which provides functions to get the area from an angle at the focal point F1 and a function to get the focal point angle for a partial ellipse area, plus some handy additional functions to draw the ellipse.

function Ellipse( a, b ) {
  // Constructor for Ellipse.
  // a = semi-major radius
  // b = semi-minor radius

  this.epsilon = 1e-13;
  this.NewtonSolverData = { MaxLoops: 1000, Result: 0, Status: '', NLoops: 0 };

  this.Size( a, b );
}

Ellipse.prototype.Size = function( a, b ) {
  // Resizes the Ellipse and calculates some usefull values from a and b

  this.a = a;
  this.b = b;
  this.e2 = a*a - b*b;
  this.e = Math.sqrt( this.e2 );
  this.a2 = a*a;
  this.b2 = b*b;
}

Ellipse.prototype.FocalRadius = function( angle ) {
  // Returns the radius from the focal point to a point on the ellipse.
  // angle = angle from the positive x axis in radian.

  return (this.a2 - this.e2) / (this.a + this.e * Math.cos(angle));
}

Ellipse.prototype.Point = function( angle ) {
  // Returns the ellipse point [ x, y ].
  // angle = angle from the positive x axis in radian.

  var r = this.FocalRadius( angle );
  var x = r * Math.cos( angle );
  var y = r * Math.sin( angle );
  return [ x, y ];
}

Ellipse.prototype.Area = function( ) {
  // Returns the area of the ellipse.

  return this.a * this.b * Math.PI;
}

Ellipse.prototype.PartialArea = function( angle ) {
  // Returns the partial area of the ellipse between 0 and angle.
  // angle = angle from the positive x axis in radian.
  // Requires 0 <= angle <= 2 pi

  if (Math.abs( angle - Math.PI ) >= this.epsilon) {
    var K = angle <= Math.PI ? 0 : 1;
    var A = this.Area();
    var tan_2 = Math.tan( angle / 2 );
    var term1 = this.a * this.b * Math.atan( (this.a - this.e) * tan_2 / this.b );
    var term2 = (this.e * this.b2 * tan_2) / ((this.a - this.e) * tan_2*tan_2 + this.a + this.e);
    return K * A + term1 - term2;
  } else {
    // because calculation above is undefined for angle around pi, we have to use this linear approx
    var A = this.Area() / 2;
    var ae = this.a + this.e;
    return (ae*ae/2) * (angle - Math.PI) + A;
  }
}

Ellipse.prototype.FocalAngleFromArea = function( areaFraction ) {
  // Returns the angle that a partial area occupies. 
  // This is the inverse function of PartialArea(angle).
  // areaFraction = fraction of the whole ellipse area.
  // Requires 0 <= areaFraction <= 1

  var area = this.Area();
  var areaPart = area * areaFraction;
  var me = this; // closure
  var angle = this.SolveWithNewton( 
    function(phi) { return me.PartialArea(phi); }, 
    areaPart, Math.PI, this.epsilon 
  );
  return angle;
}

Ellipse.prototype.SolveWithNewton = function( fn, target, guess, prec ) {
  // function fn( x )
  // In this.NewtonSolverData you can specify properties:
  //   { MaxLoops: 1000 }
  // Solver returns in aData:
  //   { Result, Status, NLoops }

  var dx, X, Y, Y1, Y2, Xnew, slope;

  data = this.NewtonSolverData;
  data.Result = 0;
  data.Status = '';
  data.NLoops = 0;
  var maxLoops = data.MaxLoops;
  var nLoops = 0;
  var eps = prec / 2.0;
  var eps2 = eps / 2;
  Xnew = guess;
  do {
    X = Xnew;
    try {
      Y = fn( X, data ) - target;
      if (Y == 0) break; // exact solution found

      Y1 = fn( X-eps2, data ) - target;
      if (Y1 == 0) { // exact solution found
        Xnew = X - eps2;
        break;
      }

      Y2 = fn( X + eps2, data ) - target;
      if (Y2 == 0) { // exact solution found
        Xnew = X + eps2;
        break;
      }

      slope = (Y2 - Y1) / eps;
      if (slope == 0) { // slope zero, cant find a solution
        data.Status = 'extremum found';
        break;
      }

      Xnew = X - (Y / slope);

    } catch(err) {
      data.Status = err.message;
      break;
    }

    nLoops++;
    dx = Math.abs( Xnew - X );

  } while ((dx > prec) && (nLoops < maxLoops));

  if (data.Status == '' && nLoops >= maxLoops) data.Status = 'max loops exceedet';

  data.NLoops = nLoops;
  data.Result = (data.Status == '') ? Xnew : 0;

  return data.Result;
}

More Page Infos / Sitemap
Created Tuesday, September 3, 2019
Scroll to Top of Page
Changed Tuesday, December 1, 2020