More tests (#27)

* Hour angle tests

* Added vertical dial and doc

* Ignore ephemeris file

* Remove reliance on equinox sun paths

* Add common dial types

* Remove unused code

* Consolidated eccentricity

* Better caching

* Fix doc

* Added orbit shape tests

* Return type hint

* Couple of TODOs

* Ignore debug notebook

---------

Co-authored-by: Russell Goyder <[email protected]>

.gitignore

@@ -134,3 +134,6 @@ dmypy.json
 .octave
 
 todo.txt
+debug.ipynb
+
+*de440s.bsp

src/analemma/geometry.py

@@ -171,6 +171,37 @@ def _within_dialface(vec, dial_length):
 
         return _within_dialface(x, self.x_length) & _within_dialface(y, self.y_length)
 
+    @classmethod
+    def equatorial(cls, latitude: float) -> "DialParameters":
+        """
+        An equatorial dial's face is parallel to the plane of the equator, and the style
+        is therefore perpendicular to both.
+        """
+        theta = (90 - latitude) / 180 * pi
+        return cls(theta=theta, iota=theta, i=theta, d=0)
+
+    @classmethod
+    def horizontal(cls, latitude: float) -> "DialParameters":
+        """
+        A horizontal dial's face is parallel to the ground, with a style aligned
+        with the planet's axis of rotation.
+        """
+        theta = (90 - latitude) / 180 * pi
+        return cls(theta=theta, iota=theta, i=0, d=0)
+
+    @classmethod
+    def vertical(cls, latitude: float) -> "DialParameters":
+        """
+        A vertical dial's face is perpendicular to the ground, and typically faces either
+        south in the northern hemisphere and north in the southern hemisphere. Instances of
+        this class represent a south-facing dial on the equator and in northern latitudes, and
+        a noth-facing dial below the equator (in southern latitudes). Its gnomon
+        is a style, aligned with the planet's axis of rotation.
+        """
+        d = pi if latitude >= 0 else 0
+        theta = (90 - latitude) / 180 * pi
+        return cls(theta=theta, iota=theta, i=pi / 2, d=d)
+
 
 def sin_sunray_dialface_angle(
     t: np.array, planet: orbit.PlanetParameters, dial: DialParameters
@@ -341,7 +372,11 @@ def find_daytime_offsets(planet: orbit.PlanetParameters, dial: DialParameters):
 
 
 class Season(Enum):
-    "Start in Winter because orbit time starts at perihelion"
+    """
+    Enum to encode the four seasons of the year as seen temperate latitudes in the northern hemisphere.
+
+    Start in Winter because orbit time starts at perihelion
+    """
 
     Winter = 0
     Spring = 1

src/analemma/orbit.py

@@ -61,6 +61,19 @@ def __post_init__(self):
         self.T_sd = self.N / (self.N + 1) * self.T_d
         self.Om = pi / self.T_y * self.a
 
+    def clone_with_eccentricity(self, e: float):
+        """
+        Return a new instance with the same parameters but a different eccentricity
+        """
+        return PlanetParameters(
+            N=self.N,
+            T_d=self.T_d,
+            rho=self.rho,
+            alpha=self.alpha,
+            a=self.a,
+            e=e,
+        )
+
     def daily_noons(self) -> np.array:
         """
         Daily time samples in seconds from noon at perihelion
@@ -117,17 +130,16 @@ def earth(cls):
 """
 
 
-def _kepler_params(planet: PlanetParameters = earth, e: float = None):
+def _kepler_params(planet: PlanetParameters = earth):
     a = planet.a
-    if not e:
-        e = planet.e
+    e = planet.e
     b = a * math.sqrt(1 - e**2)  # semi-minor axis
     A = math.sqrt((a + b) / 2)
     B = math.sqrt((a - b) / 2)
     return A, B, planet.Om, planet.T_y
 
 
-def orbital_time(s: npt.ArrayLike, planet: PlanetParameters = earth, e: float = None):
+def orbital_time(s: npt.ArrayLike, planet: PlanetParameters = earth):
     """
     Calculate orbital time given time parameter, t(s)
 
@@ -136,65 +148,47 @@ def orbital_time(s: npt.ArrayLike, planet: PlanetParameters = earth, e: float =
         planet: Planet whose orbit is being analyzed
         e: Optional override for the orbit's eccentricity
     """
-    A, B, Om, T_y = _kepler_params(planet, e)
+    A, B, Om, T_y = _kepler_params(planet)
     return (A**2 + B**2) * s + A * B / Om * sin(2 * Om * s) + T_y / 2
 
 
-_s_finegrained = np.linspace(-pi / earth.Om / 2, pi / earth.Om / 2, 10_000)
-"""
-1-d grid used when inverting the relationship between orbital and spinor time
-"""
-
-
-def _key(e: float) -> int:
-    """
-    The first four significant figures of the given number
-    """
-    return int(10_000 * e)
+_cache_max_size = 50
 
 
-_t_finegrained = {_key(earth.e): orbital_time(_s_finegrained)}
-"""
-Cache of interpolation data used when inverting the relationship between orbital and spinor time 
-"""
+@lru_cache(maxsize=_cache_max_size)
+def _finegrained_interp_points(planet: PlanetParameters):
+    s_points = np.linspace(-pi / planet.Om / 2, pi / planet.Om / 2, 10_000)
+    t_points = orbital_time(s_points, planet)
+    return s_points, t_points
 
 
-def spinor_time(t: npt.ArrayLike, planet: PlanetParameters = earth, e: float = None):
+def spinor_time(t: npt.ArrayLike, planet: PlanetParameters = earth):
     """
     Invert t(s), the relationship of orbital time t with the parameter in the spinor
     treatment of the Kepler problem, s, to give s(t).
-
-    Keep a cache of interpolants, one per eccentricity.
     """
-    if not e:
-        e = planet.e
-    k = _key(e)
-    if k not in _t_finegrained.keys():
-        _t_finegrained[k] = orbital_time(_s_finegrained, planet, e)
-    return np.interp(t, _t_finegrained[k], _s_finegrained)
+    s_points, t_points = _finegrained_interp_points(planet)
+    return np.interp(t, t_points, s_points)
 
 
-def orbital_radius(s: npt.ArrayLike, planet: PlanetParameters = earth, e: float = None):
+def orbital_radius(s: npt.ArrayLike, planet: PlanetParameters = earth):
     """
     Calculate orbital radial coordinate given spinor time parameter, r(s)
     """
-    A, B, Om, _ = _kepler_params(planet, e)
+    A, B, Om, _ = _kepler_params(planet)
     return A**2 + B**2 + 2 * A * B * cos(2 * Om * s)
 
 
-def orbital_angle(s: npt.ArrayLike, planet: PlanetParameters = earth, e=None):
+def orbital_angle(s: npt.ArrayLike, planet: PlanetParameters = earth):
     """
     Calculate orbital angular coordinate given time parameter, phi(s)
     """
-    A, B, Om, _ = _kepler_params(planet, e)
+    A, B, Om, _ = _kepler_params(planet)
     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()
@@ -209,22 +203,27 @@ def _skyfield_season_events(year: int):
 
 
 @dataclass
-class _OrbitDateAndAngle:
+class OrbitDateAndAngle:
+    """
+    Pairing of a date and the corresponding orbit angle
+    """
+
     date: datetime.date
     sigma: float
 
 
-def season_event_info(season_value: int, year: int):
+def season_event_info(season_value: int, year: int) -> OrbitDateAndAngle:
     """
-    TODO also return type
+    Return the date and orbit angle for an equinox or solstice in a given year, identified by
+    the season's value as per [analemma.geometry.Season][].
     """
     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(
+    return OrbitDateAndAngle(
         season_events[skyfield_season_value].utc_datetime().date(),
         season_event_angles[skyfield_season_value],
     )

src/analemma/plot.py

@@ -85,7 +85,7 @@ def plot_analemma_season_segment(
 
     times = _analemma_plot_sampling_times(season, hour_offset, planet, dial)
     if times.size == 0:
-        return ax
+        return []
     return _plot_analemma_segment(
         ax,
         times,
@@ -180,6 +180,14 @@ def plot_season_event_sun_path(
 
     season_event = orbit.season_event_info(season.value, year)
 
+    # for an equatorial dial on the equinoxes, the sun ray is parallel to the dial face
+    if (
+        season.name in ("Spring", "Autumn")
+        and abs(dial.d) < 1.0e-5
+        and abs(dial.theta - dial.i) < 0.25 / 180 * pi
+    ):
+        return []
+
     orbit_day = orbit.orbit_date_to_day(season_event.date)
     day_type = _determine_day_type(planet, dial, orbit_day)
     if day_type == DayType.SunNeverRises:
@@ -335,21 +343,6 @@ def _analemma_point_coordinates(
     return x, y, on_dial, above_dial_plane
 
 
-_sun_times_cache = {}
-
-
-def _get_sun_times(planet: orbit.PlanetParameters, dial: geom.DialParameters):
-    key = (id(planet), id(dial))
-    if key not in _sun_times_cache.keys():
-        _sun_times_cache[key] = [
-            geom.find_sun_rise_noon_set_relative_to_dial_face(
-                days_since_perihelion, planet, dial
-            )
-            for days_since_perihelion in np.arange(0, 365)
-        ]
-    return _sun_times_cache[key]
-
-
 def _solstice_days(planet: orbit.PlanetParameters, dial: geom.DialParameters):
     _, sines = geom.sunray_dialface_angle_over_one_year(planet, dial)
     return (np.argmax(sines), np.argmin(sines))
@@ -460,6 +453,19 @@ def hour_offset_to_oclock(hour_offset: int):
         raise Exception(f"hour_offset {hour_offset} doesn't seem to be a number")
 
 
+def _font_size(ax: Axes) -> str:
+    "Map bounding box width in inches to font size string"
+    bbox = ax.get_window_extent().transformed(
+        ax.get_figure().dpi_scale_trans.inverted()
+    )
+    if bbox.width >= 6:
+        return "small"
+    elif bbox.width >= 4:
+        return "x-small"
+    else:
+        return "xx-small"
+
+
 def annotate_analemma_with_hour(
     ax: Axes,
     hour_offset: int,
@@ -481,11 +487,34 @@ def annotate_analemma_with_hour(
                 xytext=ptext,
                 arrowprops={"arrowstyle": "-"},
                 horizontalalignment="center",
-                fontsize="small",
+                fontsize=_font_size(ax),
             )
     return None
 
 
+def _reorder_legend_info(legend_info):
+    """
+    Move Winter to the end of the list
+    """
+    if legend_info[0][1] == "Winter":
+        winter_entry = legend_info.pop(0)
+        legend_info.append(winter_entry)
+    return legend_info
+
+
+def _adjust_for_hemisphere(dial, labels):
+    if dial.theta < pi / 2:
+        return labels
+    else:
+        mapping = {
+            "Winter": "Summer",
+            "Spring": "Autumn",
+            "Summer": "Winter",
+            "Autumn": "Spring",
+        }
+        return [mapping[label] for label in labels]
+
+
 def plot_hourly_analemmas(
     ax: Axes,
     planet: orbit.PlanetParameters,
@@ -497,44 +526,44 @@ def plot_hourly_analemmas(
     """
     Plot one analemma for each hour as seen on the face of a sundial
 
-    This function plots several analemmas, one per hour of daytime. The line style shows the season. One line
-    showing the path of the shadow tip during the day for each solstice is also shown (with line style appropriate to
-    the season) and forms an envelope marking the longest shadows in Winter and the shortest shadows in Summer.
-    Similarly, the path of the shadow tip on each equinox is shown and appears as a straight line. Moreover, the two
-    straight lines fall on top of each other.
+    This function plots several analemmas, one per hour of daytime. The line style shows the season. One line showing
+    the path of the shadow tip during the day for each solstice is also shown (with line style appropriate to the
+    season) and on a horizontal dial forms an envelope marking the longest shadows in Winter and the shortest shadows in
+    Summer. Similarly, the path of the shadow tip on each equinox is shown and appears as a straight line. Moreover, the
+    two straight lines fall on top of each other.
+
+    In the legend, the seasons are labelled according to the hemisphere in which the sundial is located, so that for
+    sundials in the southern hemisphere, summer occurs in December.
 
     Parameters:
-        ax: matplotlib axes
-        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)
+        ax: matplotlib axes 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)
 
-    lines_for_legend = []
-    seasons = []
-    count = 0
+    legend_info = []
     for season in geom.Season:
+        segment_lines = []
         for hour in hour_offsets:
-            plot_analemma_season_segment(
+            segment_lines += plot_analemma_season_segment(
                 ax, season, hour, planet, dial, linewidth=0.75, **kwargs
             )
             annotate_analemma_with_hour(ax, hour, planet, dial)
-        lines = plot_season_event_sun_path(
+
+        if len(segment_lines) > 0:
+            legend_info.append((segment_lines[0], season.name))
+
+        plot_season_event_sun_path(
             ax, season, planet, dial, linewidth=0.75, year=year, **kwargs
         )
-        if len(lines) > 0:
-            lines_for_legend += lines
-            seasons += [count]
-            count += 1
 
     # put a circle at the base of the gnomon
     ax.plot(0, 0, "ok")
 
-    # reorder seasons
-    ordered_lines = [lines_for_legend[s] for s in seasons]
-    ax.legend(handles=ordered_lines)
+    handles, labels = zip(*_reorder_legend_info(legend_info))
+    labels = _adjust_for_hemisphere(dial, labels)
+    ax.legend(handles, labels, fontsize=_font_size(ax))
 
     if title:
         ax.set_title(title)

src/analemma/tests/test_geometry.py

@@ -1,2 +1,124 @@
-def test_trivial():
-    pass
+import pytest
+import numpy as np
+from sympy import utilities as util
+from analemma import geometry as geom
+from analemma.algebra import frame, result
+
+
+pi = np.pi
+hour_angle_args = ["alpha", "sigma", "psi", "iota", "theta"]
+
+
+@pytest.mark.parametrize(
+    ",".join(hour_angle_args),
+    [
+        (0, 0, 0, 0, 0),
+        (23.5 / 180 * pi, 0, 0, 0, 0),  # add Earth axis tilt
+        (pi / 2, 0, 0, 0, 0),  # 90 degree tilt
+        (
+            23.5 / 180 * pi,
+            pi / 2,
+            0.1,
+            0,
+            0,
+        ),  # September equinox (non-zero psi to avoid tan singularity)
+        (23.5 / 180 * pi, pi, 0, 0, 0),  # December solstice
+        (
+            23.5 / 180 * pi,
+            3 * pi / 2,
+            0.1,
+            0,
+            0,
+        ),  #  March equinox (non-zero psi to avoid tan singularity)
+        (23.5 / 180 * pi, 2 * pi, 0, 0, 0),  # June solstice
+        (23.5 / 180 * pi, pi / 7, 71 / 180 * pi, 0, 0),  # arbitrary psi
+        (23.5 / 180 * pi, pi / 7, 71 / 180 * pi, 0, pi / 2),  # equator
+        (23.5 / 180 * pi, pi / 7, 71 / 180 * pi, 0, pi),  # South pole
+        (
+            23.5 / 180 * pi,
+            pi / 7,
+            71 / 180 * pi,
+            0,
+            68 / 180 * pi,
+        ),  # arbitrary latitude
+        (
+            23.5 / 180 * pi,
+            pi / 7,
+            71 / 180 * pi,
+            pi / 6,
+            68 / 180 * pi,
+        ),  # gnomon with 30 degree tilt
+        (
+            23.5 / 180 * pi,
+            pi / 7,
+            71 / 180 * pi,
+            79 / 180 * pi,
+            68 / 180 * pi,
+        ),  # gnomon almost touching dial face
+    ],
+)
+def test_hour_angle(alpha, sigma, psi, iota, theta):
+    """
+    Tests for the implementation of the hour angle calculation
+
+    See [analemma.geometry.hour_angle_terms][] and [analemma.geometry.hour_angle][]
+    """
+
+    # alegbraic result
+    g = frame.gnomon("e", zero_decl=True)
+    s = frame.sunray()
+    S = result.shadow_bivector(s, g)
+    M = frame.meridian_plane()
+    sinXi_sin_mu, sinXi_cos_mu = result.hour_angle_sincos(S, M)
+
+    # convert to numerical functions
+    sin_term = util.lambdify(hour_angle_args, sinXi_sin_mu)
+    cos_term = util.lambdify(hour_angle_args, sinXi_cos_mu)
+
+    # evaluate and ensure consistent with implementation in geometry module
+    sin_val, cos_val = geom.hour_angle_terms(alpha, sigma, psi, iota - theta)
+    assert sin_term(alpha, sigma, psi, iota, theta) == pytest.approx(sin_val)
+    assert cos_term(alpha, sigma, psi, iota, theta) == pytest.approx(cos_val)
+
+    # ensure depends only on difference between iota and theta
+    arbitrary_shift = pi / 13
+    assert sin_term(alpha, sigma, psi, iota, theta) == pytest.approx(
+        sin_term(alpha, sigma, psi, iota + arbitrary_shift, theta + arbitrary_shift)
+    )
+    assert cos_term(alpha, sigma, psi, iota, theta) == pytest.approx(
+        cos_term(alpha, sigma, psi, iota + arbitrary_shift, theta + arbitrary_shift)
+    )
+
+    # test some periodicity
+    assert geom.hour_angle_terms(alpha, sigma, psi, iota - theta) == pytest.approx(
+        geom.hour_angle_terms(
+            alpha + 2 * pi, sigma + 4 * pi, psi + 10 * pi, iota - theta + 50 * pi
+        )
+    )
+
+    # combine the terms
+    assert np.tan(geom.hour_angle(alpha, sigma, psi, iota - theta)) == pytest.approx(
+        sin_val / cos_val
+    )
+
+
+@pytest.mark.parametrize(("latitude"), [0, 40, 50, 90, -10, -90, -0])
+def test_common_dial_types_basic_construction(latitude: float):
+    eq = geom.DialParameters.equatorial(latitude=latitude)
+    theta = (90 - latitude) / 180 * pi
+    assert eq.theta == pytest.approx(theta)
+    assert eq.theta == pytest.approx(eq.iota) == pytest.approx(eq.i)
+    assert pytest.approx(eq.d) == pytest.approx(0)
+
+    hor = geom.DialParameters.horizontal(latitude=latitude)
+    theta = (90 - latitude) / 180 * pi
+    assert hor.theta == pytest.approx(theta)
+    assert hor.theta == pytest.approx(hor.iota)
+
+    vert = geom.DialParameters.vertical(latitude=latitude)
+    theta = (90 - latitude) / 180 * pi
+    assert vert.theta == pytest.approx(theta)
+    assert vert.theta == pytest.approx(eq.iota)
+    assert vert.i == pytest.approx(pi / 2)
+    decl = pi if latitude >= 0 else 0
+    assert vert.d == pytest.approx(decl)

src/analemma/tests/test_orbit.py

@@ -1,6 +1,8 @@
 import pytest
 import datetime
+import math
 import numpy as np
+from numpy import sin, cos
 from skyfield import almanac
 from skyfield.api import load
 from skyfield import searchlib as sf_search
@@ -97,3 +99,54 @@ def test_season_event_day_diffs():
     diffs = [(days[0] - days[3]) % 365] + list(diffs)  # days from 353 to 77
     for diff, ans in zip(diffs, [89, 92, 94, 90]):
         assert diff == ans
+
+
+def test_circular_orbit_shape():
+    """
+    Ensure that an orbit with zero eccentricity is circular
+    """
+    circ = orbit.PlanetParameters(
+        N=100,
+        T_d=24 * 3600,
+        rho=0.0,
+        alpha=0.0,
+        a=17.0,
+        e=0.0,
+    )
+
+    t_integers = np.arange(int(circ.N))
+    t = circ.T_d * t_integers
+    s = orbit.spinor_time(t, circ)
+
+    x = orbit.orbital_radius(s, circ) * cos(orbit.orbital_angle(s, circ))
+    y = orbit.orbital_radius(s, circ) * sin(orbit.orbital_angle(s, circ))
+
+    assert np.allclose(x**2 + y**2, circ.a**2)
+
+
+def test_elliptical_orbit_shape():
+    """
+    Ensure that an orbit with non-zero eccentricity is elliptical
+    """
+    ell = orbit.PlanetParameters(
+        N=100,
+        T_d=24 * 3600,
+        rho=0.0,
+        alpha=0.0,
+        a=17.0,
+        e=0.5,
+    )
+
+    t_integers = np.arange(int(ell.N))
+    t = ell.T_d * t_integers
+    s = orbit.spinor_time(t, ell)
+
+    x = orbit.orbital_radius(s, ell) * cos(orbit.orbital_angle(s, ell))
+    y = orbit.orbital_radius(s, ell) * sin(orbit.orbital_angle(s, ell))
+
+    xx = x + ell.a * ell.e  # shift origin from focus to center of ellipse
+
+    a = ell.a
+    b = ell.a * math.sqrt(1 - ell.e**2)
+
+    assert np.allclose((xx / a) ** 2 + (y / b) ** 2, 1.0)

src/analemma/tests/test_results.py

@@ -133,7 +133,7 @@ def test_shadow_bivector_magnitude_angle_cos():
     r"""
     Compare two different ways of calculating $\cos(\Xi)$
 
-    Ensure that $s\cdot g$ as computed via symbolic is equal to an independent derivation
+    Ensure that $s\cdot g$ as computed via symbolic algebra is equal to an independent derivation
     """
 
     s = frame.sunray()
@@ -303,8 +303,6 @@ def test_gnomon_shadow_projection():
     Check that $1 + \lambda\cos(\Xi)$ is equal to the projection of $w$ onto $g$
     """
 
-    # TODO fixtures
-
     gn = frame.gnomon("n", zero_decl=True)
     p = sp.Symbol("p")