<p>Per <a href="https://astronomy.stackexchange.com/posts/comments/118753?noredirect=1">@PM 2Ring's comment</a>, I needed to use apparent RA/declination rather than GCRS (which is the default returned by <code>astropy.coordinates.get_body</code>). In astropy, this is implemented as the <a href="https://docs.astropy.org/en/stable/api/astropy.coordinates.TETE.html#astropy.coordinates.TETE" rel="noreferrer">True Equator True Equinox</a> coordinate system.</p>
<p>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.</p>
<p>Incorporating both of those helpful comments leads to this script:</p>
<pre><code>import scipy.optimize
import astropy.units as u
import astropy.time
import astropy.coordinates

astropy.coordinates.solar_system_ephemeris.set(&quot;de441_part-2&quot;)

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

print(time_solstice.iso)
</code></pre>
<p>which prints:</p>
<pre><code>2023-06-21 14:57:49.087
</code></pre>
<p>This rounds to the same value as wikipedia.</p>

CC BY-SA 4.0
2023-03-26T23:06:34.277Z
2023-03-27T06:08:15.793Z
2023-03-27T06:08:15.793Z
Roy Smart
<p>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 <code>2023-06-21 23:00 UTC</code> instead of <code>2023-06-21 14:58 UTC</code> according to <a href="https://en.wikipedia.org/wiki/Summer_solstice" rel="noreferrer">wikipedia</a>. Does anyone know where I'm going wrong?</p>
<pre><code>import matplotlib.pyplot as plt
import astropy.time
import astropy.coordinates
import astropy.visualization
import numpy as np

astropy.coordinates.solar_system_ephemeris.set(&quot;de440&quot;)

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

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

sun = astropy.coordinates.get_body(&quot;sun&quot;, time=time)

plt.figure()
plt.plot(time, sun.dec)
plt.show()
</code></pre>
<p><a href="https://i.sstatic.net/hcmtU.png" rel="noreferrer"><img src="https://i.sstatic.net/hcmtU.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-03-26T22:26:53.290Z
Horizons agrees with Wikipedia https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=YES&COMMAND=10&EPHEM_TYPE=OBSERVER&CENTER=500%40399&START_TIME=%272023-Jun-21+12%3A00%27&STOP_TIME=%272023-Jun-21+18%3A00%27&STEP_SIZE=1h&QUANTITIES=%271%2C2%2C31%27&ANG_FORMAT=DEG&APPARENT=AIRLESS&OBJ_DATA=NO You need the apparent RA & declination, rather than the astrometric.
PM 2Ring
CC BY-SA 4.0
2023-03-27T00:37:53.227Z
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.
Greg Miller
CC BY-SA 4.0
2023-03-27T17:10:20.613Z
Man it's tempting to link this question into flat earth circles
Cruncher
<observational-astronomy><the-sun><astropy><solstice>
How do I calculate the moment of the solstice?
Roy Smart
# How do I calculate the moment of the solstice?

(Asked by: username_0 on 2023-03-26)

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

**Tags:** observational-astronomy, the-sun, astropy, solstice

> **Comments:**
> **username_1**: Horizons agrees with Wikipedia https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=YES&COMMAND=10&EPHEM_TYPE=OBSERVER&CENTER=500%40399&START_TIME=%272023-Jun-21+12%3A00%27&STOP_TIME=%272023-Jun-21+18%3A00%27&STEP_SIZE=1h&QUANTITIES=%271%2C2%2C31%27&ANG_FORMAT=DEG&APPARENT=AIRLESS&OBJ_DATA=NO You need the apparent RA & declination, rather than the astrometric.  
> **username_2**: 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.  
> **username_3**: Man it's tempting to link this question into flat earth circles  

---

## 1 answers

### Answer by username_0 on 2023-03-26. Score: 7

Per @username_1'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.


---

