Introduction
In this post I’m going to walk through an algorithm for calculating the positions of the four Galilean moons of Jupiter for any given moment in time. This method was developed by Jean Meeus and published in his book Astronomical Algorithms (Willmann-Bell, 2nd Ed., 1998).
All source code is checked into my GitHub repo. You can see the working software in action here: https://moonsofjupiter.azurewebsites.net. If you’re a raw beginner and want to know what you’ve just stumbled upon before diving into this post here are some good places to start:
(1) Ernie Wright provides an introduction with some cool animations that show what Galileo saw at the eyepiece and what he sketched in his own notebook.
(2) Dominic Ford has all kinds of astronomy goodness, including a nice looking corkscrew diagram for the four Galilean moons.
(3) I published an example web site to accompany this post. It would be a good idea to play around with it a bit before proceeding so you get an idea of what we’re trying to do.
A Bit of History
About an hour after dusk on 7 Jan 1610, Galileo pointed his recently built 20X refractor at Jupiter. He saw within his field of view what he thought were three background (fixed) stars very near the planet. Two of them were to the left (East) of Jupiter and the other one was to the right (West). It might have ended there. However, one thing about these stars intrigued him: all three were in a line along Jupiter’s equator. So curiosity got the better of him and at the same time the next night he looked again. And then things got interesting! These “fixed” stars had moved relative to Jupiter so that now all three were to the right of it:
I was aroused by the question of how Jupiter could be to the east of all the said fixed stars when the day before he had been to the west of two of them. I was afraid, therefore, that perhaps, contrary to the astronomical computations, his motion was direct and that, by his proper motion, he had bypassed those stars. For this reason I waited eagerly for the next night. But I was disappointed in my hope, for the sky was everywhere covered with clouds.
Galileo was struggling to understand what he was seeing through the eyepiece. He assumed at first that Jupiter moved in its orbit to bypass these fixed stars. It was cloudy on the ninth but on 10 Jan he looked again. This time he saw only two stars (the third hidden behind Jupiter) and it hit him what was going on:
[M]oving from doubt to astonishment, I found that the observed change was not in Jupiter but in the said stars. And therefore I decided that henceforth they should be observed more accurately and diligently.
Here’s a picture of Galileo’s sketchbook showing Jupiter as a large O and each “planet” as an asterisk:
For the next two months Galileo observed these new bodies, which he would later call the “Medicean planets” in honor of his patron. On Mar 1610, he quickly published his observations in a slim book called the Siderius Nuncius, or the Starry Messenger (English PDF version). Later Johannes Kepler suggested the name “satellites” in order to distinguish them from the planets around which they orbited. The name stuck. Another astronomer named Simon Marius independently observed the moons of Jupiter and gave them the names we are familiar with today: Io, Europa, Ganymede, and Callisto.
So in honor of Galileo, let’s write software that calculates the positions of the Galilean satellites for any given date. The end goal will be to calculate two sets of X and Y coordinates for each of the four moons. The X coordinate is on the x-axis in alignment with Jupiter’s equator. The Y coordinate is the vertical offset from the X-axis:
The Ephemerides
To get started we need to know the ephemerides or the velocity, trajectory, and position of both the Earth and Jupiter as they move over time relative to each other. We need a baseline from which to start and, as with other algorithms I’ve discussed on this blog, this will be the current J2000.0 epoch. See my sample project checked into GitHub for the detailed calculations. For brevity’s sake I’ll just introduce the main ones here.
Julian Day Number
Is always $$d= {JDE – 2451545.0}$$ where JDE (Julian Ephemeris Day) is “a continuous count of days and fractions thereof from the beginning of the year -4712” (Meeus, p. 59). For a complete explanation see my previous discussion of the JDE and how it is calculated. For now I’ll assume you’re comfortable with the use of the Julian Day number in astronomical calculations.
Distance from Earth to Jupiter (Δ)
The distance Δ from Earth to Jupiter, expressed in astronomical units (AU), involves calculating the radius vector r of Jupiter, the radius vector R of Earth, and the angle K derived from the equations of center for both planets. I’ll show this in more detail later when we get to the C# code. The formula is:
$$\Delta = \space \sqrt {r^2 + R^2 – 2 r R \space cos \space K}$$
Phase Angle of Jupiter (ψ)
This is an interior angle described by Earth-Jupiter-Sun, which never exceeds 11.5 because of how much further Jupiter is from the Sun relative to the Earth. We calculate it thusly:
$$sin \space \psi = {R \over Δ} \space sin \space K$$
And derive ψ by the arcsine function.
Equation of Center for Jupiter (B)
Given the mean anomaly of Jupiter N the equation of center is
$$B = 5.555 \space sin \space N + 0.168 \space sin \space 2N$$
The Four Galilean Moons
We need to calculate an angle u for each of the four Galilean moons: Io, Europa, Ganymede, and Callisto. These angles will be u1, u2, u3, and u4, respectively. Concerning these four angles Meeus writes:
For each of the four satellites, we now calculate an angle u which is measured from the inferior conjunction with Jupiter, so that u = 0° corresponds to the satellite’s inferior conjunction, u = 90° to its greatest western elongation, u = 180° to the superior conjunction, and u = 270° to the greatest eastern elongation (Meeus, p. 302).
That’s a lot to take in unless you’re an astronomer! So let’s show what he means in this top-down view:
From earth at the bottom of the diagram we look out in the night sky to Jupiter. In orbit around it are the four Galilean satellites. When a satellite is at 0 degrees (inferior conjunction) it appears to pass in front of Jupiter from our perspective. And when it is at 180 degrees (superior conjunction) it passes behind Jupiter. At 90 degrees is the satellite’s greatest western elongation and at 270 degrees is its greatest eastern elongation. In his published book Sidereus Nuncias we see his sketches from a side view as if we were looking through the telescope with him:
Notice how the greatest western elongation is labeled “Occ.” (Occident) and the greatest eastern elongation is labeled “Ori.” (Orient). These were terms used at the time for West and East, respectively.
The functions for the four angles (u1 through u4) are:
$$u_1 = 163^\text{o} 8069 + 203^\text{o} 4058646 \space \left(d \space – {\Delta \over 173}\right) + \psi \space – B$$
$$u_2 = 358^\text{o} 4140 + 101^\text{o} 2916335 \space \left(d \space – {\Delta \over 173}\right) + \psi \space – B$$
$$u_3 = 5^\text{o} 7176+ 50^\text{o} 2345180 \space \left(d \space – {\Delta \over 173}\right) + \psi \space – B$$
$$u_4 = 224^\text{o} 8092 + 21^\text{o} 4879800 \space \left(d \space – {\Delta \over 173}\right) + \psi \space – B$$
In celestial mechanics you often deal with fast-moving moons over long periods of time. As a result, the trigonometric functions used in these calculations can produce very large angles. For example, angle u4 (Callisto) on 14 Feb 2020 at midnight UTC is 158127.30105! So we must reduce these angles down to the interval 0 to 360 degrees using a modulo operation (we’ll get to that below in the source code). In this case u4 is reduced to 87°.30105.
Corrections for Perturbations
Next Meeus recommends the following corrections to produce even more accurate results:
$$G = 331^\text{o} 18 + 50^\text{o} 310482 \space \left(d \space – {\Delta \over 173}\right)$$
$$H = 87^\text{o} 45 + 21^\text{o} 569231 \space \left(d \space – {\Delta \over 173}\right)$$
The corrections G and H are needed because Ganymede and Callisto have greater orbital eccentricities. That is to say, their orbits around Jupiter are not as circular in shape as the first two interior moons but are instead more elliptical or elongated. (See my write up and solution to Kepler’s equation for a more detailed discussion of elliptical orbits.)
Next we need to determine the distance from the center of Jupiter to each of the four satellites in their position at the moment in time in which we’re interested. This will allow us to derive the coordinates of the moon on the X- and Y-axis as in the diagram I showed above this section. These are derived by:
$$r_1 = 5.9057 – 0.0244 \space \cos \space 2 \space (u_1 – u_2)$$
$$r_2 = 9.3966 – 0.0882 \space \cos \space 2 \space (u_2 – u_3)$$
$$r_3 = 14.9883 – 0.0216 \space \cos \space G$$
$$r_4 = 26.3627 – 0.1939 \space \cos \space H$$
Finally, for each of the satellites we can derive their X,Y coordinates thusly:
$$X_1 = r_1 \space \sin \space u_1$$
$$Y_1 = -r_1 \space \cos \space u_1 \space \sin \space D_E$$
And we repeat this formula for the other three as well.
C# Implementation
Take a look at the source code provided with this post. You’ll notice that there are these three main domain classes:
// assumes midnight on the date DateTime dt = new DateTime(1992, 12, 16); Moment moment = new Moment(dt); Earth earth = new Earth(moment); Jupiter jupiter = new Jupiter(moment);
Here I’m instantiating a moment in time on 16 Dec 1992 at 0h UT (midnight) and creating an Earth and Jupiter object. The ephemeris are implemented in these objects and here we calculate them for the moment in time:
// calculate lower accuracy ephemerides for Jupiter and Earth double d = moment.DayD; Angle V = jupiter.LongPeriodTerm; Angle J = jupiter.DeltaHeliocentricLongitudeEarthAndJupiter; Angle A = earth.EquationOfCenter; Angle B = jupiter.EquationOfCenter; Angle K = new Angle(J.Degrees + A.Degrees - B.Degrees); // R, r, and delta are distances expressed in astronomical units (AU) double R = earth.RadiusVector; double r = jupiter.RadiusVector; // distance from Earth to Jupiter in AU double delta = Math.Sqrt((r * r) + (R * R) - 2 * r * R * Math.Cos(K.Radians)); // phase angle of Jupiter; psi always lies between -12 and +12 degrees double psi = Math.Asin(R / delta * Math.Sin(K.Radians)); Angle psiAngle = new Angle(psi.ToDegrees()); // Jupiter's heliocentric longitude referred to the equinox of 2000.0 double lambda = 34.35 + 0.083091 * d + 0.329 * Math.Sin(V.Radians) + B.Degrees; Angle lambdaAngle = new Angle(lambda); // Jupiter's inclination of the equator on the orbital plane is 3.12 degrees double Ds = 3.12 * Math.Sin(lambdaAngle.Radians + 0.74700091985); // Jupiter's inclination of the equator on the ecliptic is 2.22 degrees double Dx = 2.22 * Math.Sin(psiAngle.Radians) * Math.Cos(lambdaAngle.Radians + 0.38397243544); // Jupiter's inclination of the orbital plane on the ecliptic is 1.30 degrees double Dy = 1.30 * ((r - delta) / delta) * Math.Sin(lambdaAngle.Radians - 1.7540558983); // planetocentric declination De of the Earth double De = Ds - Dx - Dy;
Armed with the physical ephemeris of Jupiter we can go on to calculate the positions of the four Galilean moons for the given moment in time:
// Chap 44 algorithms // correction for light time in days; see p. 298 double t = delta / 173; double u1 = 163.8069 + 203.4058646 * (d - t) + psiAngle.Degrees - B.Degrees; double u2 = 358.4140 + (101.2916335 * (d - t)) + psiAngle.Degrees - B.Degrees; double u3 = 5.7176 + (50.2345180 * (d - t)) + psiAngle.Degrees - B.Degrees; double u4 = 224.8092 + (21.4879800 * (d - t)) + psiAngle.Degrees - B.Degrees; // reduce to interval 0 to 360 u1 = u1.CorrectDegreeRange(); u2 = u2.CorrectDegreeRange(); u3 = u3.CorrectDegreeRange(); u4 = u4.CorrectDegreeRange(); double G = 331.18 + 50.310482 * (d - t); double H = 87.45 + 21.569231 * (d - t); double r1 = 5.9057 - 0.0244 * Math.Cos(2 * (u1 - u2).ToRadians()); double r2 = 9.3966 - 0.0882 * Math.Cos(2 * (u2 - u3).ToRadians()); double r3 = 14.9883 - 0.0216 * Math.Cos(G.ToRadians()); double r4 = 26.3627 - 0.1939 * Math.Cos(H.ToRadians()); // Satellite I (Io) double x1 = r1 * Math.Sin(u1.ToRadians()); double y1 = -r1 * Math.Cos(u1.ToRadians()) * Math.Sin(De.ToRadians()); // Satellite II (Europa) double x2 = r2 * Math.Sin(u2.ToRadians()); double y2 = -r2 * Math.Cos(u2.ToRadians()) * Math.Sin(De.ToRadians()); // Satellite III (Ganymede) double x3 = r3 * Math.Sin(u3.ToRadians()); double y3 = -r3 * Math.Cos(u3.ToRadians()) * Math.Sin(De.ToRadians()); // Satellite IV (Callisto) double x4 = r4 * Math.Sin(u4.ToRadians()); double y4 = -r4 * Math.Cos(u4.ToRadians()) * Math.Sin(De.ToRadians());
I rely heavily on an Angle struct to convert between degrees and radians. This is necessary because while the C# Math library requires arguments in radians, we often want to display the results in decimal degrees (or even hours, minutes, and seconds). I mentioned earlier that planetary motion generates very large angles, which must be reduced to the interval 0 to 360 degrees. I have a double extension method called CorrectDegreeRange that handles this calculation:
/// <summary>
/// For very large angles correct to put between 0 and 360 degrees
/// </summary>
/// <param name="d"></param>
/// <returns></returns>
public static double CorrectDegreeRange(this double d)
{
d %= 360;
if (d < 0)
{
d += 360;
}
return d;
}
If I run the web app that accompanies this blog post (here) for 16 Feb 2020, I get the following result:
| Date | Io (i) | Europa (e) | Ganymede (g) | Callisto (c) |
|---|---|---|---|---|
| 2/16/2020 | (-2.96, -0.16) | (7.87, 0.15) | (-14.19, -0.15) | (20.28, -0.52) |
If Galileo were alive today to sketch this result he would have drawn something like this:
* * O * *
Play around with the app yourself and see what you get. Then grab a pair of binoculars and check Jupiter around midnight to verify the results!
