23 make date handling for equinoxes solstices and apsides more robust (…

…#24) * Add skyfield for season events * Add skyfield for perihelion * Season boundaries from skyfield * Links * Cleaner * Quieter * Changelog --------- Co-authored-by: Russell Goyder <[email protected]>
russellgoyder · Oct 28, 2024 · 31d9c59 · 31d9c59
1 parent 78f0c5c
commit 31d9c59
Show file tree

Hide file tree

Showing 12 changed files with 227 additions and 69 deletions.
diff --git a/docs/changelog.md b/docs/changelog.md
@@ -1,7 +1,11 @@
 
 # Changelog
 
-## Latest
+## v0.2
+
+Integated [Skyfield](https://rhodesmill.org/skyfield/) for season events and apsides.
+
+## v0.1
 
 Generalized `hour_offset` to float in order to plot analemma at any time of day.
 

diff --git a/docs/index.md b/docs/index.md
@@ -43,7 +43,7 @@ ap.plot_analemma(ax, earth, vertical_dial)
 See [Analemma Plots](nb/sundial_plots.md) for complete examples showing various analemmas.
 
 For the connection between the angle of the sun, the date, and the time, see [The Equation of
-Time](nb/equation_of_time.md) and [Orbit Analysis](nb/orbit_analysis.md).
+Time](nb/equation_of_time.md), [Sunrise and Sunset](nb/sunrise_and_sunset.md), and [Orbit Analysis](nb/orbit_analysis.md).
 
 ## Background
 
@@ -58,6 +58,7 @@ The analemma is the path traced by the shadow on a sundial (or the sun in the sk
  * [The Shadow Angle](nb/shadow_angle.md): unit vector parallel to the shadow and how it moves in time
  * [The Shadow Length](nb/shadow_length.md): the length of the shadow
  * [The Analemma](nb/analemma.md): parametric expression for the analemma as the coordinates of the tip of a sundial's shadow
+ * [Comparison with Rohr's book](nb/rohr_comparison.md): comparison with known results from a standard text
 
 
 ## Project Links

diff --git a/docs/nb/sundial_plots.ipynb b/docs/nb/sundial_plots.ipynb
@@ -95,7 +95,7 @@
     "ax.grid()\n",
     "ax.axis(\"equal\")\n",
     "\n",
-    "aplot.plot_hourly_analemmas(ax, earth, camdial, \"Analemmatic sundial in Cambridge, UK\")"
+    "aplot.plot_hourly_analemmas(ax, earth, camdial, year=2024, title=\"Analemmatic sundial in Cambridge, UK\")"
    ]
   },
   {
@@ -130,7 +130,7 @@
     "ax.grid()\n",
     "ax.axis(\"equal\")\n",
     "\n",
-    "aplot.plot_hourly_analemmas(ax, earth, tropicaldial, \"Tropical sundial with unusual geometry\")"
+    "aplot.plot_hourly_analemmas(ax, earth, tropicaldial, year=2024, title=\"Tropical sundial with unusual geometry\")"
    ]
   },
   {
@@ -165,7 +165,7 @@
     "ax.grid()\n",
     "ax.axis(\"equal\")\n",
     "\n",
-    "aplot.plot_hourly_analemmas(ax, earth, arcticdial, \"Arctic sundial\")"
+    "aplot.plot_hourly_analemmas(ax, earth, arcticdial, year=2024, title=\"Arctic sundial\")"
    ]
   }
  ],

diff --git a/docs/nb/sunrise_and_sunset.ipynb b/docs/nb/sunrise_and_sunset.ipynb
@@ -52,7 +52,7 @@
     "earth = orbit.PlanetParameters.earth()\n",
     "camdial = geom.DialParameters(theta=37.5/180*pi, iota=37.5/180*pi, i=0, d=0) # Analemmatic dial in Cambridge, UK\n",
     "\n",
-    "aplot.plot_sunrise_sunset(ax, geom.orbit_day_to_date(50), earth, camdial)"
+    "aplot.plot_sunrise_sunset(ax, orbit.orbit_day_to_date(50), earth, camdial)"
    ]
   },
   {

diff --git a/pyproject.toml b/pyproject.toml
@@ -24,6 +24,7 @@ keywords = [
 dependencies = [
     "numpy",
     "scipy",
+    "skyfield",
     "matplotlib",
     "sympy",
     "galgebra",

diff --git a/scripts/clean.sh b/scripts/clean.sh
@@ -4,3 +4,4 @@ git checkout docs/nb/*.ipynb # undo scripts/dollar_dollar.py
 rm -rf docs/nb/*_files
 rm -f docs/nb/*.md
 find . -name .DS_Store -delete
+rm -f *.bsp docs/nb/*.bsp # ephemeris files from skyfield
diff --git a/scripts/run_tests.sh b/scripts/run_tests.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+pytest -n auto src --nbmake docs/nb/*.ipynb
diff --git a/src/analemma/geometry.py b/src/analemma/geometry.py
@@ -7,7 +7,6 @@
 from scipy import optimize as sci_opt
 from typing import Tuple
 from dataclasses import dataclass
-import datetime
 from enum import Enum
 from analemma import orbit
 
@@ -302,30 +301,6 @@ def find_sun_rise_noon_set_relative_to_dial_face(
     return SunTimes(t_sunrise, t_noon, t_sunset, days_since_perihelion)
 
 
-def orbit_day_to_date(orbit_day: int, year=2024) -> datetime.date:
-    """
-    Convert from the number of days since perihelion to the date
-
-    Note that this varies from year to year and this implementation is only exact for 2024
-    """
-    perihelion_date = datetime.date.fromisoformat(
-        f"{year}-01-03"
-    )  # approximately true for other years
-    return perihelion_date + datetime.timedelta(days=int(orbit_day))
-
-
-def orbit_date_to_day(the_date: datetime.date, year=2024) -> int:
-    """
-    Convert from the date to the number of days since perihelion
-
-    Note that this varies from year to year and this implementation is only exact for 2024
-    """
-    perihelion_date = datetime.date.fromisoformat(
-        f"{year}-01-03"
-    )  # approximately true for other years
-    return (the_date - perihelion_date).days
-
-
 def sunray_dialface_angle(
     planet: orbit.PlanetParameters,
     dial: DialParameters,

diff --git a/src/analemma/orbit.py b/src/analemma/orbit.py
@@ -6,6 +6,11 @@
 import numpy as np
 from numpy import sin, cos, typing as npt
 from dataclasses import dataclass, field
+import datetime
+from functools import lru_cache
+from skyfield import almanac
+from skyfield.api import load
+from skyfield import searchlib as sf_search
 
 
 pi = math.pi
@@ -185,3 +190,79 @@ def orbital_angle(s: npt.ArrayLike, planet: PlanetParameters = earth, e=None):
     tanSigY = (A**2 - B**2) * sin(2 * Om * s)
     tanSigX = (A**2 + B**2) * cos(2 * Om * s) + 2 * A * B
     return np.arctan2(tanSigY, tanSigX) + pi
+
+
+_cache_max_size = 50
+
+
+@lru_cache(maxsize=_cache_max_size)
+def _skyfield_ephemeris():
+    return load("de440s.bsp"), load.timescale()
+
+
+def _skyfield_season_events(year: int):
+    eph, ts = _skyfield_ephemeris()
+    jan1 = ts.utc(year, 1, 1)
+    dec31 = ts.utc(year, 12, 31)
+    event_times, _ = almanac.find_discrete(jan1, dec31, almanac.seasons(eph))
+    return event_times
+
+
+@dataclass
+class _OrbitDateAndAngle:
+    date: datetime.date
+    sigma: float
+
+
+def season_event_info(season_value: int, year: int):
+    """
+    TODO also return type
+    """
+    season_events = _skyfield_season_events(year)
+    # S S A W (seasons)
+    # 1 2 3 0 (Season enum)
+    # 0 1 2 3 (Skyfield)
+    skyfield_season_value = (season_value + 3) % 4
+    season_event_angles = [pi / 2, 0, 3 * pi / 2, pi]
+    return _OrbitDateAndAngle(
+        season_events[skyfield_season_value].utc_datetime().date(),
+        season_event_angles[skyfield_season_value],
+    )
+
+
+def _earth_distance(time):
+    eph, _ = _skyfield_ephemeris()
+    earth = eph["earth"]
+    sun = eph["sun"]
+    e = earth.at(time)
+    s = e.observe(sun)
+    return s.distance().km
+
+
+_earth_distance.step_days = 1
+
+
+def _perihelion_date(year: int) -> datetime.date:
+    _, ts = _skyfield_ephemeris()
+    start_time = ts.utc(year, 1, 1)
+    end_time = ts.utc(year + 1, 1, 1)
+    perihelion = sf_search.find_minima(start_time, end_time, _earth_distance)
+    return perihelion[0][0].utc_datetime().date()
+
+
+def orbit_day_to_date(orbit_day: int, year: int = None) -> datetime.date:
+    """
+    Convert from the number of days since perihelion to the date
+    """
+    if not year:
+        year = datetime.date.today().year
+    return _perihelion_date(year) + datetime.timedelta(days=int(orbit_day))
+
+
+def orbit_date_to_day(the_date: datetime.date, year: int = None) -> int:
+    """
+    Convert from the date to the number of days since perihelion
+    """
+    if not year:
+        year = datetime.date.today().year
+    return (the_date - _perihelion_date(year)).days
diff --git a/src/analemma/plot.py b/src/analemma/plot.py
@@ -3,7 +3,6 @@
 """
 
 import numpy as np
-from dataclasses import dataclass
 import datetime
 from enum import Enum
 from typing import TypeVar
@@ -34,7 +33,7 @@ def _analemma_plot_sampling_times(
     planet: orbit.PlanetParameters,
     dial: geom.DialParameters,
 ):
-    # season lengths are [89, 91, 94, 91] (Winter Spring Summer Autumn)
+    # season lengths are [89, 92, 94, 90] (Winter Spring Summer Autumn) in 2024
     # place equinoxes and solstices in the middle for plotting
     season_boundaries = [0, 44, 135, 229, 320]
     if season != geom.Season.Winter:
@@ -129,28 +128,6 @@ def plot_analemma(
     )
 
 
-@dataclass
-class _OrbitDateAndAngle:
-    date: datetime.date
-    sigma: float
-
-
-_equinox_or_solstice_info = {
-    geom.Season.Summer.value: _OrbitDateAndAngle(
-        datetime.date.fromisoformat("2024-06-20"), 0
-    ),
-    geom.Season.Spring.value: _OrbitDateAndAngle(
-        datetime.date.fromisoformat("2024-03-20"), pi / 2
-    ),
-    geom.Season.Winter.value: _OrbitDateAndAngle(
-        datetime.date.fromisoformat("2024-12-21"), pi
-    ),
-    geom.Season.Autumn.value: _OrbitDateAndAngle(
-        datetime.date.fromisoformat("2024-09-22"), 3 * pi / 2
-    ),
-}
-
-
 class DayType(Enum):
     SunNeverRises = 0
     SunNeverSets = 1
@@ -177,20 +154,33 @@ def _determine_day_type(
         return DayType.SunRisesAndSets
 
 
-def plot_special_sun_path(
+def plot_season_event_sun_path(
     ax: Axes,
     season: geom.Season,
     planet: orbit.PlanetParameters,
     dial: geom.DialParameters,
+    year: int = None,
     **kwargs,
 ):
     """
     Plot the path of the sun across the dial on the equinox or solstice in the given season
+
+    Parameters:
+        ax: matplotlib axes
+        season: The given season
+        planet: The planet on which the dial is located
+        dial: The orientation and location of the sundial
+        year: The year in which the seasons events fall (defaults to current year)
     """
 
     num_times = 1000
 
-    orbit_day = geom.orbit_date_to_day(_equinox_or_solstice_info[season.value].date)
+    if not year:
+        year = datetime.today().year
+
+    season_event = orbit.season_event_info(season.value, year)
+
+    orbit_day = orbit.orbit_date_to_day(season_event.date)
     day_type = _determine_day_type(planet, dial, orbit_day)
     if day_type == DayType.SunNeverRises:
         return []
@@ -213,7 +203,7 @@ def plot_special_sun_path(
 
     psis = planet.rotation_angle(times)
 
-    sigma = _equinox_or_solstice_info[season.value].sigma
+    sigma = season_event.sigma
     x_raw, y_raw = geom.shadow_coords_xy(
         planet.alpha, sigma, psis, dial.iota, dial.theta, dial.i, dial.d
     )
@@ -248,7 +238,7 @@ def plot_sunrise_sunset(
         planet: The planet on which the dial is located
         dial: The orientation and location of the sundial
     """
-    orbit_day = geom.orbit_date_to_day(date)
+    orbit_day = orbit.orbit_date_to_day(date)
     day_type = _determine_day_type(planet, dial, orbit_day)
     if not day_type == DayType.SunRisesAndSets:
         raise Exception(
@@ -501,6 +491,7 @@ def plot_hourly_analemmas(
     planet: orbit.PlanetParameters,
     dial: geom.DialParameters,
     title: str = None,
+    year: int = None,
     **kwargs,
 ):
     """
@@ -517,6 +508,7 @@ def plot_hourly_analemmas(
         planet: The planet on which the dial is located
         dial: The orientation and location of the sundial
         title: Title to add to the axes
+        year: Year for which the plot hourly analemmas (defaults to current year)
     """
     hour_offsets = geom.find_daytime_offsets(planet, dial)
 
@@ -529,8 +521,8 @@ def plot_hourly_analemmas(
                 ax, season, hour, planet, dial, linewidth=0.75, **kwargs
             )
             annotate_analemma_with_hour(ax, hour, planet, dial)
-        lines = plot_special_sun_path(
-            ax, season, planet, dial, linewidth=0.75, **kwargs
+        lines = plot_season_event_sun_path(
+            ax, season, planet, dial, linewidth=0.75, year=year, **kwargs
         )
         if len(lines) > 0:
             lines_for_legend += lines
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,3 @@
		#!/bin/bash

		pytest -n auto src --nbmake docs/nb/*.ipynb