I need to compute the date and time of the vernal equinox for a given year. I have the following code:
from astropy.time import Time
from astropy.coordinates import get_sun, AltAz, EarthLocation
import numpy as np
import astropy.units as u
def find_vernal_equinox(year):
# Make an array of times around the expected equinox (March 19-21)
times = Time(f"{year}-03-19 00:00:00") + np.linspace(0, 3, 1000) * u.day
# Get the declination of the Sun at those times
sun_coords = get_sun(times)
declinations = sun_coords.dec.deg
# The equinox occurs when declination crosses 0
zero_crossings = np.where(np.diff(np.sign(declinations)))[0]
if not len(zero_crossings):
raise ValueError("No zero crossing found!")
idx = zero_crossings[0]
# Interpolate to get a more precise time
t1, t2 = times[idx], times[idx+1]
d1, d2 = declinations[idx], declinations[idx+1]
# Linear interpolation
frac = -d1 / (d2 - d1)
equinox_time = t1 + frac * (t2 - t1)
return equinox_time.utc.iso
# Example usage:
print(find_vernal_equinox(2025))
Strangely, the results I get are about 8 or 9 hours off (all times UTC, actual data from https://data.giss.nasa.gov/modelE/ar5plots/srvernal.html):
Computed 2025-03-20 17:34:03.719 – Actual 2025 3/20 9:00
Computed 2026-03-20 23:39:44.858 – Actual 2026 3/20 14:49
Computed 2027-03-21 05:41:50.943 – Actual 2027 3/20 20:38
Moreover, the declination of the Sun at the time of the equinox is not zero:
print(get_sun(Time(f"2025-03-20 09:00:00")).dec.deg)
# returns -0.14103777488837554
What am I doing wrong?