How do I compute equinox/solstice moments?

Question

What algorithms or formulas are available for computing the equinoxes and solstices? I found one of these a few years ago and implemented it, but the precision was not great: the time of day seemed to be assumed at 00:00, 06:00, 12:00, and 18:00 UTC depending on which equinox or solstice was computed. Wikipedia gives these computed out to the minute, so something more exact must be possible. Libraries for my favorite programming language also come out to those hardcoded times, so I assume they are using the same or a similar algorithm as the one I implemented.

I also once tried using a library that gave me the solar longitude and implementing a search routine to zero in on the exact moments of 0, 90, 180, and 270 degrees; this worked down to the second but did not agree with the times in Wikipedia, so I assume there was something wrong with this approach. I am, however, pleasantly surprised to discover that Maimonides (medieval Jewish scholar) proposed an algorithm using the exact same idea a millenium ago.

Answer 1

A great source for the (complex!) underlying formulas and algorithms is Astronomical Algorithms by Jean Meeus.

Using the PyMeeus implementation of those algorithms, and the code below, you can get the following values for the 2018 winter solstice (where "winter" refers to the northern hemisphere).

winter solstice for 2018 in Terrestrial Time is at:
 (2018, 12, 21, 22, 23, 52.493725419044495)

winter solstice for 2018 in UTC, if last leap second was (2016, 12):
 (2018, 12, 21, 22, 22, 43.30972542127711)

winter solstice for 2018 in local time, if last leap second was (2016, 12)
 and local time offset is -7.00 hours:
 (2018, 12, 21, 15, 22, 43.30973883232218)

i.e. 2018-12-21T15:22:43.309725-07:00

Of course, the answer is not accurate down to microseconds, but I also wanted to show how to do high-precision conversions with arrow.

Code:

from pymeeus.Sun import Sun
from pymeeus.Epoch import Epoch

year = 2018  # datetime.datetime.now().year
target="winter"

# Get terrestrial time of given solstice for given year
solstice_epoch = Sun.get_equinox_solstice(year, target=target)

print("%s solstice for %d in Terrestrial Time is at:\n %s" %
      (target, year, solstice_epoch.get_full_date()))

print("%s solstice for %d in UTC, if last leap second was %s:\n %s" %
 (target, year, Epoch.get_last_leap_second()[:2], solstice_epoch.get_full_date(utc=True)))

solstice_local = (solstice_epoch + Epoch.utc2local()/(24*60*60))
print("%s solstice for %d in local time, if last leap second was %s\n"
 " and local time offset is %.2f hours:\n %s" %
 (target, year, Epoch.get_last_leap_second()[:2],
  Epoch.utc2local() / 3600., solstice_local.get_full_date(utc=True)))

Using the very cool more ISO and TZ aware module Arrow: better dates and times for Python, that can be printed more nicely:

import arrow
import math

slutc = solstice_epoch.get_full_date(utc=True)
frac, whole = math.modf(slutc[5])

print("i.e. %s" % arrow.get(*slutc[:5], int(whole), round(frac * 1e6)).to('local'))

Answer 2

I'm not sure if this is an accurate enough solution for you, but I found a NASA website that has some code snippets for calculating the vernal equinox as well as some other astronomical-type information. I've also found some references to a book called Astronomical Algorithms which may have the answers you need if the info somehow isn't available online.

Answer 3

I know you're looking for something that'll paste into an answer here, but I have to mention SPICE, a toolkit produced by NAIF at JPL, funded by NASA. It might be overkill for Farmer's Almanac stuff, but you mentioned interest in precision and this toolkit is routinely used in planetary science.