find_minima for sun declination doesn't match almanac seasons

[ordemann]

The times of the seasons using find_minima/maxima for the sun declination don't match what comes from the skyfield almanac.

I used Skyfield almanac to find the earth seasons and find_minima/maxima to find the earth apsides. The results match what's on the US Naval Observatory site for 2022 [https://aa.usno.navy.mil/calculated/seasons?year=2022&tz=0&tz_sign=-1&tz_label=false&dst=false&submit=Get+Data].

I looked up the sun declination for the time skyfield almanac gives for equinoxes and was surprised that the angles weren't exactly 0 degrees. I proceeded to find the equinoxes and solstices using find_minima/maxima of the sun declination via radec() and found a discrepancy.

I must be doing something wrong. I've read through the skyfield documentation and tried a bunch of things, but I can't get the two approaches to finding the seasons to match. Any help is appreciated.

Thanks!

Here's my code:
`
"""
Given a year, find the earth seasons and apsides.
Check the seasons against min/max sun declination
"""

import sys
import skyfield
from skyfield.api import load
from skyfield import almanac
from skyfield.searchlib import find_minima, find_maxima

print('Python Version:', sys.version)
print('Skyfield Version:', skyfield.VERSION)
print('\n\n')

ephemeris = load('de440s.bsp')

def earth_distance(time):
earth = ephemeris['earth']
sun = ephemeris['sun']
e = earth.at(time)
s = e.observe(sun)
return s.distance().km
earth_distance.step_days = 1

def sun_declination(time):
earth = ephemeris['earth']
sun = ephemeris['sun']
e = earth.at(time)
s = e.observe(sun)
ra, declin, dist = s.radec(time)
return abs(declin.degrees)
sun_declination.step_days = 1

def seasons(year):
"""
Find the earth seasons and apsides.
Use the Skyfield almanac and check against the
min/max sun declinations which define the seasons
"""
### Set up Skyview and initialize data ###
event_time = 0
dist = angle = 1
timescale = load.timescale()
start_time = timescale.utc(year, 1, 1)
end_time = timescale.utc(year+1, 1, 1)

# Find seasons and corresponding the Sun declination angle
# Note: the declination is displayed as an absolute value
print('#### SEASONS ####')
eq_times, eq_codes = almanac.find_discrete(
    start_time, end_time, almanac.seasons(ephemeris))
for eq_code, eq_time in zip(eq_codes, eq_times):
    print(f'{almanac.SEASON_EVENTS[eq_code]:18} ' +
          f'{eq_time.utc_datetime()}  ' +
          f'{sun_declination(eq_time):3.5f}°')

# Find the apsides
print('\n#### APSIDES ####')
aphelion = find_maxima(start_time, end_time, earth_distance)
perihelion = find_minima(start_time, end_time, earth_distance)
print(f'{"Aphelion":12} ' +
      f'{aphelion[event_time][0].utc_datetime()}  ' +
      f'{aphelion[dist][0]:.0f} km')
print(f'{"Perihelion":12} ' +
      f'{perihelion[event_time][0].utc_datetime()}  ' +
      f'{perihelion[dist][0]:.0f} km')

# Find dates corresponding to Min/Max Sun declination.
# The dates and angles should match those of the seasons.
print('\n#### MIN/MAX DECLINATIONS ####')
min_decl = find_minima(start_time, end_time, sun_declination)
max_decl = find_maxima(start_time, end_time, sun_declination)
print(f'{min_decl[event_time][0].utc_datetime()} ' +
      f'{min_decl[angle][0]:3.5f}°')
print(f'{max_decl[event_time][0].utc_datetime()} ' +
      f'{max_decl[angle][0]:3.5f}°')
print(f'{min_decl[event_time][1].utc_datetime()} ' +
      f'{min_decl[angle][1]:3.5f}°')
print(f'{max_decl[event_time][1].utc_datetime()} ' +
      f'{max_decl[angle][1]:3.5f}°')

`
And, here's the output:

`
Python Version: 3.9.12 (main, May 8 2022, 17:57:49)
[Clang 13.1.6 (clang-1316.0.21.2)]
Skyfield Version: (1, 42)

SEASONS

Vernal Equinox 2022-03-20 15:33:24.902557+00:00 0.00228°
Summer Solstice 2022-06-21 09:13:51.057466+00:00 23.43780°
Autumnal Equinox 2022-09-23 01:03:42.256731+00:00 0.00206°
Winter Solstice 2022-12-21 21:48:13.006813+00:00 23.43814°

APSIDES

Aphelion 2022-07-04 07:10:39.145961+00:00 152098463 km
Perihelion 2022-01-04 06:54:25.954538+00:00 147105046 km

MIN/MAX DECLINATIONS

2022-03-20 15:25:06.501680+00:00 0.00000°
2022-06-21 09:10:39.667275+00:00 23.43780°
2022-09-23 00:56:05.255710+00:00 0.00000°
2022-12-21 21:44:23.066157+00:00 23.43814°

`

[brandon-rhodes]

Try editing your answer to triple-backtick your code; it should fix the formatting problems, which make it a little difficult to wade into right now. But at first glance: it looks like you're determining declination from the astrometric position, whereas I think the Naval Observatory uses the apparent position? Try s = e.observe(sun).apparent() instead. The differences are discussed here:

https://rhodesmill.org/skyfield/positions.html#barycentric-astrometric-apparent

Let us know if that makes a difference in your results.

[ordemann]

Using 'apparent()' gets the equinox times within 60 seconds of the almanac time, but the solstice times are still 4+ minutes off.

Vernal Equinox     2022-03-20 15:33:24.902557+00:00  0.00001°
Summer Solstice    2022-06-21 09:13:51.057466+00:00  23.43780°
Autumnal Equinox   2022-09-23 01:03:42.256731+00:00  0.00020°
Winter Solstice    2022-12-21 21:48:13.006813+00:00  23.43814°

#### APSIDES ####
Aphelion     2022-07-04 07:10:39.145961+00:00  152098463 km
Perihelion   2022-01-04 06:54:25.954538+00:00  147105046 km

#### MIN/MAX DECLINATIONS ####
2022-03-20 15:33:23.765247+00:00  0.00000°
2022-06-21 09:19:06.534652+00:00  23.43780°
2022-09-23 01:04:25.720547+00:00  0.00000°
2022-12-21 21:52:34.460796+00:00  23.43814°

Here's the original result, i.e. not using 'apparent()'

Python Version: 3.9.12 (main, May  8 2022, 17:57:49) 
[Clang 13.1.6 (clang-1316.0.21.2)]
Skyfield Version: (1, 42)



#### SEASONS ####
Vernal Equinox     2022-03-20 15:33:24.902557+00:00  0.00001°
Summer Solstice    2022-06-21 09:13:51.057466+00:00  23.43780°
Autumnal Equinox   2022-09-23 01:03:42.256731+00:00  0.00020°
Winter Solstice    2022-12-21 21:48:13.006813+00:00  23.43814°

#### APSIDES ####
Aphelion     2022-07-04 07:10:39.145961+00:00  152098463 km
Perihelion   2022-01-04 06:54:25.954538+00:00  147105046 km

#### MIN/MAX DECLINATIONS ####
2022-03-20 15:33:23.765247+00:00  0.00000°
2022-06-21 09:19:06.534652+00:00  23.43780°
2022-09-23 01:04:25.720547+00:00  0.00000°
2022-12-21 21:52:34.460796+00:00  23.43814°

And, here's the code again with proper formatting:

"""
Given a year, find the earth seasons and apsides.
Check the seasons against min/max sun declination
"""

import sys
import skyfield
from skyfield.api import load
from skyfield import almanac
from skyfield.searchlib import find_minima, find_maxima


print('Python Version:', sys.version)
print('Skyfield Version:', skyfield.VERSION)
print('\n\n')


ephemeris = load('de440s.bsp')

def earth_distance(time):
    earth = ephemeris['earth']
    sun = ephemeris['sun']
    e = earth.at(time)
    s = e.observe(sun)
    return s.distance().km
earth_distance.step_days = 1


def sun_declination(time):
    earth = ephemeris['earth']
    sun = ephemeris['sun']
    e = earth.at(time)
    s = e.observe(sun).apparent()   # only line changed between two results
    ra, declin, dist = s.radec(time)
    return abs(declin.degrees)
sun_declination.step_days = 1


def seasons(year):
    """
    Find the earth seasons and apsides.
    Use the Skyfield almanac and check against the
    min/max sun declinations which define the seasons
    """
    ### Set up Skyview and initialize data ###
    event_time = 0
    dist = angle = 1
    timescale = load.timescale()
    start_time = timescale.utc(year, 1, 1)
    end_time = timescale.utc(year+1, 1, 1)

    # Find seasons and corresponding the Sun declination angle
    # Note: the declination is displayed as an absolute value
    print('#### SEASONS ####')
    eq_times, eq_codes = almanac.find_discrete(
        start_time, end_time, almanac.seasons(ephemeris))
    for eq_code, eq_time in zip(eq_codes, eq_times):
        print(f'{almanac.SEASON_EVENTS[eq_code]:18} ' +
              f'{eq_time.utc_datetime()}  ' +
              f'{sun_declination(eq_time):3.5f}°')

    # Find the apsides
    print('\n#### APSIDES ####')
    aphelion = find_maxima(start_time, end_time, earth_distance)
    perihelion = find_minima(start_time, end_time, earth_distance)
    print(f'{"Aphelion":12} ' +
          f'{aphelion[event_time][0].utc_datetime()}  ' +
          f'{aphelion[dist][0]:.0f} km')
    print(f'{"Perihelion":12} ' +
          f'{perihelion[event_time][0].utc_datetime()}  ' +
          f'{perihelion[dist][0]:.0f} km')

    # Find dates corresponding to Min/Max Sun declination.
    # The dates and angles should match those of the seasons.
    print('\n#### MIN/MAX DECLINATIONS ####')
    min_decl = find_minima(start_time, end_time, sun_declination)
    max_decl = find_maxima(start_time, end_time, sun_declination)
    print(f'{min_decl[event_time][0].utc_datetime()}  ' +
          f'{min_decl[angle][0]:3.5f}°')
    print(f'{max_decl[event_time][0].utc_datetime()}  ' +
          f'{max_decl[angle][0]:3.5f}°')
    print(f'{min_decl[event_time][1].utc_datetime()}  ' +
          f'{min_decl[angle][1]:3.5f}°')
    print(f'{max_decl[event_time][1].utc_datetime()}  ' +
          f'{max_decl[angle][1]:3.5f}°')

seasons(2022)

[brandon-rhodes]

Thanks for sharing both your results and your scripts—and, I can read the formatting very easily this time, thanks!

Is it possible that the moment of the solstice simply doesn't coincide with the moment of maximum or minimum declination? You might try looking up the definitions to be sure, but I think that the solstice is the moment at which the ecliptic longitude reaches a certain angle (maybe 90° and 270°?), regardless of whether the declination happens to be at a maximum or minimum at that moment.

[gyth]

I do not know if it is useful: the definition of vernal equinox as given by IMCCE (French Institute of Celestial Mechanics and Calculation of Ephemerides) is "the apparent geocentric longitude of the center of the Sun is zero degrees".
sorry article in french but it have a interesting drawing interestingly enough they give 15 h 33 min 26,14 s UTC for the equinox 20th march 2022. Most probably based on INPOP planetary planetary ephemeris : not de440s.bsp

[brandon-rhodes]

Yes, thanks, that's what I was remembering: it's based on the Sun's longitude in the sky. If the vernal equinox is at 90° ecliptic longitude, then the solstices must indeed be at 90° and 270°.

[ordemann]

Fantastic, thanks Brandon and Gyth for your quick replies. Indeed, it seems I have an incorrect definition of solstice burned into my brain. Also, I'll give the INPOP ephemeris a try.