Eclipse Functionality

[Bernmeister]

I noticed in the file TODO.rst adding the functionality for lunar/solar eclipses.

I am not sure if you intend to calculate this information or rather draw from a table. For what it's worth, below is a code sample I use in my application to extract eclipse information from tables I took from NASA. If any of the code is useful, please use it in Skyfield.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


# Lunar/solar eclipse tables up to 2099.
# Courtesy of Fred Espenak and Jean Meeus.


import datetime


# Eclipse Types.
ECLIPSE_TYPE_ANNULAR = "A"
ECLIPSE_TYPE_HYBRID = "H"
ECLIPSE_TYPE_PARTIAL = "P"
ECLIPSE_TYPE_PENUMBRAL = "N"
ECLIPSE_TYPE_TOTAL = "T"


# Gets the upcoming eclipse, lunar or solar.
#
# dateTimeUTC - Current date/time in UTC as a python datetime.datetime object.
# isLunar - If True, finds the next lunar eclipse; otherwise finds the next solar eclipse.
#
# Returns a tuple of strings describing the next eclipse:
#    dateTime (format of YYYY-MM-DD HH:MM:SS)
#    eclipseType ("A" = Annular, "H" = Hybrid, "N" = Penumbral, "P" = Partial, "T" = Total)
#    latitude (south is negative)
#    longitude (east is negative).
def getEclipse( dateTimeUTC, isLunar ):
    if isLunar:
        eclipseData = __lunarEclipseData

    else:
        eclipseData = __solarEclipseData

    eclipseInfo = None
    for eclipse in eclipseData:
        dateTime = datetime.datetime.strptime( eclipse[ 0 ] + ", " + eclipse[ 1 ] + ", " + eclipse[ 2 ] + ", " + eclipse[ 3 ], "%Y, %m, %d, %H:%M:%S" )
        dateTime = dateTime - datetime.timedelta( seconds = int( eclipse[ 4 ] ) ) # Need to subtract delta T (https://eclipse.gsfc.nasa.gov/LEcat5/deltat.html).
        if dateTimeUTC <= dateTime:
            latitude = eclipse[ 6 ][ 0 : len( eclipse[ 6 ] ) - 1 ]
            if eclipse[ 6 ][ -1 ] == "S":
                latitude = "-" + latitude

            longitude = eclipse[ 7 ][ 0 : len( eclipse[ 7 ] ) - 1 ]
            if eclipse[ 7 ][ -1 ] == "E":
                longitude = "-" + longitude

            eclipseInfo = ( str( dateTime ), eclipse[ 5 ], latitude, longitude )
            break

    return eclipseInfo


# https://eclipse.gsfc.nasa.gov/5MCLE/5MKLEcatalog.txt
#
# Title: Five Millennium Catalog of Lunar Eclipses: -1999 to +3000 (2000 BCE to 3000CE)
# Authors: Fred Espenak and Jean Meeus
# Source: NASA Technical Publication NASA/TP-2009-214173
# On Web: https://eclipse.gsfc.nasa.gov/SEpubs/5MKLE.html
# Catalog Key: https://eclipse.gsfc.nasa.gov/LEcat5/LEcatkey.html
# Date: 2011 May 23
#
#      Year   Month  Day   HH:MM:SS    DT   Type  Lat   Long
__lunarEclipseData = [
    [ "2020", "11", "30", "09:44:01", "72", "N", "21N", "148W" ],
    [ "2021", "05", "26", "11:19:53", "72", "T", "21S", "170W" ],
    [ "2021", "11", "19", "09:04:06", "73", "P", "19N", "139W" ],
    # <snip>
    [ "2099", "09", "29", "10:36:38", "202", "N", "3N", "161W" ] ]


# https://eclipse.gsfc.nasa.gov/SEpubs/5MKSE.html
# https://eclipse.gsfc.nasa.gov/5MCSE/5MKSEcatalog.txt
#
# Title: Five Millennium Catalog of Solar Eclipses: -1999 to +3000 (2000 BCE to 3000CE)
# Authors: Fred Espenak and Jean Meeus
# Source: NASA Technical Publication NASA/TP-2008-214170
# Catalog Key: https://sunearth.gsfc.nasa.gov/eclipse/SEcat5/catkey.html
# Date: 2008 10 07
#
#      Year   Month  Day   HH:MM:SS    DT   Type  Lat      Long
__solarEclipseData = [
    [ "2020", "12", "14", "16:14:39", "72", "T", "40.3S", "67.9W" ],
    [ "2021", "06", "10", "10:43:07", "72", "A", "80.8N", "66.8W" ],
    [ "2021", "12", "04", "07:34:38", "73", "T", "76.8S", "46.2W" ],
    # <snip>
    [ "2099", "09", "14", "16:57:53", "202", "T", "23.4N", "62.8W" ] ]

[brandon-rhodes]

Thanks for the sample code! It might help folks who find this issue from a Google search and are interested in reading the tables themselves, and also having this issue open can replace that old todo item.

The goal of the todo item, however, is for Skyfield to compute the eclipses itself, not to draw from a table; so while your code will be useful to folks who need it in the meantime, it's not how I expect the Skyfield internals will ultimately look.

[Bernmeister]

I'll have a crack at it...

[Bernmeister]

I've had a bash at calculating lunar/solar eclipse date/time. Outstanding is determining the eclipse type...need help with that one. See code below; free to use in Skyfield or otherwise.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


from enum import Enum

from skyfield import almanac
from skyfield.api import load

import calendar, datetime, math


class EclipseType( Enum ):
    CENTRAL = 0
    NONE = 1
    NOT_CENTRAL = 2
    TOTAL = 3
    ANNULAR = 4
    ANNULAR_TOTAL = 5
    PARTIAL = 6


def getMoonDatesForOneYear( utcNow, timeScale, ephemeris ):
    t0 = timeScale.utc( utcNow.year, utcNow.month, utcNow.day )
    t1 = timeScale.utc( utcNow.year, utcNow.month, utcNow.day + 366 )
    t, y = almanac.find_discrete( t0, t1, almanac.moon_phases( ephemeris ) )
    moonPhases = [ almanac.MOON_PHASES[ yi ] for yi in y ]
    moonPhaseDateTimes = t.utc_datetime()
    nextNewMoonDateTimes = moonPhaseDateTimes[ moonPhases.index( almanac.MOON_PHASES[ 0 ] ) : : 4 ].tolist() # New moon.
    nextFullMoonDateTimes = moonPhaseDateTimes[ moonPhases.index( almanac.MOON_PHASES[ 2 ] ) : : 4 ].tolist() # Full moon.
    return nextFullMoonDateTimes, nextNewMoonDateTimes


# Chapter 5
# Practical Astronomy with your Calculator or Spreadsheet, Fourth Edition
# Peter Duffett-Smith, Jonathan Zwart
def julianDayToGregorian( julianDay ):
    F, Z = math.modf( julianDay + 0.5 )
    if Z < 2299161:
        A = Z

    else:
        alpha = int( ( Z - 1867216.25 ) / 36524.25 )
        A = Z + 1 + alpha - int( alpha / 4 )

    B = A + 1524
    C = int( ( B - 122.1 ) / 365.25 )
    D = int( 365.25 * C )
    E = int( ( B - D ) / 30.6001 )

    dayOfMonth = B - D - int( 30.6001  * E ) + F

    if E < 14:
        month = E - 1

    else:
        month = E - 13

    if month > 2:
        year = C - 4716

    else:
        year = C - 4715

    return datetime.datetime( year, month, 1 ) + datetime.timedelta( days = dayOfMonth - 1 )


# This function is unfinished...and likely very wrong!
# Cannot make sense of the way to determine eclipse type.
#
# Sample implementations...not sure which is correct, if any.
# https://github.com/pavolgaj/AstroAlgorithms4Python/blob/master/eclipse.py
# https://github.com/soniakeys/meeus/blob/master/v3/eclipse/eclipse.go
# https://github.com/soniakeys/meeus/blob/master/v3/eclipse/eclipse_test.go
def getType( gamma, u ):
    if math.fabs( gamma ) < 0.9972:
        eclipseType = EclipseType.CENTRAL
        if u < 0:
            eclipseType = EclipseType.TOTAL

        elif u > 0.0047:
            eclipseType = EclipseType.ANNULAR

        else:
            if u < 0.00464 * math.sqrt( 1 - gamma * gamma ):
                eclipseType = EclipseType.ANNULAR_TOTAL

            else:
                eclipseType = EclipseType.ANNULAR

    elif math.fabs( gamma ) > 1.5433 + u:
        eclipseType = EclipseType.NONE

    else: # Not central...partial.
        eclipseType = EclipseType.NOT_CENTRAL
        if math.fabs( gamma ) > 0.9972 and math.fabs( gamma ) < 1.0260:
            pass # Not sure what to do here.

        else:
            eclipseType = EclipseType.PARTIAL

    return eclipseType


# Chapter 54
# Astronomical Algorithms, Second Edition
# Jean Meeus
def getEclipse( utcNow, isSolar, upcomingMoonDates ):
    eclipseDate = None
    date = utcNow
    moonIndex = 0
    while moonIndex <= len( upcomingMoonDates ):
        daysInYear = 365
        if calendar.isleap( date.year ):
            daysInYear = 366

        yearAsDecimal = date.year + ( date.timetuple().tm_yday / daysInYear )
        k = ( yearAsDecimal - 2000 ) * 12.3685 # Equation 49.2
        k = round( k ) # Must be an integer for a New Moon (solar eclipse).
        if not isSolar:
            k += 0.5 # Must be an integer increased by 0.5 for a Full Moon (lunar eclipse).

        T =  k / 1236.85 # Equation 49.3

        JDE =   2451550.09766 \
              + 29.530588861 * k \
              + 0.00015437 * T * T \
              - 0.000000150 * T * T * T \
              + 0.00000000073 * T * T * T * T # Equation 49.1

        E = 1 - 0.002516 * T - 0.0000074 * T * T # Equation 47.6

        M = (  2.5534 \
             + 29.10536570 * k \
             - 0.0000014 * T * T \
             - 0.00000011 * T * T * T ) % 360 # Equation 49.4

        Mdash = (  201.5634 \
                 + 385.81693528 * k \
                 + 0.0107582 * T * T \
                 + 0.00001238 * T * T * T \
                 - 0.000000058 * T * T * T * T ) % 360 # Equation 49.5

        F = (  160.7108 \
             + 390.67050284 * k \
             - 0.0016118 * T * T \
             - 0.000000227 * T * T * T \
             + 0.000000011 * T * T * T * T ) % 360 # Equation 49.6

        noEclipse = math.fabs( math.sin( math.radians( F ) ) ) > 0.36
        if noEclipse:
            date = upcomingMoonDates[ moonIndex ] #TODO Maybe add one day?
            moonIndex += 1
            continue

        omega = (  124.7746 \
                 - 1.56375588 * k \
                 + 0.0020672 * T * T \
                 + 0.00000215 * T * T * T ) % 360 # Equation 49.7

        F1 = F - 0.02665 * math.sin( math.radians( omega ) )

        A1 = 299.77 + 0.107408 * k - 0.009173 * T * T

        # Equation 54.1
        if isSolar:
            correctedJDE = JDE \
                           - 0.4075 * math.sin( math.radians( Mdash ) ) \
                           + 0.1721 * E * math.sin( math.radians( M ) )

        else:
            correctedJDE = JDE \
                           - 0.4065 * math.sin( math.radians( Mdash ) ) \
                           + 0.1727 * E * math.sin( math.radians( M ) )

        correctedJDE +=   0.0161 * math.sin( math.radians( 2 * Mdash ) ) \
                        - 0.0097 * math.sin( math.radians( 2 * F1 ) ) \
                        + 0.0073 * E * math.sin( math.radians( Mdash - M ) ) \
                        - 0.0050 * E * math.sin( math.radians( Mdash + M ) ) \
                        - 0.0023 * math.sin( math.radians( Mdash - 2 * F1 ) ) \
                        + 0.0021 * E * math.sin( math.radians( 2 * M ) ) \
                        + 0.0012 * math.sin( math.radians( Mdash + 2 * F1 ) ) \
                        + 0.0006 * E * math.sin( math.radians( 2 * Mdash + M ) ) \
                        - 0.0004 * math.sin( math.radians( 3 * Mdash ) ) \
                        - 0.0003 * E * math.sin( math.radians( M + 2 * F1 ) ) \
                        + 0.0003 * math.sin( math.radians( A1 ) ) \
                        - 0.0002 * E * math.sin( math.radians( M - 2 * F1 ) ) \
                        - 0.0002 * E * math.sin( math.radians( 2 * Mdash - M ) ) \
                        - 0.0002 * math.sin( math.radians( omega ) )

        eclipseDate = julianDayToGregorian( correctedJDE )

        P =   0.2070 * E * math.sin( math.radians( M ) ) \
            + 0.0024 * E * math.sin( math.radians( 2 * M ) ) \
            - 0.0392 * math.sin( math.radians( Mdash ) ) \
            + 0.0116 * math.sin( math.radians( 2 * Mdash ) ) \
            - 0.0073 * E * math.sin( math.radians( Mdash + M ) ) \
            + 0.0067 * E * math.sin( math.radians( Mdash - M ) ) \
            + 0.0118 * math.sin( math.radians( 2 * F1 ) )

        Q =   5.2207 \
            - 0.0048 * E * math.cos( math.radians( M ) ) \
            + 0.0020 * E * math.cos( math.radians( 2 * M ) ) \
            - 0.3299 * math.cos( math.radians( Mdash ) ) \
            - 0.0060 * E * math.cos( math.radians( Mdash + M ) ) \
            + 0.0041 * E * math.cos( math.radians( Mdash - M ) )

        W = math.fabs( math.cos( math.radians( F1 ) ) )

        gamma = ( P * math.cos( math.radians( F1 ) ) + Q * math.sin( math.radians( F1 ) ) ) * ( 1 - 0.0048 * W )
        
        u =   0.0059 \
            + 0.0046 * E * math.cos( math.radians( M ) ) \
            - 0.0182 * math.cos( math.radians( Mdash ) ) \
            + 0.0004 * math.cos( math.radians( 2 * Mdash ) ) \
            - 0.0005 * math.cos( math.radians( M + Mdash ) )

        eclipseType = getType( gamma, u )

        break

    return eclipseDate, eclipseType


# Verify
#    lunar: https://eclipse.gsfc.nasa.gov/5MCLE/5MKLEcatalog.txt
#    solar: https://eclipse.gsfc.nasa.gov/5MCSE/5MKSEcatalog.txt
now = datetime.datetime.utcnow()
timeScale = load.timescale( builtin = True )
ephemeris = load( "de421.bsp" )
fullMoonDates, newMoonDates = getMoonDatesForOneYear( now, timeScale, ephemeris )
print( "Meeus solar eclipse:", getEclipse( now, True, newMoonDates ) )
print( "Meeus lunar eclipse:", getEclipse( now, False, fullMoonDates ) )

[brandon-rhodes]

@Bernmeister — Thanks for sharing your work, I’ll read this over!

Did you happen to see the commit listed just above your comment? There is now a lunar eclipse routine in Skyfield. You’re welcome to try it out before the next release with:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

Here are a few thoughts, looking at your code for the first time:

1. It had not occurred to me to use the existing moon phases routine to learn the dates and times to check for eclipses!

2. If you would like to use it (though of course sometimes it’s nicer to be in control of how a calculation is performed), you might find that skyfield.timelib.compute_calendar_date() would return the same numbers you are computing in julianDayToGregorian(). (The more official interface in Skyfield, I suppose, is something more like ymdhms = ts.tt_jd(my_jd).tt_calendar(), but that produces the time too.)

3. Is the Meeus routine doing its own calculation of where the Moon and Sun are in the sky? I don’t see Skyfield loading an ephemeris. My guess is that this would create a confusing mismatch for folks using a modern high-accuracy ephemeris, since the exact positions of the Sun and Moon that Skyfield gives them would not quite match the numbers here — do you know if this is using VSOP87, or another theory of the Moon and Sun’s positions? I don’t recognize the numbers offhand.

4. Do you have a copy of the Meeus book? In recently writing the lunar eclipse support I mentioned above, I used a book, and I found its diagrams essential in understanding what the different quantities were. You might find that the book your code comes from would also itself have helpful diagrams, that might clear up the tests the code is doing to classify the different kinds of eclipse.

I am planning on releasing Skyfield in the next couple of days, with the new lunar eclipse logic you can see in the commit above, but without solar eclipse logic yet — that’s in the next section of the book I have, and I haven’t had time to write it up yet.

Did you write the code above yourself from a book, or were you able to cut and paste it from the GitHub projects you link to?

[Bernmeister]

@brandon-rhodes

Did you happen to see the commit listed just above your comment?

Only noticed it after you've pointed it out...oh well.

3. Is the Meeus routine doing its own calculation of where the Moon and Sun are in the sky?

No idea how the calculation works under the covers. Did my best to "understand" the underlying algorithm, but ultimately, I just implemented things as best as I could and verified the final answer of eclipse date/time rather than verifying the algorithm itself.

When it came to the eclipse type (from Meeus) I found the description (differences between solar and lunar eclipse types) ambiguous and confusing at best (fault is on my part, not Meeus).

Did you write the code above yourself from a book, or were you able to cut and paste it from the GitHub projects you link to?

All code written by me but based exclusively on Meeus (the Duffett-Smith/Zwart book was used for the Julian stuff). I tried to find existing code to verify results and/or help me understand the algorithm but in the end, I wrote code to implement the formulae in the book(s).

[brandon-rhodes]

Only noticed it after you've pointed it out...oh well.

If you were following the issue but GitHub didn‘t light up your notifications, I'll be happy in the future to make a comment after a commit to try to draw attention to it.

… in the end, I wrote code to implement the formulae in the book(s).

My approach with the commit above was similar: reading, followed by just converting their formulae into Python, with a focus on getting the right answer rather than understanding immediately. 🙂

Only once it worked did I review the text, figure out what the symbols meant, and substituted in some cases meaningful names for what had been plain variables in the original formulae.

Hopefully your code helps other folks who need eclipse data before solar eclipses are implemented in Skyfield, which probably won’t be for a few weeks, even if I put it at the front of my to-do list! My guess is that it would be hard to incorporate into official Skyfield code, because I think many of those equations are approximations for things Skyfield does at high precision, like TDB time, and the Sun and Moon positions — but to untangle which are which, we would need to hunt down the source of all of the coefficients in your example code. Of which, wow, there are a lot! You did a lot of work to successfully bring them into Python, and I'm glad for the moment that they provide a good match for at least the times of both lunar and solar eclipses.

[lkangas]

In preparation of the 2024-04-08 eclipse, I've written some code utilizing Skyfield for calculating the Besselian elements for general and local eclipse circumstances as described in the well known Explanatory supplement by Urban and Seidelmann. For example accurate calculation of contact times for a fixed location works well. Using Skyfield, it's easy and fast to accurately calculate the actual elements as functions of time, no need for the usual polynomial tables.

After this years eclipse, I expect to be diving into even more calculations like shadow path calculations and also more general searching of eclipse occurrences.

So just asking, what is the status of the solar eclipse calculations, and would the formulae I've implemented be useful here?

[brandon-rhodes]

@lkangas — I would be interested in seeing your approach. All that I myself have had time to do can be seen in this file in the repository:

https://github.com/skyfielders/python-skyfield/blob/master/design/solar_eclipse.py

As you can see from the resulting image:

— all that I've managed to do so far is the bare-bones step of drawing the path of the shadow across the surface. As yet there is no work on things like size of the shadow, length of totality, or percent-sun-covered for areas outside of totality.

I haven't yet decided which code to write first: whether to work further on solar eclipses right now, or whether to take the routine I just wrote line_and_ellipsoid_intersection() and get it working more generally, so that people can use it with other planets besides Earth, and for other situations besides eclipses. Like: someone was asking about the concept of "geocentric altitude" over in another issue, which I think might be solvable with the right call to line_and_ellipsoid_intersection(). So material from you showing further eclipse calculations could get those available to Skyfield users much sooner than if we wait around for me to eventually work on them!

[lkangas]

Roger, once the eclipse is over, I will put together some proof of concept and share it!

Line-ellipsoid intersection is surely useful in general. For solar eclipses specifically, using Bessel's approach with a fundamental plane perpendicular to the shadow axis makes everything super easy. The intersection is found simply by applying two rotation matrices, and converting geodetic latitude into geocentric latitude.

The local circumstances (contact times, umbral duration, percent coverage, etc) are also best carried out using the Besselian elements. I found that Skyfields find_discrete and find_minima work excellent for these purposes.

I would also suppose that the generation of Besselian elements themselves would be interesting for almanac applications.