8
$\begingroup$

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()

enter image description here

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Horizons agrees with Wikipedia ssd.jpl.nasa.gov/api/… You need the apparent RA & declination, rather than the astrometric. $\endgroup$
    – PM 2Ring
    Commented Mar 26, 2023 at 22:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Greg Miller
    Commented Mar 27, 2023 at 0:37
  • $\begingroup$ Man it's tempting to link this question into flat earth circles $\endgroup$
    – Cruncher
    Commented Mar 27, 2023 at 17:10

1 Answer 1

Reset to default
7
$\begingroup$

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.

Incorporating both of those helpful comments leads to this script:

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.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .