Looking ahead to this new year, some of us here in the United States will have several opportunities to see the Earth and Moon get in each other’s way, so to speak. There will be a lunar eclipse on 14-15 April 2014, and another in October. There will also be a partial solar eclipse on 23 October 2014… and a *total* solar eclipse three years from now, on 21 August 2017.

For example, the following simple animation shows what the 2017 eclipse will look like from my dad’s house, which is slightly south of the path of totality.

Predicting eclipses requires knowing where the Sun and Moon will be in relation to the Earth at any given time. (Actually, I guess it really only requires knowing which web sites have prediction/animation tools. NASA has a good site for predicting eclipses at arbitrary latitude/longitude locations, but without corresponding visualizations. And timeanddate.com has some nice animated visualizations of eclipses, but only for canned locations.)

Computing the location of the Sun and/or Moon presents some interesting challenges. In particular, the complexity of the solution depends a lot on how much accuracy is required. The references at the end of this post provide equations with accuracies that are suitable for the kinds of applications I am interested in… the trick is simply stitching together those equations, coordinate frames, etc.

The result is a couple of functions that can be handy in many different exercises, from recreational to educational to work-related: predicting eclipses, predicting the visible illumination of satellite passes (e.g., the International Space Station), evaluating the accuracy of trying to use an analog watch as a compass, and extending the accuracy of ballistically propagated satellite orbits.

First, following is Python code for computing the position of the Sun, in ECEF coordinates, with an accuracy of “approximately 0.01 degree.” Time is specified in seconds since 1970-01-01T00:00:00Z.

import math def sun_position(t): """Return ECEF position (m) of the sun at time t (sec. since 1970). """ # Reference: Meeus, Jean, Astronomical Algorithms (2nd Ed.). Richmond: # Willmann-Bell, Inc., 2009, Ch. 25. d = t / 86400.0 - 10957.5 T = d / 36525 L = math.radians(280.46646 + T * (36000.76983 + 0.0003032 * T)) M = math.radians(357.52911 + T * (35999.05029 - 0.0001537 * T)) e = 0.016708634 - T * (0.000042037 + 0.0000001267 * T) C = math.radians((1.914602 - T * (0.004817 + 0.000014 * T)) * math.sin(M) + (0.019993 - 0.000101 * T) * math.sin(2 * M) + 0.000289 * math.sin(3 * M)) s = L + C v = M + C R = 149597870000 * (1.000001018 * (1 - e * e) / (1 + e * math.cos(v))) omega = math.radians(125.04 - 1934.136 * T) lon = math.radians(math.degrees(s) - 0.00569 - 0.00478 * math.sin(omega)) epsilon = math.radians(23.0 + (26.0 + (21.448 - T * (46.8150 + T * (0.00059 - 0.001813 * T))) / 60) / 60 + 0.00256 * math.cos(omega)) x = R * math.cos(lon) y = R * math.sin(lon) * math.cos(epsilon) z = R * math.sin(lon) * math.sin(epsilon) theta = math.radians(280.46061837 + 360.98564736629 * d + T * T * (0.000387933 - T / 38710000)) s = math.sin(theta) c = math.cos(theta) return (x * c + y * s, -x * s + y * c, z)

Next is the corresponding function for the Moon, with an accuracy of “several arcminutes and about 500 kilometers in the lunar distance.”

import math def moon_position(t): """Return ECEF position (m) of the moon at time t (sec. since 1970). """ # Reference: Montenbruck and Gill, Satellite Orbits. Berlin: Springer, # 2005, Ch. 3.3.2. d2000 = t / 86400.0 - 10957.5 T = d2000 / 36525 L0 = math.radians(218.31617 + 481267.88088 * T - 1.3972 * T) l = math.radians(134.96292 + 477198.86753 * T) lp = math.radians(357.52543 + 35999.04944 * T) F = math.radians(93.27283 + 483202.01873 * T) D = math.radians(297.85027 + 445267.11135 * T) # Longitude with respect to the equinox and ecliptic of the year 2000. lon = L0 + math.radians( (22640 * math.sin(l) + 769 * math.sin(2 * l) - 4586 * math.sin(l - 2 * D) + 2370 * math.sin(2 * D) - 668 * math.sin(lp) - 412 * math.sin(2 * F) - 212 * math.sin(2 * l - 2 * D) - 206 * math.sin(l + lp - 2 * D) + 192 * math.sin(l + 2 * D) - 165 * math.sin(lp - 2 * D) + 148 * math.sin(l - lp) - 125 * math.sin(D) - 110 * math.sin(l + lp) - 55 * math.sin(2 * F - 2 * D)) / 3600) # Lunar latitude. beta = math.radians( (18520 * math.sin(F + lon - L0 + math.radians( (412 * math.sin(2 * F) + 541 * math.sin(lp)) / 3600)) - 526 * math.sin(F - 2 * D) + 44 * math.sin(l + F - 2 * D) - 31 * math.sin(-l + F - 2 * D) - 25 * math.sin(-2 * l + F) - 23 * math.sin(lp + F - 2 * D) + 21 * math.sin(-l + F) + 11 * math.sin(-lp + F - 2 * D)) / 3600) # Distance from the center of the Earth. r = (385000 - 20905 * math.cos(l) - 3699 * math.cos(2 * D - l) - 2956 * math.cos(2 * D) - 570 * math.cos(2 * l) + 246 * math.cos(2 * l - 2 * D) - 205 * math.cos(lp - 2 * D) - 171 * math.cos(l + 2 * D) - 152 * math.cos(l + lp - 2 * D)) * 1000 # Convert ecliptic to equatorial Cartesian coordinates. lon = lon + math.radians(1.3972 * T) x = r * math.cos(lon) * math.cos(beta) y = r * math.sin(lon) * math.cos(beta) z = r * math.sin(beta) epsilon = math.radians(23.43929111) s = math.sin(-epsilon) c = math.cos(-epsilon) y, z = (y * c + z * s, -y * s + z * c) # Convert to ECEF. theta = math.radians(280.46061837 + 360.98564736629 * d2000) s = math.sin(theta) c = math.cos(theta) return (x * c + y * s, -x * s + y * c, z)

**References:**

- Meeus, Jean, Astronomical Algorithms (2nd Ed.). Richmond: Willmann-Bell, Inc., 2009, Ch. 25.
- Montenbruck and Gill, Satellite Orbits. Berlin: Springer, 2005, Ch. 3.3.2.

Where does 0.00256 factor come from in the variable epsilon?

The epsilon value is the obliquity of the ecliptic, and the 0.00256*cos(omega) term is from equation (25.8) in Meeus, used to correct for the *apparent* position of the Sun.

Would it be alright to use some or all of this code in a project I’m working on? If so, are there any restrictions in case I publish the source code (I’d probably use GPL v2 or MIT license).

Thanks!

No problem, the source code may be considered in the public domain, although be sure to reference/credit Meeus and Montenbruck/Gill for the math heavy-lifting.

Anything become of this project? I’ve an application for this sort of code, but I’m far too innumerate to actually figure it out myself :-).