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.
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 | ||||||||||||||||
where' |
|
A point on the ellipse is described by the following:
(2) |
| |||||||||
where' |
|
To get the area of a segment from the x axis to an angle
The infinitesimal small triangle area is:
(3) | ||||||||||||||||
where' |
|
The total area of the segment from the x axis to the angle
(4) |
| |||||||||
where' |
|
This integral can be solved analytically. I used Geogebra, see ellipse area - A(angle).zip . After restricting the angle to
(5) |
| ||||||||||||||||||
with | |||||||||||||||||||
where' |
|
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' |
|
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
N: Number of sub-areas of equal size Ao / N.
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; }