<jupyter_start><jupyter_text>Description This notebook computes the dates of the equinox on solar system planets using the Orekit Python wrapper.The seasons are computed by calculating the angle between the Sun vector and the North pole of the planet.  Calculations  The date range for the equinox search must be specified. Be careful of making sure that at least one equinox occurs within this date range.<jupyter_code>from datetime import datetime
date_start = datetime(2020, 1, 1)  # Start of search range
date_end = datetime(2025, 1, 1)  # End of search range
dt = 86400.0  # Time step for the equinox search

planets_list = [
    'MERCURY',
    'VENUS',
    'EARTH',
    'MARS',
    'JUPITER',
    'SATURN',
    'URANUS',
    'NEPTUNE'
]

planet_name = 'MARS'

if planet_name not in planets_list:
    print(f'Error: {planet_name} not a valid planet name')
import orekit
orekit.initVM()

# Modified from https://gitlab.orekit.org/orekit-labs/python-wrapper/blob/master/python_files/pyhelpers.py
from java.io import File
from org.orekit.data import DataProvidersManager, DirectoryCrawler
from orekit import JArray

orekit_data_dir = 'orekit-data'
DM = DataProvidersManager.getInstance()
datafile = File(orekit_data_dir)
if not datafile.exists():
    print('Directory :', datafile.absolutePath, ' not found')

crawler = DirectoryCrawler(datafile)
DM.clearProviders()
DM.addProvider(crawler)
from org.orekit.frames import FramesFactory
from org.orekit.utils import IERSConventions 
ecliptic = FramesFactory.getEcliptic(IERSConventions.IERS_2010)

from org.orekit.time import TimeScalesFactory
utc = TimeScalesFactory.getUTC()

from org.orekit.bodies import CelestialBodyFactory
from org.orekit.utils import PVCoordinatesProvider
planet = CelestialBodyFactory.getBody(planet_name)
body_frame = planet.getBodyOrientedFrame()

sun = CelestialBodyFactory.getSun()
sun_pv_provider = PVCoordinatesProvider.cast_(sun)
import math
# Modified from https://gitlab.orekit.org/orekit-labs/python-wrapper/-/blob/943cae0fe44bedab78137c9fd1263267631e93f5/python_files/pyhelpers.py
# To avoid an exception during a leap second
def absolutedate_to_datetime_no_leap_seconds(orekit_absolutedate):
    """ Converts from orekit.AbsoluteDate objects
    to python datetime objects (utc)"""

    utc = TimeScalesFactory.getUTC()
    or_comp = orekit_absolutedate.getComponents(utc)
    or_date = or_comp.getDate()
    or_time = or_comp.getTime()
    seconds = or_time.getSecond()
    seconds_int = int(math.floor(seconds))
    microseconds = int(1000000.0 * (seconds - math.floor(seconds)))
    if seconds_int > 59:  # This can take the value 60 during a leap second
        seconds_int = 59
        microseconds = 999999  # Also modifying microseconds to ensure that the time flow stays monotonic

    return datetime(or_date.getYear(),
                    or_date.getMonth(),
                    or_date.getDay(),
                    or_time.getHour(),
                    or_time.getMinute(),
                    seconds_int,
                    microseconds)
from orekit.pyhelpers import datetime_to_absolutedate, absolutedate_to_datetime
import pandas as pd
import numpy as np

from org.hipparchus.geometry.euclidean.threed import Vector3D 
planet_pole = Vector3D.PLUS_K

df = pd.DataFrame(columns=['pole_sun_angle_deg'])
df.index.name = 'datetime_utc'

date_current = datetime_to_absolutedate(date_start)
date_end_orekit = datetime_to_absolutedate(date_end)

while date_current.compareTo(date_end_orekit) <= 0:
    sun_position_bf = sun_pv_provider.getPVCoordinates(date_current, body_frame).getPosition()
    sun_pole_angle = np.rad2deg(Vector3D.angle(sun_position_bf, planet_pole))
    df.loc[absolutedate_to_datetime_no_leap_seconds(date_current)] = [sun_pole_angle]
    date_current = date_current.shiftedBy(dt)
from scipy.signal import find_peaks
i_equinoxes, _ = find_peaks(- abs(df['pole_sun_angle_deg'] - 90))
i_spring_equinoxes = []
i_autumn_equinoxes = []
for i in i_equinoxes:
    if i == 0:
        diff = df['pole_sun_angle_deg'][1] - df['pole_sun_angle_deg'][0]
    else:
        diff = df['pole_sun_angle_deg'][i] - df['pole_sun_angle_deg'][i-1]
        
    if diff < 0.0:
        i_spring_equinoxes.append(i)
    else:
        i_autumn_equinoxes.append(i)

print('The spring equinoxes are:')
display(df.iloc[i_spring_equinoxes])

i_summer_solstices, _ = find_peaks(- df['pole_sun_angle_deg'])
print('The summer solstices are:')
display(df.iloc[i_summer_solstices])

print('The autumn equinoxes are:')
display(df.iloc[i_autumn_equinoxes])

i_winter_solstices, _ = find_peaks(df['pole_sun_angle_deg'])
print('The winter solstices are:')
display(df.iloc[i_winter_solstices])<jupyter_output>The spring equinoxes are:
{'fraction_non_alphanumeric': 0.09208572034797369, 'fraction_numerical': 0.016974326331423723, 'mean_word_length': 4.565525383707202, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '39804', 'license': 'Apache-2.0', 'n_tokens_mistral': 1690, 'n_tokens_neox': 1602, 'n_words': 354, 'path': 'seasons-calculation.ipynb', 'repo_name': 'GorgiAstro/planetary-timemachine', 'size': '4713'}