Trinary Minor Planets

Appendix I: Definitions

Planetary and satellite orbits are defined by a minimum of six parameters, explained in the table below. Following the table is an interactive diagram.

Symbol	Parameter	Definition
	Epoch	The moment for which the values below are given.
a	Semimajor axis	Half the maximum dimension of the orbital ellipse.
e	Eccentricity	A measure of the “stretching” of the orbital ellipse. It is defined as $ \displaystyle e = \sqrt { 1 - \frac {b^2}{a^2} } $, where b is half the minimum dimension of the orbital ellipse (also known as the semi-minor axis).
i	Inclination	Tilt of the orbital ellipse with respect to a defined reference plane. In this document, this plane is that of the Earth’s equator at epoch J2000.0. Other reference planes used in astronomical literature include the ecliptic plane, defined as that of the Earth’s orbit around the Sun); the invariable plane (usually that of the Solar System, but may also be of a given planet/satellites or minor planet / satellites systems), defined as passing through the system’s barycentre and perpendicular to its angular momentum vector; or a main body’s equatorial plane.
Ω	Longitude of the ascending node	The angle, measured along the reference plane, between the reference point (usually, the vernal equinox; symbol ♈︎) and the ascending node. The ascending node is the point where the object crosses the reference plane from “below” to “above”—from south to north, usually. The symbol for the ascending node is ☊, which must not be confused with Ω; notice the shape of the bottom.
ω	Argument of the periastron	The angle, measured along the plane of the body’s orbit, between the ascending node and the periastron. The periastron (most often called “periapsis,” although this term’s etymology is wrong) is the point along the body’s orbit which is the closest to the main body. In the case of a body orbiting the Sun, this point is the perihelion; if orbiting the Earth, it’s the perigee. Sometimes, the longitude of periastron (symbol: ϖ or π) is given. It is equal to Ω + ω, even though these two angles are measured on two different planes (for historical reasons). Care must be taken not to confuse ω, ϖ, and π—and then, the latter with the mathematical constant π = 3.14159265…
And at least one of the following, which can be calculated from each other:
Tₚ	Time of periastron passage	The moment when the body is at periastron. It may be a calendar date with decimal day or day and hours, or a Julian Day (see Appendix II).
P	Orbital period	The time it takes the body to complete one full orbit. Its unit of measurement may depend on the context: for example, measuring the orbital period of a star around its galaxy’s centre in days is pointless.
ν	True anomaly	This is the angle, measured along the body’s orbit, between the periastron and the body’s current position. It is sometimes noted θ or f.
M	Mean anomaly	The fraction of the orbit’s period that has elapsed since the passing of the periastron, expressed as an angle. This angle is that of a fictitious “mean body,” that would travel at uniform speed along a fictitious circular orbit whose radius is equal to the semimajor axis of the ellipse.
n	Mean motion	The average angular speed of the body—also that of the mean fictitious body used in calculating M. Its value is $ \displaystyle n = \frac{360°}{P} $.
L	Mean longitude	The longitude at which the body would be located if its orbit were circular and free of perturbations. It does not correspond to any real angle. Its value is $ \displaystyle L = Ω + ω + M $.

Because orbits are subject to gravitational perturbations from other bodies in the system or the “parent” system (for planetary satellites, the system is that of the main planet and its satellites, whereas the “parent” system is the Solar System with the Sun, the other planets and their satellites, the minor planets, etc.), the first five elements will change slowly over long periods of time (whereas ν, M, and L change on generally much shorter timescales). A rate of change is sometimes given for these elements; for example, we have given ΔΩ and Δω for several of the minor planet satellites mentioned in the present document. These are sometimes noted instead by $ \displaystyle \dot{Ω} $ and $ \displaystyle \dot{ω} $. The epoch of the orbit must be known in order to derive accurate positions for the bodies being studied.

Because of the time-dependent nature and evolution of orbital parameters, we do not recommend using the starting values or the results of the present document to compute positions for practical purposes such as celestial observing or navigation. We recommend using accurate and up-to-date orbital elements from the JPL HORIZONS Web-Interface.

Appendix II: On the Perspective Projection of Ellipses, with Worked Examples

Converting observed positions into orbital elements and vice versa is a rather complex procedure, which is detailed in Jeremy Tatum’s Course Notes from when he was teaching at the University of Victoria, Canada. There, he was “the chief founder and is the driving force behind the [University of Victoria’s] Climenhaga Observatory’s program of astrometry of minor planets and comets, which was the only one of its kind in Canada.” Consequently, an asteroid was named in his honor: (3748) Tatum (1981 JQ)—its number is incidently just one less than that of Balam from our trinary minor planets section above.

For the following, we use the notation found on Wikipedia, not that found on Pr. Tatum’s website; the rest of the procedure is left unaltered.

Parts of the below procedures are from the author of the present document.

Equation of an Ellipse

An ellipse can described in various ways. An orbit generally being described by its semimajor axis a and its eccentricity e, we can find its semiminor axis b through $ \displaystyle e = \sqrt{1 - \frac{b^2}{a^2}} $. If it is not centered on the origin, and that it is rotated (as is usually the case for orbits), its angle of rotation may be called $ \theta $ and its centre point being $ \displaystyle \left( x_0, y_0 \right) $, we have for general form of the ellipse:

$$ \displaystyle \frac{\left[ (x - x_0) \text{ cos }\theta + (y - y_0) \text{ sin }\theta \right]^2}{a^2} + \frac{\left[ -(X - x_0) \text{ sin }\theta + (Y - y_0) \text{ cos }\theta \right]^2}{b^2} = 1 $$

Another way to describe an ellipse mathematically is by its parametric form:

$$ \displaystyle Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 $$

Conversion back and forth between these two forms is possible:

$$ \begin{array}{cl} A & = & a^2\text{ sin}^2\ \theta + b^2\text{ cos}^2\ \theta \\ \\ B & = & 2\ (b^2 - a^2)\text{ sin }\theta\text{ cos }\theta \\ \\ C & = & a^2\text{ cos}^2\ \theta + b^2\text{ sin}^2\ \theta \\ \\ D & = & -2Ax_0 - By_0 \\ \\ E & = & -Bx_0 - 2Cy_0 \\ \\ F & = & Ax_0^2 + Bx_0y_0 + Cy_0^2 - a^2b^2 \\ \\ a, b & = & \displaystyle \frac {-\sqrt{2 \left( AE^2 + CD^2 - BDE + \left( B^2 - 4AC \right) F \right) \left( \left( A + C \right) ± \sqrt{\left( A - C \right)^2 + B^2 } \right)}}{B^2 - 4AC} \\ \\ x_0 & = & \displaystyle \frac{2CD - BE}{B^2 - 4AC} \\ \\ y_0 & = & \displaystyle \frac{2AE - BD}{B^2 - 4AC} \\ \\ \theta & = & \begin{cases} \begin{array}{ll} \displaystyle \text{arctan }\left( \frac{1}{B} \left( C- A - \sqrt{\left( A - C \right)^2 + B^2} \right) \right) & \text{for }B \neq 0 \\ 0 & \text{for }B = 0, A \lt C \\ 90° & \text{for }B = 0, A \gt C \\ \end{array} \end{cases} \end{array} $$

When the equation of the ellipse is unknown, it is possible to find it with at least five points; while there is a way to solve for five equations with five unknowns with a matrix, Tatum presents “a better way.”

Let us first designate our five points as A, B, C, D, and E. Between them are imaginary straight lines AB, whose equation we’ll call α = 0; CD, to be β = 0; AC, to be γ = 0; and BD, to be δ = 0 (we’ll ignore point E for the time being). With αβ = 0 being lines AB and CD, and γδ = 0 being lines AC and BD, we have αβ + λγδ = 0, where λ is an arbitrary constant, to define any conic section that passes through A, B, C, and D. Adding E in the lot allows us to “find the value of λ that forces the equation to go through all five points.”

Example

Let us now define points A = (1, 8); B = (4, 9); C = (5, 2); D = (7, 6); and E = (8, 4). Finding α thus boils down to finding the slope and origin of AB, and so on. We thus have $ \displaystyle \mathrm{_{slope}AB} = \frac {B_y - A_y} { B_x - A_x } = \frac {9 - 8} {4 - 1} = \frac {1}{3} $, so $ \displaystyle y = \frac {x}{3} + d $. Inputting values in there, we have $ \displaystyle A_y = \frac {A_x}{3} + d \Rightarrow 8 = \frac {1}{3} + d \Rightarrow d = \frac {23}{3} $, hence $ \displaystyle y = \frac {x}{3} + \frac {23}{3} = \frac {x + 23}{3} $, which can be rewritten to $ \displaystyle \alpha = x - 3 y + 23 = 0 $. In a similar way, we find the other lines:

$$ \begin{array}{lr} \alpha = 0: & x - 3y + 23 = 0\\ \\ \beta = 0: & 2x - y - 8 = 0\\ \\ \gamma = 0: & 3x + 2y - 19 = 0\\ \\ \delta = 0: & x + y - 13 = 0 \end{array} $$

We then find the pairs of lines by multiplying two equations together. For example:

$$ \alpha\beta = (x - 3y + 23) (2x - y - 8) = 0 \Rightarrow 2x^2 - 7xy + 3y^2 + 38x + y - 184 = 0 $$

Thus:

$$ \begin{array}{lr} \alpha\beta = 0: & 2x^2 - 7xy + 3y^2 + 38x + y - 184 = 0\\ \\ \gamma\delta = 0: & 3x^2 + 5xy + 2y^2 - 58x - 45y + 247 = 0\\ \\ \alpha\beta + \lambda\gamma\delta = 0: & (2 + 3 \lambda)x^2 - (7 - 5 \lambda) xy + (3 + 2 \lambda) y^2 + (38 - 58 \lambda) x + (1 - 45 \lambda) y - 184 + 246 \lambda = 0 \end{array} $$

Now we substitute the coordinates of E in lieu of x and y to find the value of λ:

$$ \begin{array}{rl} (2 + 3 \lambda) \cdot 8^2 - (7 - 5 \lambda) \cdot 8 \cdot 4 + (3 + 2 \lambda) \cdot 4^2 + (38 - 58 \lambda) \cdot 8 + (1 - 45 \lambda) \cdot 4 - 184 + 247 \lambda & = & 0\\ \\ \displaystyle \lambda & = & \frac {76}{13} \end{array} $$

And we get the final equation of the ellipse by putting λ back in our equation:

$$ \displaystyle 508x^2 + 578xy + 382y^2 - 7827x - 6814y + 32760 = 0 $$

Coming back to our parametric equation $ Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 $, we thus have the following values:

$$ \begin{array}{cr} A & = & 508 \\ B & = & 578 \\ C & = & 382 \\ D & = & -7828 \\ E & = & -6814 \\ F & = & 32760 \end{array} $$

An additional test is needed here to know if the conic is a circle, an ellipse, a parabola, or a hyperbola, but we’ll skip it by assuming all observed object positions are indeed elliptical in nature.

We then calculate the following cofactors:

$$ \begin{array}{ccr} \bar C & = & \displaystyle AC - \left( \frac{B}{2} \right)^2 & = & 110\,535 \\ \\ \bar D & = & \displaystyle \frac{DB}{4} - \frac{AE}{2} & = & 599\,610 \\ \\ \bar E & = & \displaystyle \frac{BE}{4} - \frac{CD}{2} & = & 510\,525 \\ \end{array} $$

The centre of the ellipse is at $ \displaystyle \left( \frac {\bar E}{\bar C}, \frac {\bar D}{\bar C} \right) = \left( \frac {510\,525}{110\,535}, \frac {599\,610}{110\,535} \right) = \left( 4.61867, 5.42462 \right) $, and its major axis is inclined at an angle $ \displaystyle \text{tan }2\theta = \frac{B}{A - C} = \frac {578} {508 - 382} = \frac {578}{126} = 4.58730 \Rightarrow 2\theta = \text{atan } 4.58730 = 77.70231° \Rightarrow \theta = 38.85115° = 38°\ 51^′\ 4.15^″ $ to the y axis (or 128° 51′ 4.15″ to the x axis).

Another, less precise way is to draw the observations and the fitted ellipse, and to measure its semimajor and semiminor axes and tilt on the resulting graph with a ruler and a protractor.

It is common practice, when dealing with orbits, to set the position of the primary body as (0, 0), with North up and East corresponding to 90°.

Conversion of Apparent Orbit to Orbital Elements

We first need to establish the apparent orbit of the body. This is done by measuring, at various moments, its separation and position angle from the primary. In the case of minor planets and other Solar System objects, compensation must be made for the relative movement of the system and the Earth around the Sun; this is the subject for a possible future article. Once the apparent orbit’s equation is known (see above procedure), we can move to finding the actual orbital elements.

For this example, we’ll be using the values from Tatum’s §17.4 (while still using usual equation parameters, and not those used by Tatum):

$$ \displaystyle 14x^2 - 23xy + 18y^2 - 3x - 31y - 100 = 0 $$

So we have:

$$ \begin{array}{} A = 14\text{, }& B = -23\text{, }& C = 18\text{, }& D = -3\text{, }& E = -31\text{, }& F = -100 \end{array}$$

This gives us a centre at:

$$ \begin{array}{ccrcrcr} \bar C & = & \displaystyle 14 \cdot 18 - \left(\frac{-23}{2}\right)^2 & = & 252 - 132.25 & = & 119.75 \\ \\ \bar D & = & \displaystyle \frac{(-3)\cdot(-23)}{4} - \frac{14 \cdot (-31)}{2} & = & 17.25 - (-217) & = & 234.25 \\ \\ \bar E & = & \displaystyle \frac{(-23) \cdot (-31)}{4} - \frac{18 \cdot (-3)}{2} & = & 178.25 - (-27) & = & 205.25 \\ \\ x & = & & & \displaystyle \frac {205.25}{119.75} & = & 1.71399 \\ \\ y & = & & & \displaystyle \frac {234.25}{119.75} & = & 1.95616 \end{array} $$

The slope of the major axis is $ \displaystyle m = \frac {\bar D}{\bar E} = \frac {234.25}{205.25} = 1.14129 $; any point on the major axis of the apparent ellipse will be along the line $ \displaystyle y = 1.14129\ x $ (there is no + d here as it passes by the origin). The apparent ellipse has two points in common with the actual ellipse: periastron and apoastron. If we replace y in the parametric equation of our ellipse, we get the coordinates of these two points, by solving the following two equations:

$$ (A + Bm + Cm^2)\ x^2 + (D + Em)\ x + F = 0 $$

and

$ (C + Bn + An^2)\ y^2 + (E + Dn)\ y + F = 0 $,

where $ \displaystyle n = \frac {1}{m} = 0.87620 $. This gives us:

$$ \begin{array}{rcl} (14 + (-23) \cdot (1.14129) + 18 \cdot (1.14129)^2)\ x^2 + ((-3) + (-31) \cdot (1.14129))\ x + (-100) & = & 0 \\ \\ 11.19612\ x^2 + (-38.38002)\ x + (-100) & = & 0 \\ \\ x & = & 5.15919\text{ or }-1.73121 \end{array}$$ $$ \begin{array}{rcl} (18 + (-23) \cdot (0.87620) + 14 \cdot (0.87620)^2)\ y^2 + ((-31) + (-3) \cdot (0.87620))\ y + (-100) & = & 0 \\ \\ 8.59557\ y^2 + (-33.62860)\ y + (-100) & = & 0 \\ \\ y & = & 5.88814\text{ or }-1.97582 \end{array}$$

“If m is positive the larger solution for y corresponds to the larger solution for x; If m is negative the larger solution for y corresponds to the smaller solution for x.”

There are two ways, not detailed here, for knowing which is the periastron and which is the apoastron: Either draw the graph, or calculate the distance between each point and the origin. Here we find the periastron to be located at (-1.73121, -1.97582).

The eccentricity of the actual orbit is calculated from the ratio of distances focus–centre (or origin–centre) to periastron–centre (semimajor axis). This gives us:

$$ \displaystyle e = \frac { \sqrt { 1.71399^2 + 1.95616^2 }} { \sqrt { (1.71399 - (-1.73121))^2 + (1.95616 - (-1.97582))^2 }} = 0.49750 $$

We then need to calculate an auxiliary ellipse. The first step is calculating a lengthening factor $ \displaystyle k = \frac {1}{\sqrt{1 - e^2}} = \frac {1}{\sqrt{1 - 0.49750^2}} = 1.15279 $, then the slope of the projected latus rectum—a chord perpendicular to the major axis and passing by the focus—of the actual orbit $ \displaystyle m = \frac{-\frac{D}{2}}{\frac{E}{2}} = \frac{-\frac{-3}{2}}{\frac{-31}{2}} = -0.09677 $. We find the two points where the latus rectum intersects the apparent ellipse the same way we found the position of the apsides, except with different values of m and $ \displaystyle n = \frac{1}{m} = -10.33333 $:

$$ \begin{array}{rcl} (14 + (-23) \cdot (-0.09677) + 18 \cdot (-0.09677)^2)\ x^2 + ((-3) + (-31) \cdot (-0.09677))\ x + (-100) & = & 0 \\ \\ 16.39438\ x^2 + (-100) & = & 0 \\ \\ x & = & 2.46975\text{ or }-2.46975 \end{array}$$ $$ \begin{array}{rcl} (18 + (-23) \cdot (-10.33378) + 14 \cdot (-10.33378)^2)\ y^2 + ((-31) + (-3) \cdot (-10.33378))\ y + (-100) & = & 0 \\ \\ 1750.55556\ y^2 + (-100) & = & 0 \\ \\ y & = & 0.23901\text{ or }-0.23901 \end{array}$$

Again, “[i]f m is positive the larger solution for y corresponds to the larger solution for x; If m is negative the larger solution for y corresponds to the smaller solution for x,” so we have M = (-2.46975, 0.23901) and N = (2.46975, -0.23901).

We also need to find (at least one of) the point(s) on the apparent ellipse where a line parallel to the latus rectum crosses the centre. This is a slightly more complex set of equations, but nothing dramatic:

$$ \begin{array}{rcl} (A + Bm + Cm^2)\ x^2 + (Bd + 2Cmd + D + Em)\ x + Cd^2 + Ed + F & = & 0 \\ \\ (C + Bn + An^2)\ y^2 + (Be + 2Ane + E + Dn)\ y + Ae^2 + De + F & = & 0 \end{array} $$

where $ \displaystyle d = \bar y - m \bar x $ and $ \displaystyle e = \frac{-d}{m} $ (not to be confused with the eccentricity!), where (x̄, ȳ) are the coordinates of the centre point. We thus have:

$$ \begin{array}{rcl} d & = & 1.95616 - (-0.09677)(1.71399) & = & 2.12203 \\ \\ e & = & \displaystyle \frac{-2.12203}{-0.09677} & = & 21.92764 \\ \\ 16.39438\ x^2 + (-56.19957)\ x + (-84.72873) & = & 0 \\ \\ x & = & -1.13310 \text{ or }4.56108 \\ \\ 1750.55556\ y^2 + (-6848.73393)\ y + 6565.71866 & = & 0 \\ \\ y & = & 1.68064 \text{ or }2.23168 \end{array} $$

Because, once again, “[i]f m is positive the larger solution for y corresponds to the larger solution for x; If m is negative the larger solution for y corresponds to the smaller solution for x,” so we have K = (-1.13310, 2.23168).

Points on the auxiliary ellipse are found by multiplying segment lengths MF, FN, and KC (centre) by k to get points and segments M′F, FN′, and K′C. For M′F and FN′, this is easy enough, as MF (hence M′N′) passes by the origin, so we have $ M^′(x\text{, }y) = M(x \cdot k\text{, }y \cdot k) = (-2.84710, 0.27552) $ and $ N^′(x\text{, }y) = N(x \cdot k\text{, }y \cdot k) = (2.84710, -0.27552) $. For K′, we replace x and y by x̄ + k(x - x̄) and ȳ + k(y - ȳ), so we get $ x = 1.71399 + 1.15279\ ((-1.13310) - 1.71399) = -1.56810 $ and $ y = 1.95616 + 1.15279\ (2.23168 - 1.95616) = 2.27378 $.

We now have five points on the auxiliary ellipse, which means we can derive its parametric formula from the procedure in the previous section.

$$ \begin{array}{rcl} A & = & (+5.15919, +5.88814) \\ P & = & (-1.73121, -1.97582) \\ M^\prime & = & (-2.84709, +0.27552) \\ N^\prime & = & (+2.84709, -0.27552) \\ K^\prime & = & (-1.56810, +2.27378) \\ \end{array} $$ $$ 10.5518\ x^2 - 16.9575\ xy + 15.3528\ y^2 - 3\ x - 31\ y - 100 = 0 $$ $$ \begin{array}{rcccccl} \text{tan }2\theta & = & \displaystyle \frac{B}{A-C} & = & \displaystyle \frac{-16.9575}{10.5518-15.3528} & = & 3.53208 \\ \\ & & \theta & = & 37.09607 \\ \\ & & m & = & \text{tan }\theta & = & 0.75619 \\ \\ & & n & = & \displaystyle \frac{1}{m} & = & 1.32243 \end{array}$$

The semimajor axis of the auxiliary ellipse is its only axis that has not been foreshortened by projection. Its length is thus equal to the semimajor axis of the true orbit:

$$ \begin{array}{rcl} a & = & \frac {-\sqrt{2 \left( AE^2 + CD^2 - BDE + \left( B^2 - 4AC \right) F \right) \left( \left( A + C \right) + \sqrt{\left( A - C \right)^2 + B^2 } \right)}}{B^2 - 4AC} \\ \\ & = & 5.66544 \end{array} $$

Considering the original observations must have been in arcseconds, the answer is also in arcseconds. Once the distance to the object is known, it is simple to convert angular dimensions to physical dimensions.

The inclination of the true orbit to the plane of the sky is found by:

$$ \begin{array}{rcl} i & = & \displaystyle \text{acos } \frac {b}{a} \\ & = & \displaystyle \text{acos } \frac {2.47102}{5.66544} \\ & = & 64.14108° \end{array}$$

The longitude of the node is given by:

$$ \begin{array}{rcl} \Omega & = & \text{arctan }m \\ & = & 37.09607° \end{array}$$

from the x axis or 127.09607° from the y axis (the position angle used in astronomy). In strict mathematical terms, it is impossible to know which of the nodes this is, the ascending or the descending one, but observational data may allow answering this question—whether it’s drawing the apparent orbit and verifying the direction of motion, or using spectroscopically derived radial velocities.

Finally, we know the position of the node and of the periastron, but the longitude of the latter is not the apparent angle between the two, as it’s affected by perspective. Let’s instead call the measured angle λ. While it can be measured on a graph, it can also be calculated. Let’s first introduce angle θ:

$$ \begin{array}{rcl} \theta & = & \displaystyle \text{atan2 } \frac {\bar y}{\bar x} \\ & = & \displaystyle \text{atan } \frac {1.95616}{1.71399} \\ & = & \text{atan } 1.14129 \\ & = & 48.77511° & \text{or} & 228.77511° \end{array} $$

Deciding which value is good is done by looking at the sign of $ P_x $: if it is negative, then the second angle is good—looking at the graph would also give the answer. Angle λ is then given by subtracting the lower values of θ and Ω (the one measured from the x axis):

$$ \begin{array}{rcl} \lambda & = & \theta - \Omega \\ & = & 48.77511° - 37.09607° \\ & = & 11.67904° \end{array} $$

The penultimate piece of our puzzle is the argument of periastron, given by:

$$ \begin{array}{rcl} \text{tan } \omega & = & \text{tan } \lambda \text{ sec } i \\ \text{tan } \omega & = & \text{tan } 11.67904 \text{ sec } 64.14108 \\ \text{tan } \omega & = & 0.47393 \\ \omega & = & 25.35776° \end{array} $$

The last bit of information that is required cannot be obtained by mathematical devices but by observation alone: Determining the orbital period P and the moment of periastron passage Tₚ.

Conversion of Orbital Elements to Apparent Orbit

Orbital Elements of Iapetus
Semimajor axis	a	3564524.725648914 km
Eccentricity	e	0.02908967274787569
Ascending node	Ω	253.4568570764956°
Perichron	ω	349.4625601591317°
Inclination	i	15.75897128219843°
Mean anomaly	M	3.149688132606785°
True anomaly	ν	3.339731721713938°
Period	P	6865675.917013768 seconds
Mean motion	n	0.00005243474995781366°/d
Time of perichron	Tₚ	2459273.804760225117
Epoch		2459274.5 TDB (2021‑03‑01 00:00:00 TDB)

While it is easy to find various apparent positions of a body knowing its orbital elements, the author of the present document wanted a quick and efficient way of drawing the apparent ellipse with graphic-design software such as Adobe Illustrator, or merely coding it in SVG. The easiest solution found was to determine five random apparent points from the orbital elements, then to use the above-mentioned procedure to find an ellipse’s parametric equation from five points. We’ll use the orbital elements of Iapetus, a high-inclination moon of Saturn, for our example (table to the right). These elements, taken from the JPL Horizons Web‑Interface, are relative to Saturn’s equator and node of date for their epoch. (Note that they differ from those given, e.g., on Wikipedia, as these are for epoch J2000.0, while those we obtained from JPL Horizons are for March 1, 2021.)

The five points may be taken at random, but the author of the present document normally uses the apsides, the nodes, and the current position of the object. One way or another, we will end up with points A, B, C, D, and E and can determine the apparent ellipse from them. The procedure here is usable for any orbit seen from any position and is the same one used in the interactive diagram above. However, we do not discuss here the determination of the viewing angle, which will be left to the viewer.

At perichron (yes, that is the word for “closest point of the orbit to Saturn”), true anomaly ν = 0°. At apochron, ν = 180°. At the ascending and descending nodes, ν = (−ω) and ν = (−ω) + 180°, respectively. Finally, let’s use a different value of ν than the one in the table, just to get a more precise ellipse. We’ll use instead ν = 45°. From here, and for each of the five points, we compute the following:

$$ \begin{array}{l|r} & \text{Perichron} & \text{Apochron} & \text{Asc. node} & \text{Desc. node} & \nu = 45° \\ \hline u = \omega + \nu & 349.463 & 169.463 & 0 & 180 & 34.463 \\ x = r\ (\text{cos }\Omega\text{ cos }u - \text{sin }\Omega\text{ sin }u\text{ cos }i ) & -1552716.430 & 1645759.034 & -985898.280 & 1043950.104 & 1002563.872 \\ y = r\ (\text{sin }\Omega\text{ cos }u + \text{cos }\Omega\text{ sin }u\text{ cos }i ) & -3088186.458 & 3273238.220 & -3319157.020 & 3514596.166 & -3299306.003 \\ z = r\text{ sin }i\text{ sin }u & -171892.848 & 182193.092 & 0 & 0 & 536317.311 \end{array} $$

Where r is the radius vector of the object—its distance to the main body—and is found from:

$$ \displaystyle r = \frac {a\ (1 - e^2)}{1 + e \text{ cos }\nu} $$

Each of the $ (x, y, z) $ values (measured in the same units as a and r) must then be multiplied successively by each of the following matrices:

$$ \begin{vmatrix} \text{cos }\gamma & \text{sin }\gamma & 0 \\ -\text{sin }\gamma & \text{cos }\gamma & 0 \\ 0 & 0 & 1 \end{vmatrix} $$

where $ \gamma $ is the roll angle, i.e. the left–right angle between the reference point and the observer. It corresponds to the observer’s longitude measured in the same reference frame as the orbital elements.

$$ \begin{vmatrix} 1 & 0 & 0 \\ 0 & \text{cos }\beta & \text{sin }\beta \\ 0 & -\text{sin }\beta & \text{cos }\beta \end{vmatrix} $$

where $ \beta $ is the yaw angle, i.e. the up–down angle between the reference point and the observer. It corresponds to the observer’s latitude measured in the same reference frame as the orbital elements.

$$ \begin{vmatrix} \text{cos }\alpha & 0 & \text{sin }\alpha \\ 0 & 1 & 0 \\ -\text{sin }\alpha & 0 & \text{cos }\alpha \end{vmatrix} $$

where $ \alpha $ is the pitch angle, i.e. the apparent tilt of the reference frame’s z axis from the point of view of the observer. It corresponds to the observer-centric (geocentric) position angle of the North pole of the main body, when the reference frame is the main body’s equator, such as is the case in our example.

As a reminder, the procedure to multiply a 1 × 3 matrix (the point’s coordinates) by a 3 × 3 matrix (the transformation one) is the following:

$$ (x, y, z) \begin{vmatrix}a & b & c \\ d & e & f \\ g & h & i \end{vmatrix} = (ax + by + cz, dx + ey + fz, gx + hy + iz) $$

(Proceeding with three successive matrices allows us to avoid long, complicated equations. For example, after just the first matrix transformation, x is replaced by $ x \text{ cos }\gamma + y\ (-\text{sin }\gamma) $.)

After transformations using $ \beta = 71.47459, \alpha = 92.32287 $ (respectively, 90° minus Earth’s saturnicentric declination and 90° plus Earth’s saturnicentric right ascension, we now have the following values:

$$ \begin{array}{l|r} & \text{Perichron} & \text{Apochron} & \text{Asc. node} & \text{Desc. node} & \nu = 45° \\ \hline x & −3022716.461 & −3276470.681 & −3337229.418 & 3203845.1 & 3469396.362 \\ y & 369714.338 & 355730.0725 & 232735.8188 & −391868.5345 & −376676.2287 \\ z & −1644344.113 & −1061600.085 & 993440.7518 & 1742877.275 & 1124109.394 \end{array} $$

The z values are useless to us from now on, but may be used after the procedure is completed to know which of the nodes is the ascending and which is the descending. Using the (x, y) values to find the parametric equation of the apparent ellipse (procedure not shown), we have:

$$ 0.008473082\ x^2 + 0.127903821\ xy + 0.726626093\ y^2 + (−117.9145641)\ x + 4514.294782\ y + (−35826127649) = 0 $$

With its centre at (90564.31908, −11077.09827) and tilt angle $ \theta = 84.95072355° $ from the y axis. We have also seen how to find its dimensions, which are (432437.7422, 3563145.572). As the distance of Saturn on that day is 1,616,001,830.5 km, we can find angular distances from:

$$ \displaystyle \delta = 2\text{ arctan } \frac{dim}{2\times Dist} $$

to get an apparent ellipse of 0.126369013 × 0.007845615° centered at (0.003210982, −0.000392741°) relative to Saturn. This is illustrated in the image to the right. Note that Saturn’s angular dimension at the time is 15.38509″. The respective sizes of the orbit ellipse and Saturn’s globe (rings not shown) are proportional.

Note: We have not attempted to correct for Earth’s tilt relative to Saturn, so the orientation of the image does not match the actual orientation of Iapetus’s orbit with respect to the celestial meridian.

Appendix III: Julian Day

The Julian day system is used by astronomers to measure time in a locale- and culture-independent way. Its unit is the day, and its era (beginning moment) is noon on January 1, 4713 BCE (−4712) of the proleptic Julian calendar. This system was created by Joseph Scaliger in 1583 and is based on the calendrical solar cycle (28 years), lunar cycle (19 years), and indiction cycle (15 years). The current Julian day is .

For dates in the Gregorian calendar (used by most countries in the world), the Julian day can be calculated with the formula:

$$ \displaystyle JD = \frac{ 1461 \times \left( Y + 4800 + \frac { \left( M - 14 \right) } { 12 } \right) } { 4 } + \frac { 367 \times \left( M - 2 - 12 \times \frac {M - 14}{12} \right) } { 12 } - \frac { 3 \times \frac { Y + 4900 + \frac { M - 14 } { 12 } } { 100 } } { 4 } + D - 32075 $$

Meeus instead offers the following algorithm, valid for any date of any calendar:

Let Y be the year, M the month number (1 for January, 2 for February, etc., to 12 for December), and D the day of the month (with decimals, if any) of the given calendar date.

• If M > 2, leave Y and M unchanged.
  If M = 1 or 2, replace Y by Y − 1, and M by M + 12.
  In other words, if the date is in January or February, it is considered to be in the 13th or 14th month of the preceding year.
• In the Gregorian calendar, calculate $$ \begin{array}{llllllllllll} \displaystyle A = \text{INT }\left( \frac { Y } { 100 } \right) & & & & & & & & & & & \displaystyle B = 2 - A + \text{INT }\left( \frac { A } { 4 } \right) \end{array} $$
  In the Julian calendar, take B = 0.
• The required Julian Day is then
  $$ \displaystyle JD = \text{INT }\left( 365.25\ \left( Y + 4716 \right) \right) + \text{INT }\left( 30.6001\ \left( M + 1 \right) \right) + D + B - 1524.5 $$

Some computer software have built-in or add-in functions to compute the Julian day: In PHP, one may use the function gregoriantojd (int $month, int $day, int $year); in JavaScript, one may use the NPM package julian; in SQLite, there is a built-in julianday() function; etc.