# How do I calculate the moment of the solstice?

I'm trying to calculate the moment of the upcoming summer solstice using astropy, and I can't seem to get an answer that matches the accepted values that I see online. I thought that the moment of the solstice would be when the Sun has maximum/minimum declination, but according to my plot below that moment is at about 2023-06-21 23:00 UTC instead of 2023-06-21 14:58 UTC according to wikipedia. Does anyone know where I'm going wrong?

import matplotlib.pyplot as plt
import astropy.time
import astropy.coordinates
import astropy.visualization
import numpy as np

astropy.coordinates.solar_system_ephemeris.set("de440")

astropy.visualization.quantity_support()
astropy.visualization.time_support()

time = np.linspace(astropy.time.Time("2023-06-21"), astropy.time.Time("2023-06-22"), num=1000)

sun = astropy.coordinates.get_body("sun", time=time)

plt.figure()
plt.plot(time, sun.dec)
plt.show()


• Horizons agrees with Wikipedia ssd.jpl.nasa.gov/api/… You need the apparent RA & declination, rather than the astrometric. Mar 26, 2023 at 22:26
• Also, the instant of the solstice is when the Sun's ecliptic longitude is 90 degrees. It will be slightly different than when the maximun declination is reached. Mar 27, 2023 at 0:37
• Man it's tempting to link this question into flat earth circles Mar 27, 2023 at 17:10

Per @PM 2Ring's comment, I needed to use apparent RA/declination rather than GCRS (which is the default returned by astropy.coordinates.get_body). In astropy, this is implemented as the True Equator True Equinox coordinate system.

Also, @GregMiller points out that the instant of the summer solstice is the moment that the Sun's ecliptic longitude is 90 degrees, not when declination reaches a maximum.

import scipy.optimize
import astropy.units as u
import astropy.time
import astropy.coordinates

astropy.coordinates.solar_system_ephemeris.set("de441_part-2")

time_solstice = astropy.time.Time(
scipy.optimize.root_scalar(
f=lambda t: (astropy.coordinates.get_body("sun", time=astropy.time.Time(t, format="unix")).tete.ra - 90 * u.deg).value,
bracket=[
astropy.time.Time("2023-06-20").unix,
astropy.time.Time("2023-06-23").unix,
],
).root,
format="unix",
)

print(time_solstice.iso)


which prints:

2023-06-21 14:57:49.087


This rounds to the same value as wikipedia.