find syzygy / equinox / solstice with astropy

syzygy_astropy.py

# Find the nearest syzygy (solstice or equinox) to a given date.
# Solution for https://stackoverflow.com/questions/55838712/finding-equinox-and-solstice-times-with-astropy
# TODO: need to ensure we're using the right sun position functions for the syzygy definition....

import math
from astropy.time import Time, TimeDelta
import astropy.coordinates
from scipy.optimize import brentq
from astropy import units as u


# We'll usually find a zero crossing if we look this many days before and after
# given time, except when it is is within a few days of a cross-quarter day.
# But going over 50,000 years back, season lengths can vary from 85 to 98 days!
#  https://individual.utoronto.ca/kalendis/seasons.htm#seasons
delta = 44.0


def mjd_to_time(mjd):
    "Return a Time object corresponding to the given Modified Julian Date."

    return Time(mjd, format='mjd', scale='utc')


def sunEclipticLongitude(mjd):
    "Return ecliptic longitude of the sun in degrees at given time (MJD)"

    t = mjd_to_time(mjd)
    sun = astropy.coordinates.get_body('sun', t)
    # TODO: Are these the right functions to call for the definition of syzygy? Geocentric? True?
    eclipticOfDate = astropy.coordinates.GeocentricTrueEcliptic(equinox=t)
    sunEcliptic = sun.transform_to(eclipticOfDate)
    return sunEcliptic.lon.deg


def linearize(angle):
    """Map angle values in degrees near the quadrants of the circle
    into smooth linear functions crossing zero, for root-finding algorithms.
    Note that for angles near 90 or 270, increasing angles yield decreasing results

    >>> linearize(5) > 0 > linearize(355)
    True
    >>> linearize(95) > 0 > linearize(85)
    False
    """

    return math.sin(math.radians(angle * 2))


def map_syzygy(t):
    "Map times into linear functions crossing zero at each syzygy"

    return linearize(sunEclipticLongitude(t))


def find_nearest_syzygy(t):
    """Return the precise Time of the nearest syzygy to the given Time,
    which must be within 43 days of one syzygy"
    """

    syzygy_mjd = brentq(map_syzygy, t.mjd - delta, t.mjd + delta)

    syzygy = mjd_to_time(syzygy_mjd)
    syzygy.format = 'isot'
    return syzygy


if __name__ == '__main__':
    import doctest
    doctest.testmod()

    t0 = Time('2019-09-23T07:50:10', format='isot', scale='utc')

    td = TimeDelta(1.0 * u.day)

    seq = t0 + td * range(0, 365, 15)

    for t in seq:
        try:
            syzygy = find_nearest_syzygy(t)
        except ValueError as e:
            print(f'{e=}, {t.value=}, {t.mjd-delta=}, {map_syzygy(t.mjd-delta)=}')
            continue
        print(f'{t.value=}, {syzygy.value=}, {sunEclipticLongitude(syzygy)=}')