Computing positions and eclipses of the Sun and Moon

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.

Partial solar eclipse on 21 August 2017 as seen from my dad’s house. Times are UTC.

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 apparent position of the Sun in inertial coordinates (including a function for rotating to ECEF), with an accuracy of “approximately 0.01 degree.”  Time is specified as a modified Julian date UTC.

from math import radians, sin, cos

def sun_position(mjd_utc, dt_tt=69.184):
    """ECI position (meters) of the Sun at the given MJD (UTC)."""
    # Reference: Meeus, Jean, Astronomical Algorithms (2nd Ed.). Richmond:
    # Willmann-Bell, Inc., 2009.
    T = (mjd_utc - 51544.5 + dt_tt / 86400) / 36525                  # (25.1)
    L0_deg = 280.46646 + T * (36000.76983 + T * 0.0003032)           # (25.2)
    M_deg = 357.52911 + T * (35999.05029 - T * 0.0001537)            # (25.3)
    M = radians(M_deg)
    e = 0.016708634 - T * (0.000042037 + T * 0.0000001267)           # (25.4)
    C_deg = ((1.914602 - T * (0.004817 + T * 0.000014)) * sin(M)
           + (0.019993 - T * 0.000101) * sin(2 * M)
            + 0.000289 * sin(3 * M))
    true_lon_deg = L0_deg + C_deg
    nu = radians(M_deg + C_deg)
    R = 149597870000 * 1.000001018 * (1 - e * e) / (1 + e * cos(nu)) # (25.5)
    Omega = radians(125.04 - 1934.136 * T)
    lon = radians(true_lon_deg - 0.00569 - 0.00478 * sin(Omega))
    epsilon = radians(23 + (26 + (21.448                             # (22.2)
        - T * (46.8150 + T * (0.00059 - T * 0.001813))) / 60) / 60
        + 0.00256 * cos(Omega))                                      # (25.8)
    return (R * cos(lon),                                            # (26.1)
            R * sin(lon) * cos(epsilon),
            R * sin(lon) * sin(epsilon))

def eci_to_ecef(pos, mjd_utc):
    """Convert ECI position (meters) to ECEF at the given MJD (UTC)."""
    d = mjd_utc - 51544.5                                            # (12.1)
    T = d / 36525
    theta0 = radians(280.46061837 + 360.98564736629 * d              # (12.4)
        + T * T * (0.000387933 - T / 38710000))
    s, c = sin(theta0), cos(theta0)
    return (pos[0] * c + pos[1] * s, -pos[0] * s + pos[1] * c, pos[2])

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

def moon_position(mjd_utc, dt_tt=69.184):
    """ECI position (meters) of the Moon at the given MJD (UTC)."""
    # Reference: Montenbruck and Gill, Satellite Orbits. Berlin: Springer,
    # 2005, Chapter 3.3.2.
    T = (mjd_utc - 51544.5 + dt_tt / 86400) / 36525
    L0 = radians(218.31617 + 481267.88088 * T)                       # (3.47)
    l = radians(134.96292 + 477198.86753 * T)
    lp = radians(357.52543 + 35999.04944 * T)
    F = radians(93.27283 + 483202.01873 * T)
    D = radians(297.85027 + 445267.11135 * T)
    dL = radians((22640 * sin(l) + 769 * sin(2 * l)                  # (3.48)
        - 4586 * sin(l - 2 * D) + 2370 * sin(2 * D)
        - 668 * sin(lp) - 412 * sin(2 * F)
        - 212 * sin(2 * l - 2 * D) - 206 * sin(l + lp - 2 * D)
        + 192 * sin(l + 2 * D) - 165 * sin(lp - 2 * D)
        + 148 * sin(l - lp) - 125 * sin(D)
        - 110 * sin(l + lp) - 55 * sin(2 * F - 2 * D)) / 3600)
    lon = L0 + dL;
    beta = radians((18520 * sin(F + dL + radians(                    # (3.49)
            (412 * sin(2 * F) + 541 * sin(lp)) / 3600))
        - 526 * sin(F - 2 * D) + 44 * sin(l + F - 2 * D)
        - 31 * sin(-l + F - 2 * D) - 25 * sin(-2 * l + F)
        - 23 * sin(lp + F - 2 * D) + 21 * sin(-l + F)
        + 11 * sin(-lp + F - 2 * D)) / 3600)
    r = 1000 * (385000 - 20905 * cos(l) - 3699 * cos(2 * D - l)      # (3.50)
         - 2956 * cos(2 * D) - 570 * cos(2 * l) + 246 * cos(2 * l - 2 * D)
         - 205 * cos(lp - 2 * D) - 171 * cos(l + 2 * D)
         - 152 * cos(l + lp - 2 * D))
    x = r * cos(lon) * cos(beta)                                     # (3.51)
    y = r * sin(lon) * cos(beta)
    z = r * sin(beta)
    epsilon = radians(23.43929111)                                   # (3.45)
    s, c = sin(-epsilon), cos(-epsilon)
    return (x, y * c + z * s, -y * s + z * c)

References:

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