Test satellite properties

CHANGELOG.md

@@ -4,6 +4,10 @@ All notable changes to this project will be documented here.
 
 # Unreleased
 - Certain orbital parameters can be displayed in the plot
+- Implemented automatic testing
+- Tests architecture rewritten
+- Wrote invariants unit tests
+- PerifocalOrbitEq refactored to live in the Satellite class
 
 # 0.4.2 - 08/02/2026
 - Refactored config builder and controller into variable, property and display builders/controllers

src/orbit_visualiser/core/orbit.py

@@ -20,7 +20,6 @@ def y(self) -> Callable[[float], float]:
         return self._y
 
 
-# TODO: self.orbital_angles should return just a tuple, not the linspace. The linspace is only for plotting so logic for calculating it should be in OrbitFigure.
 class Orbit():
 
 
@@ -86,12 +85,11 @@ def is_closed(self) -> bool:
     def orbit_eq(self) -> PerifocalOrbitEq:
         return PerifocalOrbitEq(self._e, self._p)
 
-    def orbital_angles(self):
+    def orbital_angles(self) -> tuple[float]:
         if self._e < 1:
-            return np.linspace(0, 2*pi, 100_000)
+            return 0, 2*pi
 
-        delta = 0.0001
-        return np.linspace(-self._t_asymp + delta, self._t_asymp - delta, 100_000)
+        return -self._t_asymp, self._t_asymp
 
     def update_orbital_properties(self):
         e, rp = self._e, self._rp

src/orbit_visualiser/core/satellite.py

@@ -95,7 +95,7 @@ def update_satellite_properties(self) -> None:
         self._h = h
 
         # Helper quantities
-        den = 1 + e*np.cos(nu)
+        den = 1 + e*np.cos(nu) # Denominator of the orbit equation
         mu_over_h = mu/h
 
         # Geometry

src/orbit_visualiser/ui/figure/orbit_figure.py

@@ -3,12 +3,15 @@
 from matplotlib.figure import Figure
 from matplotlib.axes import Axes
 from matplotlib.patches import Circle
+import numpy as np
+from numpy.typing import NDArray
 from orbit_visualiser.core import Orbit, CentralBody, Satellite
 
 # TODO: Fix bug where scroll zoom doesn't register as changing the view so the native matplotlib home button has unexpected (and often undesirable) behaviour.
 class OrbitFigure():
 
     DISPLAY_TEXT_OFFSET = (1.5, 1.5)
+    NUM_POINTS = 100_000
 
     def __init__(
             self,
@@ -75,7 +78,7 @@ def _configure_axes(self) -> None:
 
     def _initialise_plot(self) -> None:
         # Plot the initial orbit
-        t = self._orbit.orbital_angles()
+        t = self._get_anomaly_data(self._orbit.orbital_angles())
         orbit_eq = self._orbit.orbit_eq
         x, y = orbit_eq.x, orbit_eq.y
         self._line, = self._ax.plot(x(t) , y(t), color = "#2F2F2F", alpha = 0.5, linewidth = 1.5)
@@ -105,8 +108,17 @@ def _build_toolbar(self) -> None:
         toolbar.update()
         toolbar.pack(side = "bottom", fill = "x")
 
+    def _get_anomaly_data(self, angles: tuple[float]) -> NDArray[np.float64]:
+        lower_lim, upper_lim = angles
+
+        # Check if angles represent closed or open orbit (open orbits will have negative angles).
+        # If the orbit is open then we need a small offset so the plot doesn't evaluate to infinity
+        # and cause runtime errors or unusual graphical artifacts.
+        delta = 0.0001 if lower_lim < 0 else 0
+        return np.linspace(lower_lim + delta, upper_lim - delta, OrbitFigure.NUM_POINTS)
+
     def redraw_orbit(self) -> None:
-        t = self._orbit.orbital_angles()
+        t = self._get_anomaly_data(self._orbit.orbital_angles())
         orbit_eq = self._orbit.orbit_eq
         self._line.set_data(orbit_eq.x(t), orbit_eq.y(t))
 
@@ -131,7 +143,13 @@ def plot_periapsis_point(self) -> None:
         self._rp_point, = self._ax.plot(
             self._orbit.rp, 0, ms = 3, marker = "o", zorder = 9, color = "#502BF2", label = "$r_p$"
         )
-        self._rp_annotation = self._ax.annotate("$r_p$", xy = (self._orbit.rp, 0), xycoords = "data", xytext = self.DISPLAY_TEXT_OFFSET, textcoords = "offset points")
+        self._rp_annotation = self._ax.annotate(
+            "$r_p$",
+            xy = (self._orbit.rp, 0),
+            xycoords = "data",
+            xytext = OrbitFigure.DISPLAY_TEXT_OFFSET,
+            textcoords = "offset points"
+        )
 
     @staticmethod
     def _zoom_factory(ax: Axes, base_scale = 2.):

tests/conftest.py

@@ -0,0 +1,42 @@
+import pytest
+import numpy as np
+from numpy.typing import NDArray
+from math import pi
+from typing import Callable
+from orbit_visualiser.core import Satellite, Orbit, CentralBody
+
+# ---------- True anomaly grids --------------------
+@pytest.fixture(scope = "session")
+def closed_anomaly_grid() -> NDArray[np.float64]:
+    return np.linspace(0, 2*pi, 20)
+
+@pytest.fixture(scope = "session")
+def open_anomaly_grid() -> Callable[[Orbit, int], NDArray[np.float64]]:
+    def _create_grid(orbit: Orbit, num: int = 50) -> NDArray[np.float64]:
+        orbit_angles = orbit.orbital_angles()
+        return np.linspace(orbit_angles[0], orbit_angles[1], num)
+    return _create_grid
+
+# ---------- Orbital object factories --------------
+@pytest.fixture
+def orbit_factory() -> Callable[[float, float], Orbit]:
+    def _create(e: float = 0.0, rp: float = 50_000.0) -> Orbit:
+        return Orbit(e, rp)
+    return _create
+
+@pytest.fixture
+def central_body_factory() -> Callable[[float], CentralBody]:
+    def _create(mu: float = 398_600.0) -> CentralBody:
+        return CentralBody(mu)
+    return _create
+
+@pytest.fixture
+def satellite_factory(
+    orbit_factory: Callable[[float, float], Orbit],
+    central_body_factory: Callable[[float], CentralBody]
+) -> Callable[[float, float, float], Satellite]:
+    def _create(e: float = 0.0, rp: float = 50_000.0, mu: float = 398_600.0) -> Satellite:
+        orbit: Orbit = orbit_factory(e, rp)
+        central_body: CentralBody = central_body_factory(mu)
+        return Satellite(orbit, central_body)
+    return _create

tests/invariants/circular_orbit_test.py

@@ -0,0 +1,77 @@
+import pytest
+import numpy as np
+from numpy.typing import NDArray
+from typing import Callable
+from orbit_visualiser.core import Orbit, Satellite
+from tests.test_cases import rp_test_cases, rp_mu_test_cases
+
+@pytest.mark.parametrize("rp", rp_test_cases)
+def test_circular_orbit_distance_parameters(orbit_factory, rp: float):
+    """
+    For circular orbits the radius of periapsis, radius of apoapsis, semi-major axis, semi-minor axis
+    and orbital parameter should all be equivalent.
+    """
+    orbit: Orbit = orbit_factory(rp = rp)
+
+    periapsis = orbit.rp
+    distances = [
+        orbit.ra,
+        orbit.a,
+        orbit.b,
+        orbit.p
+    ]
+    assert np.allclose(distances, periapsis)
+
+@pytest.mark.parametrize("rp, mu", rp_mu_test_cases)
+def test_circular_orbit_flight_angle(
+    satellite_factory: Callable[[float, float, float], Satellite],
+    closed_anomaly_grid: NDArray[np.float64],
+    rp: float,
+    mu: float):
+    """
+    For circular orbits the flight angle of a satellite at any true anomaly should be 0.
+    """
+    satellite: Satellite = satellite_factory(rp = rp, mu = mu)
+
+    for nu in closed_anomaly_grid:
+        satellite.nu = nu
+        satellite.update_satellite_properties()
+
+        gam = satellite.gam
+        assert np.isclose(gam, 0)
+
+@pytest.mark.parametrize("rp, mu", rp_mu_test_cases)
+def test_circular_orbit_radial_velocity(
+    satellite_factory: Callable[[float, float, float], Satellite],
+    closed_anomaly_grid: NDArray[np.float64],
+    rp: float,
+    mu: float):
+    """
+    For circular orbits the radial velocity of a satellite at any true anomaly should be 0.
+    """
+    satellite: Satellite = satellite_factory(rp = rp, mu = mu)
+
+    for nu in closed_anomaly_grid:
+        satellite.nu = nu
+        satellite.update_satellite_properties()
+
+        v_rad = satellite.v_radial
+        assert np.isclose(v_rad, 0)
+
+@pytest.mark.parametrize("rp, mu", rp_mu_test_cases)
+def test_circular_orbit_anomalies(
+    satellite_factory: Callable[[float, float, float], Satellite],
+    closed_anomaly_grid: NDArray[np.float64],
+    rp: float,
+    mu: float):
+    """
+    For circular orbits the eccentric, mean and true anomalies should always be equivalent.
+    """
+    satellite: Satellite = satellite_factory(rp = rp, mu = mu)
+
+    for nu in closed_anomaly_grid:
+        satellite.nu = nu
+        satellite.update_satellite_properties()
+
+        anomalies = [satellite.m_anomaly, satellite.e_anomaly]
+        assert np.allclose(anomalies, nu)

tests/invariants/energy_test.py

@@ -0,0 +1,64 @@
+import pytest
+import numpy as np
+from typing import Callable
+from numpy.typing import NDArray
+from orbit_visualiser.core import Orbit, Satellite
+from tests.test_cases import standard_test_cases
+
+
+@pytest.mark.parametrize("e, rp, mu, closure, orbit_type", standard_test_cases)
+def test_specific_energy_conservation(
+    satellite_factory: Callable[[float, float, float], Satellite],
+    closed_anomaly_grid: NDArray[np.float64],
+    open_anomaly_grid: Callable[[Orbit, int], NDArray[np.float64]],
+    e: float,
+    rp: float,
+    mu: float,
+    closure: str,
+    orbit_type: str
+):
+    satellite: Satellite = satellite_factory(e, rp, mu)
+
+    if closure == "closed":
+        anomaly_grid = closed_anomaly_grid
+    else:
+        anomaly_grid = open_anomaly_grid(satellite._orbit)
+
+    for nu in anomaly_grid:
+        satellite.nu = nu
+        satellite.update_satellite_properties()
+
+        specific_energy = satellite.eps
+
+        specific_kin_energy = 0.5*satellite.v**2
+        specific_pot_energy = -mu/satellite.r
+        vis_viva_energy = specific_kin_energy + specific_pot_energy
+
+        assert np.isclose(specific_energy, vis_viva_energy)
+
+@pytest.mark.parametrize("e, rp, mu, closure, orbit_type", standard_test_cases)
+def test_specific_and_characteristic_energy_relation(
+    satellite_factory: Callable[[float, float, float], Satellite],
+    closed_anomaly_grid: NDArray[np.float64],
+    open_anomaly_grid: Callable[[Orbit, int], NDArray[np.float64]],
+    e: float,
+    rp: float,
+    mu: float,
+    closure: str,
+    orbit_type: str
+):
+    satellite: Satellite = satellite_factory(e, rp, mu)
+
+    if closure == "closed":
+        anomaly_grid = closed_anomaly_grid
+    else:
+        anomaly_grid = open_anomaly_grid(satellite._orbit)
+
+    for nu in anomaly_grid:
+        satellite.nu = nu
+        satellite.update_satellite_properties()
+
+        specific_energy = satellite.eps
+        characteristic_energy = satellite.c3
+
+        assert np.isclose(2*specific_energy, characteristic_energy)

tests/invariants/parabolic_orbit_test.py

@@ -0,0 +1 @@
+# TODO: Write test that gam = nu/2 for parabolic trajectory

tests/orbit_test.py

@@ -1,58 +1,25 @@
-from pytest import Subtests
+import pytest
 import numpy as np
+from typing import Callable
 from orbit_visualiser.core import Orbit
+from tests.test_cases import e_tagged_test_cases
 
+# TODO: implement sanity tests (use different formulae for the same value and check they are equal)
 
-
-def test_circular_orbit_distance_parameters(subtests: Subtests):
-    """
-    For circular orbits the radius of periapsis, radius of apoapsis, semi-major axis, semi-minor axis
-    and orbital parameter should all be equivalent.
-
-    circular_orbit_test_cases is the list of initial radius of periapsis we set for our orbits.
-    """
-    circular_orbit_test_cases = [
-        0.000000001,
-        1,
-        6789,
-        50000,
-        1000000,
-        4095384905805829034
-    ]
-
-    for i, periapsis in enumerate(circular_orbit_test_cases):
-        with subtests.test("Circular orbit distance parameters test cases", i = i):
-            orbit = Orbit(rp = periapsis)
-            periapsis = orbit.rp
-            distances = [
-                orbit.ra,
-                orbit.a,
-                orbit.b,
-                orbit.p
-            ]
-            assert np.allclose(distances, periapsis)
-
-def test_orbit_type(subtests: Subtests):
+@pytest.mark.parametrize("e, closure, orbit_type", e_tagged_test_cases)
+def test_orbit_type_sanity(
+    orbit_factory: Callable[[float, float], Orbit],
+    e: float,
+    closure: str,
+    orbit_type: str):
     """
     If the eccentricity is 0<= e < 1 then the orbit is closed (True) and either circular (e = 0) or
     elliptical (0 < e < 1). If the eccentricity is e >= 1 then the orbit is not closed (False) and
     is either parabolic (e = 1) or hyperbolic (e > 1).
-
-    eccentricity_test_cases is a dictionary giving the eccentricity test value and the expected
-    orbit type and is_closed boolean value.
     """
-    eccentricity_test_cases = {
-        0 : ("circular", True),
-        0.0000000001: ("elliptical", True),
-        0.5: ("elliptical", True),
-        0.9999999999: ("elliptical", True),
-        1: ("parabolic", False),
-        1.00000000001: ("hyperbolic", False),
-        1.5: ("hyperbolic", False),
-        10000000: ("hyperbolic", False)
-    }
+    orbit: Orbit = orbit_factory(e = e)
+
+    test_orbit_closure = "closed" if orbit.is_closed else "open"
+    test_orbit_type = orbit.orbit_type
 
-    for i, (e, output) in enumerate(eccentricity_test_cases.items()):
-        with subtests.test("Orbit type test cases", i = i):
-            orbit = Orbit(e = e)
-            assert orbit.orbit_type == output[0] and orbit.is_closed is output[1]
+    assert test_orbit_closure == closure and test_orbit_type == orbit_type