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.

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 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:

  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.
This entry was posted in Uncategorized. Bookmark the permalink.

4 Responses to Computing positions and eclipses of the Sun and Moon

  1. Pete says:

    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.

  2. Fred Barclay says:

    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!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s