From 52f7cab09d0f4ba463b9b1778940b89afdeae070 Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 12 Aug 2025 15:50:00 +1000 Subject: [PATCH 1/2] Added a comprehensive test suite for KeplerianElements. Adds at least one test case for all methods and property getters and setters. Tests specifically address edge cases and cases that are handled differently by the underlying code (e.g. inclined vs non-inclined, circular vs non-circular orbits). Also introduces a new low-level case in test_utilities for the complicated function eccentric_anomaly_from_mean. Many of these cases are either expecting failure or commented out due to known bugs, which have been filed on GitHub and linked from the appropriate tests. The commented-out assertions demonstrate the correct expected behaviour once those bugs are fixed. --- tests/test_orbital.py | 1070 ++++++++++++++++++++++++++++++++++++++- tests/test_utilities.py | 27 + 2 files changed, 1092 insertions(+), 5 deletions(-) diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 702c2f6..ab76868 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -1,17 +1,410 @@ import unittest +import warnings +from math import sqrt, tau +import numpy as np +import numpy.testing +from astropy import time from numpy import radians +from numpy.linalg import norm from scipy.constants import kilo from orbital import KeplerianElements, earth, venus +from orbital.utilities import ConvergenceError, OrbitalWarning, Position, Velocity + +J2000 = time.Time("J2000", scale="utc") class TestOrbitalElements(unittest.TestCase): + def test_circular(self): + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + self.assertAlmostEqual(orbit.a, RADIUS) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + + numpy.testing.assert_almost_equal(orbit.r, Position(RADIUS, 0, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(0, sqrt(earth.mu / RADIUS), 0) + ) + + # Manually calculate angular velocity and period of a circular orbit. + self.assertAlmostEqual(orbit.n, sqrt(earth.mu / RADIUS**3)) + self.assertAlmostEqual(orbit.T, tau * sqrt(RADIUS**3 / earth.mu)) + self.assertAlmostEqual(orbit.fpa, 0.0) + + self.assertAlmostEqual(orbit.apocenter_radius, RADIUS) + self.assertAlmostEqual(orbit.pericenter_radius, RADIUS) + self.assertAlmostEqual(orbit.apocenter_altitude, RADIUS - earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_altitude, RADIUS - earth.mean_radius) + + numpy.testing.assert_almost_equal(orbit.U, np.array([1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + def test_circular_arg_pe(self): + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=radians(45), M0=0.0, body=earth + ) + self.assertAlmostEqual(orbit.arg_pe, radians(45)) + # The arg_pe gives the angle at the epoch, so these vectors should be + # rotated 45°. + numpy.testing.assert_almost_equal( + orbit.r, Position(RADIUS * 0.5 * sqrt(2), RADIUS * 0.5 * sqrt(2), 0) + ) + numpy.testing.assert_almost_equal( + orbit.v, + Velocity( + -sqrt(earth.mu / RADIUS) * 0.5 * sqrt(2), + sqrt(earth.mu / RADIUS) * 0.5 * sqrt(2), + 0, + ), + ) + + numpy.testing.assert_almost_equal( + orbit.U, np.array([0.5 * sqrt(2), 0.5 * sqrt(2), 0]) + ) + numpy.testing.assert_almost_equal( + orbit.V, np.array([-0.5 * sqrt(2), 0.5 * sqrt(2), 0]) + ) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + def test_anomaly_at_time(self): - orbit = KeplerianElements.with_period(90 * 60, e=0, body=earth) - self.assertAlmostEqual(orbit.M, 0) - orbit.t += orbit.T - self.assertAlmostEqual(orbit.M, 0) + RADIUS = 10000000.0 + orbit = KeplerianElements(a=RADIUS, M0=radians(90), body=earth) + # Test all of the properties that change with time. + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + # For a circular orbit, all three anomalies are the same. + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(90)) + self.assertAlmostEqual(orbit.f, radians(90)) + numpy.testing.assert_almost_equal(orbit.r, Position(0, RADIUS, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(-sqrt(earth.mu / RADIUS), 0, 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([-1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # 1/4 of the way around. + orbit.t = orbit.T * 0.25 + self.assertAlmostEqual( + orbit.epoch, J2000 + time.TimeDelta(orbit.T * 0.25, format="sec") + ) + self.assertAlmostEqual(orbit.M, radians(180)) + self.assertAlmostEqual(orbit.E, radians(180)) + self.assertAlmostEqual(orbit.f, radians(180)) + numpy.testing.assert_almost_equal(orbit.r, Position(-RADIUS, 0, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(0, -sqrt(earth.mu / RADIUS), 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([-1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, -1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # 1/2 of the way around. + orbit.t = orbit.T * 0.5 + self.assertAlmostEqual( + orbit.epoch, J2000 + time.TimeDelta(orbit.T * 0.5, format="sec") + ) + self.assertAlmostEqual(orbit.M, radians(270)) + self.assertAlmostEqual(orbit.E, radians(270)) + self.assertAlmostEqual(orbit.f, radians(270)) + numpy.testing.assert_almost_equal(orbit.r, Position(0, -RADIUS, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(sqrt(earth.mu / RADIUS), 0, 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([0, -1, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # A full revolution around. + orbit.t = orbit.T + self.assertAlmostEqual( + orbit.epoch, J2000 + time.TimeDelta(orbit.T, format="sec") + ) + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(90)) + self.assertAlmostEqual(orbit.f, radians(90)) + numpy.testing.assert_almost_equal(orbit.r, Position(0, RADIUS, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(-sqrt(earth.mu / RADIUS), 0, 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([-1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # 5/4 of the way around. + orbit.t = orbit.T * 1.25 + self.assertAlmostEqual( + orbit.epoch, J2000 + time.TimeDelta(orbit.T * 1.25, format="sec") + ) + self.assertAlmostEqual(orbit.M, radians(180)) + self.assertAlmostEqual(orbit.E, radians(180)) + self.assertAlmostEqual(orbit.f, radians(180)) + numpy.testing.assert_almost_equal(orbit.r, Position(-RADIUS, 0, 0)) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(0, -sqrt(earth.mu / RADIUS), 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([-1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, -1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + def test_zero(self): + orbit = KeplerianElements( + a=0.0, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + self.assertAlmostEqual(orbit.a, 0.0) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + + numpy.testing.assert_almost_equal(orbit.r, Position(0, 0, 0)) + # XXX Arguably, v should be the zero vector. + self.assertRaises(ZeroDivisionError, lambda: orbit.v) + + # XXX Arguably, n should be infinity. + self.assertRaises(ZeroDivisionError, lambda: orbit.n) + # XXX If n is infinity, T should be 0.0. + self.assertRaises(ZeroDivisionError, lambda: orbit.T) + self.assertAlmostEqual(orbit.fpa, 0.0) + + self.assertAlmostEqual(orbit.apocenter_radius, 0.0) + self.assertAlmostEqual(orbit.pericenter_radius, 0.0) + self.assertAlmostEqual(orbit.apocenter_altitude, -earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_altitude, -earth.mean_radius) + + numpy.testing.assert_almost_equal(orbit.U, np.array([1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + def test_inclined(self): + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, + e=0.0, + i=radians(45), + raan=radians(90), + arg_pe=0.0, + M0=0.0, + body=earth, + ) + self.assertAlmostEqual(orbit.a, RADIUS) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, radians(45)) + self.assertAlmostEqual(orbit.raan, radians(90)) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + + numpy.testing.assert_almost_equal(orbit.r, Position(0, RADIUS, 0)) + numpy.testing.assert_almost_equal( + orbit.v, + Velocity( + -sqrt(earth.mu / RADIUS) * 0.5 * sqrt(2), + 0, + sqrt(earth.mu / RADIUS) * 0.5 * sqrt(2), + ), + ) + + # n, T, apsides do not change when an inclination is introduced. + self.assertAlmostEqual(orbit.n, sqrt(earth.mu / RADIUS**3)) + self.assertAlmostEqual(orbit.T, tau * sqrt(RADIUS**3 / earth.mu)) + self.assertAlmostEqual(orbit.fpa, 0.0) + + self.assertAlmostEqual(orbit.apocenter_radius, RADIUS) + self.assertAlmostEqual(orbit.pericenter_radius, RADIUS) + self.assertAlmostEqual(orbit.apocenter_altitude, RADIUS - earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_altitude, RADIUS - earth.mean_radius) + + numpy.testing.assert_almost_equal(orbit.U, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal( + orbit.V, np.array([-0.5 * sqrt(2), 0, 0.5 * sqrt(2)]) + ) + numpy.testing.assert_almost_equal( + orbit.W, np.array([0.5 * sqrt(2), 0, 0.5 * sqrt(2)]) + ) + self.assertUVWMatches(orbit) + + # Advance time: 1/4 of the way around. + orbit.t = orbit.T * 0.25 + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(90)) + self.assertAlmostEqual(orbit.f, radians(90)) + numpy.testing.assert_almost_equal( + orbit.r, Position(-RADIUS * 0.5 * sqrt(2), 0, RADIUS * 0.5 * sqrt(2)) + ) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(0, -sqrt(earth.mu / RADIUS), 0) + ) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal( + orbit.U, np.array([-0.5 * sqrt(2), 0, 0.5 * sqrt(2)]) + ) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, -1, 0])) + numpy.testing.assert_almost_equal( + orbit.W, np.array([0.5 * sqrt(2), 0, 0.5 * sqrt(2)]) + ) + self.assertUVWMatches(orbit) + + def test_elliptical(self): + A = 10000000.0 + orbit = KeplerianElements( + a=A, e=0.75, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + self.assertAlmostEqual(orbit.a, A) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + + numpy.testing.assert_almost_equal(orbit.r, Position(2500000, 0, 0)) + numpy.testing.assert_almost_equal(orbit.v, Velocity(0, 16703.901013, 0)) + + # Angular velocity and period are the same as for a circular orbit. + self.assertAlmostEqual(orbit.n, sqrt(earth.mu / A**3)) + self.assertAlmostEqual(orbit.T, tau * sqrt(A**3 / earth.mu)) + self.assertAlmostEqual(orbit.fpa, 0.0) + + self.assertAlmostEqual(orbit.apocenter_radius, 17500000.0) + self.assertAlmostEqual(orbit.pericenter_radius, 2500000.0) + self.assertAlmostEqual(orbit.apocenter_altitude, 17500000.0 - earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_altitude, 2500000.0 - earth.mean_radius) + + numpy.testing.assert_almost_equal(orbit.U, np.array([1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + def test_elliptical_extreme_e(self): + A = 500.0 + e = 0.99999999 + orbit = KeplerianElements( + a=A, e=e, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + self.assertAlmostEqual(orbit.a, A) + self.assertAlmostEqual(orbit.e, e) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.epoch, J2000) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + + # Problem case: e is very close to 1.0 and M is very close to 360°. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.M = radians(359.9) + self.assertAlmostEqual(orbit.M, radians(359.9)) + # XXX This should produce the results below (it does if + # utilities.MAX_ITERATIONS is set to 100000), but instead it fails to + # converge in 100 iterations. + self.assertRaises(ConvergenceError, lambda: orbit.E) + # self.assertAlmostEqual(orbit.E, radians(347.454759)) + # self.assertAlmostEqual(orbit.f, radians(180.073718)) + + def test_elliptical_times(self): + A = 10000000.0 + orbit = KeplerianElements(a=A, e=0.75, body=earth) + # Test all of the properties that change with time. + # At periapsis. + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, 0.0) + self.assertAlmostEqual(orbit.E, 0.0) + self.assertAlmostEqual(orbit.f, 0.0) + numpy.testing.assert_almost_equal(orbit.r, Position(2500000, 0, 0)) + numpy.testing.assert_almost_equal(orbit.v, Velocity(0, 16703.901013, 0)) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, 1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # 1/4 of the period has elapsed (more than a quarter of the way around). + orbit.t = orbit.T * 0.25 + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(125.140095)) + self.assertAlmostEqual(orbit.f, radians(157.802569)) + numpy.testing.assert_almost_equal( + orbit.r, Position(-13255776.4031414, 5408888.899183, 0) + ) + numpy.testing.assert_almost_equal( + orbit.v, Velocity(-3606.1267047, -1678.8615886, 0) + ) + self.assertAlmostEqual(orbit.fpa, radians(42.837854)) + numpy.testing.assert_almost_equal(orbit.U, np.array([-0.9258875, 0.3777993, 0])) + numpy.testing.assert_almost_equal( + orbit.V, np.array([-0.3777993, -0.9258875, 0]) + ) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) + + # 1/2 of the period has elapsed, half way around (at apoapsis). + orbit.t = orbit.T * 0.5 + self.assertAlmostEqual(orbit.M, radians(180)) + self.assertAlmostEqual(orbit.E, radians(180)) + self.assertAlmostEqual(orbit.f, radians(180)) + numpy.testing.assert_almost_equal(orbit.r, Position(-17500000, 0, 0)) + numpy.testing.assert_almost_equal(orbit.v, Velocity(0, -2386.2715733, 0)) + self.assertAlmostEqual(orbit.fpa, 0.0) + numpy.testing.assert_almost_equal(orbit.U, np.array([-1, 0, 0])) + numpy.testing.assert_almost_equal(orbit.V, np.array([0, -1, 0])) + numpy.testing.assert_almost_equal(orbit.W, np.array([0, 0, 1])) + self.assertUVWMatches(orbit) def test_apsides(self): orbit = KeplerianElements.with_period(90 * 60, body=earth) @@ -52,10 +445,82 @@ def test_mean_motion(self): orbit.propagate_anomaly_to(M=radians(40)) self.assertAlmostEqual(orbit.M, radians(40)) - def test_period(self): + def test_with_altitude(self): + # TODO: This isn't the best example, as most of the cases below get a + # negative periapsis altitude because it's so close to the earth. + ALTITUDE = 10000.0 + + # Circular orbit. + orbit = KeplerianElements.with_altitude(ALTITUDE, body=earth) + self.assertAlmostEqual(norm(orbit.r), ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_radius, ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_radius, ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_altitude, ALTITUDE) + self.assertAlmostEqual(orbit.pericenter_altitude, ALTITUDE) + self.assertAlmostEqual(orbit.a, ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # Elliptical orbit. + orbit = KeplerianElements.with_altitude(ALTITUDE, e=0.75, body=earth) + self.assertAlmostEqual(norm(orbit.r), ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_radius, 38296000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_radius, ALTITUDE + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_altitude, 38296000.0) + self.assertAlmostEqual(orbit.pericenter_altitude, ALTITUDE) + self.assertAlmostEqual(orbit.a, 25524000.0) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + # Elliptical orbit, nonzero M0. + orbit = KeplerianElements.with_altitude( + ALTITUDE, e=0.75, M0=radians(35), body=earth + ) + self.assertAlmostEqual(orbit.M, radians(35)) + self.assertAlmostEqual(norm(orbit.r), ALTITUDE + earth.mean_radius) + self.assertAlmostEqual( + orbit.apocenter_radius, 7094311.533422537 + earth.mean_radius + ) + self.assertAlmostEqual( + orbit.pericenter_radius, -4447384.066653923 + earth.mean_radius + ) + self.assertAlmostEqual(orbit.apocenter_altitude, 7094311.533422537) + self.assertAlmostEqual(orbit.pericenter_altitude, -4447384.066653923) + self.assertAlmostEqual(orbit.a, 7694463.733384307) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, radians(35)) + self.assertAlmostEqual(orbit.t, 0.0) + + def test_with_period(self): orbit = KeplerianElements.with_period(2 * 60 * 60, M0=radians(35), body=earth) self.assertAlmostEqual(orbit.T, 2 * 60 * 60) + self.assertAlmostEqual(orbit.a, 8058997.3045416) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, radians(35)) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # TODO: These should be part of separate tests, as they have nothing to + # do with with_period. # Propagate to set t nonzero to test M0 fix when a is set. orbit.propagate_anomaly_by(M=radians(10)) @@ -69,6 +534,601 @@ def test_period(self): orbit.propagate_anomaly_to(M=radians(40)) self.assertAlmostEqual(orbit.M, radians(40)) + def test_with_apside_altitudes(self): + # Circular orbit. + orbit = KeplerianElements.with_apside_altitudes(10000.0, 10000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_radius, 10000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_radius, 10000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_altitude, 10000.0) + self.assertAlmostEqual(orbit.pericenter_altitude, 10000.0) + self.assertAlmostEqual(orbit.a, 10000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # Elliptical orbit. + orbit = KeplerianElements.with_apside_altitudes(10000.0, 38296000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_radius, 38296000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.pericenter_radius, 10000.0 + earth.mean_radius) + self.assertAlmostEqual(orbit.apocenter_altitude, 38296000.0) + self.assertAlmostEqual(orbit.pericenter_altitude, 10000.0) + self.assertAlmostEqual(orbit.a, 25524000.0) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + # Elliptical orbit, arguments reversed (should make no difference). + orbit = KeplerianElements.with_apside_altitudes(38296000.0, 10000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_altitude, 38296000.0) + self.assertAlmostEqual(orbit.pericenter_altitude, 10000.0) + self.assertAlmostEqual(orbit.a, 25524000.0) + self.assertAlmostEqual(orbit.e, 0.75) + + def test_with_apside_radii(self): + # Circular orbit. + orbit = KeplerianElements.with_apside_radii(10000000.0, 10000000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_radius, 10000000.0) + self.assertAlmostEqual(orbit.pericenter_radius, 10000000.0) + self.assertAlmostEqual(orbit.apocenter_altitude, 10000000.0 - earth.mean_radius) + self.assertAlmostEqual( + orbit.pericenter_altitude, 10000000.0 - earth.mean_radius + ) + self.assertAlmostEqual(orbit.a, 10000000.0) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # Elliptical orbit. + orbit = KeplerianElements.with_apside_radii(10000000.0, 20000000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_radius, 20000000.0) + self.assertAlmostEqual(orbit.pericenter_radius, 10000000.0) + self.assertAlmostEqual(orbit.apocenter_altitude, 20000000.0 - earth.mean_radius) + self.assertAlmostEqual( + orbit.pericenter_altitude, 10000000.0 - earth.mean_radius + ) + self.assertAlmostEqual(orbit.a, 15000000.0) + self.assertAlmostEqual(orbit.e, 1.0 / 3.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + # Elliptical orbit, arguments reversed (should make no difference). + orbit = KeplerianElements.with_apside_radii(20000000.0, 10000000.0, body=earth) + self.assertAlmostEqual(orbit.apocenter_radius, 20000000.0) + self.assertAlmostEqual(orbit.pericenter_radius, 10000000.0) + self.assertAlmostEqual(orbit.a, 15000000.0) + self.assertAlmostEqual(orbit.e, 1.0 / 3.0) + + def test_from_state_vector_circular(self): + # Circular orbit. + RADIUS = 10000000.0 + R = Position(RADIUS, 0, 0) + V = Velocity(0, sqrt(earth.mu / RADIUS), 0) + + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, RADIUS) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # 1/4 of the way around. + R = Position(0, RADIUS, 0) + V = Velocity(-sqrt(earth.mu / RADIUS), 0, 0) + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, RADIUS) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + # t is always 0; instead it is expected to set M0. + self.assertAlmostEqual(orbit.M0, radians(90)) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, radians(90)) + + def test_from_state_vector_zero(self): + # Degenerate orbit with r=0 and v=0. + R = Position(0, 0, 0) + V = Velocity(0, 0, 0) + + # This case currently violates an assertion due to internal values being + # NaN. + # XXX An AssertionError indicates an unexpected case. This should + # probably be explicitly detected and raise another kind of error. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + self.assertRaises( + AssertionError, KeplerianElements.from_state_vector, R, V, body=earth + ) + + def test_from_state_vector_circular_retrograde(self): + # Circular orbit, retrograde. + # Regression test for https://github.com/RazerM/orbital/issues/18. + RADIUS = 10000000.0 + R = Position(RADIUS, 0, 0) + V = Velocity(0, -sqrt(earth.mu / RADIUS), 0) + + # XXX This erroneously generates warnings and raises an assertion, due + # to internal values being NaN. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + self.assertRaises( + AssertionError, KeplerianElements.from_state_vector, R, V, body=earth + ) + # The expected values, after this bug is fixed: + # numpy.testing.assert_almost_equal(orbit.r, R) + # numpy.testing.assert_almost_equal(orbit.v, V) + # self.assertAlmostEqual(orbit.a, RADIUS) + # self.assertAlmostEqual(orbit.e, 0.0) + # self.assertAlmostEqual(orbit.i, radians(180)) + # self.assertAlmostEqual(orbit.raan, 0.0) + # self.assertAlmostEqual(orbit.arg_pe, 0.0) + # self.assertAlmostEqual(orbit.M0, 0.0) + + # self.assertAlmostEqual(orbit.ref_epoch, J2000) + # self.assertEqual(orbit.body, earth) + # self.assertAlmostEqual(orbit.t, 0.0) + + def test_from_state_vector_inclined(self): + # Inclined circular orbit, 1/4 of the way around. + # Regression test for https://github.com/RazerM/orbital/issues/38. + RADIUS = 10000000.0 + R = Position(-RADIUS * 0.5 * sqrt(2), 0, RADIUS * 0.5 * sqrt(2)) + V = Velocity(0, -sqrt(earth.mu / RADIUS), 0) + + with warnings.catch_warnings(): + # XXX This has a warning for dividing by zero. + warnings.simplefilter("ignore", category=RuntimeWarning) + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + # XXX: r, v and arg_pe are nan due to a bug in from_state_vector. + # This happens for a perfect circle, but doesn't happen in + # test_from_state_vector_circular for some reason. + # numpy.testing.assert_almost_equal(orbit.r, R) + # numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, RADIUS) + self.assertAlmostEqual(orbit.e, 0.0) + self.assertAlmostEqual(orbit.i, radians(45)) + self.assertAlmostEqual(orbit.raan, radians(90)) + # self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, radians(90)) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + def test_from_state_vector_elliptical(self): + # Elliptical orbit, at periapsis. + R = Position(2500000, 0, 0) + V = Velocity(0, 16703.9010129, 0) + + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, 10000000.0, places=3) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + # Elliptical orbit, 1/4 of the period elapsed. + R = Position(-13255776.4031414, 5408888.899183, 0) + V = Velocity(-3606.1267047, -1678.8615886, 0) + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + numpy.testing.assert_almost_equal(orbit.r, R, decimal=4) + numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, 10000000.0, places=3) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, radians(90)) + self.assertAlmostEqual(orbit.t, 0.0) + self.assertAlmostEqual(orbit.M, radians(90)) + + def test_from_state_vector_elliptical_arg_pe_gt_180(self): + # Elliptical orbit, at periapsis, with arg_pe > 180°. + # Regression test for https://github.com/RazerM/orbital/issues/39. + R = Position(0, -2500000, 0) + V = Velocity(16703.9010129, 0, 0) + + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + # XXX These do not match (they are 180° out, due to arg_pe). + # numpy.testing.assert_almost_equal(orbit.r, R) + # numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.a, 10000000.0, places=3) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + # XXX This is incorrectly calculated as 90°. + # self.assertAlmostEqual(orbit.arg_pe, radians(270.0)) + self.assertAlmostEqual(orbit.M0, 0.0) + + self.assertAlmostEqual(orbit.ref_epoch, J2000) + self.assertEqual(orbit.body, earth) + self.assertAlmostEqual(orbit.t, 0.0) + + def test_from_state_vector_f_at_periapsis(self): + # Elliptical orbit, inclined, at periapsis. + # Regression test for https://github.com/RazerM/orbital/issues/40. + # In this particular case, the above bug would cause f to be nan due to + # floating point rounding errors. + R = Position(0, -1767766.952966369, -1767766.952966369) + V = Velocity(16703.901013, 0, 0) + + # XXX Currently crashes due to the above bug. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + self.assertRaises( + AssertionError, KeplerianElements.from_state_vector, R, V, body=earth + ) + # numpy.testing.assert_almost_equal(orbit.r, R) + # numpy.testing.assert_almost_equal(orbit.v, V) + # self.assertAlmostEqual(orbit.a, 10000000.0, places=2) + # self.assertAlmostEqual(orbit.e, 0.75) + # self.assertAlmostEqual(orbit.i, radians(45)) + # self.assertAlmostEqual(orbit.raan, 0.0) + # self.assertAlmostEqual(orbit.arg_pe, radians(270.0)) + # self.assertAlmostEqual(orbit.M0, 0.0) + + # self.assertAlmostEqual(orbit.ref_epoch, J2000) + # self.assertEqual(orbit.body, earth) + # self.assertAlmostEqual(orbit.t, 0.0) + + def test_from_state_vector_iss(self): + # ISS (Zarya) from 2008-09-20 12:25:40 + # Source: SGP4 parsed TLE example from + # https://en.wikipedia.org/wiki/Two-line_element_set + # (This is the same state vector used by test_from_tle below.) + # Used as a real-world smoke test, just ensure that the vectors + # round-trip. + R = Position(4083902.4635207, -993631.9996058, 5243603.6653708) + V = Velocity(2512.8372952, 7259.888525, -583.7785365) + orbit = KeplerianElements.from_state_vector(R, V, body=earth) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + + def test_from_state_vector_roundtrips(self): + # Simple round-trip tests for a lot of edge cases in + # utilities.elements_from_state_vector. Just confirm that r and v come + # back relatively unchanged for all of these inputs (without + # hand-verifying all the elements). + + # List of (r, v) pairs. + # XXX Commented-out cases are failing. + CASES = [ + # Zero. + # (Position(0, 0, 0), Velocity(0, 0, 0)), + # Circular flat, prograde. + (Position(10000000, 0, 0), Velocity(0, 6313.4811435530555, 0)), + # Circular flat, f > 180°. + (Position(0, -10000000, 0), Velocity(6313.4811435530555, 0, 0)), + # Circular flat, retrograde. + # (Position(10000000, 0, 0), Velocity(0, -6313.4811435530555, 0)), + # Circular inclined, prograde. + # (Position(10000000, 0, 0), Velocity(0, 4464.305329499764, 4464.305329499764)), + # Circular inclined, prograde (raan 90°, M0 90°). + # (Position(-7071067.811865476, 0, 7071067.811865476), Velocity(0, -6313.4811435530555, 0)), + # Circular inclined, raan > 180°. + # (Position(0, -10000000, 0), Velocity(4464.305329499764, 0, 4464.305329499764)), + # Circular inclined, f > 180°. + # (Position(0, -7071067.811865476, -7071067.811865476), Velocity(6313.4811435530555, 0, 0)), + # Circular polar. + # (Position(10000000, 0, 0), Velocity(0, 0, 6313.4811435530555)), + # Elliptical flat (e=0.75), prograde. + (Position(2500000, 0, 0), Velocity(0, 16703.901013, 0)), + # Elliptical flat (e=0.75), arg_pe > 180°. + # (Position(0, -2500000, 0), Velocity(16703.901013, 0, 0)), + # Elliptical flat (e=0.75), retrograde. + # (Position(2500000, 0, 0), Velocity(0, -16703.901013, 0)), + # Elliptical flat (e=0.75), f > 180°. + ( + Position(-13255776.4031414, -5408888.899183, 0), + Velocity(3606.1267047, -1678.8615886, 0), + ), + # Elliptical inclined. + ( + Position(2500000, 0, 0), + Velocity(0, 11811.441678561141, 11811.441678561141), + ), + # Elliptical inclined, arg_pe > 180°. + # (Position(0, -1767766.952966369, -1767766.952966369), Velocity(16703.901013, 0, 0)), + ] + + for i, (r, v) in enumerate(CASES): + orbit = KeplerianElements.from_state_vector(r, v, body=earth) + numpy.testing.assert_almost_equal( + orbit.r, r, decimal=4, err_msg=f"Case #{i:d}: {orbit}" + ) + numpy.testing.assert_almost_equal( + orbit.v, v, decimal=4, err_msg=f"Case #{i:d}: {orbit}" + ) + + def test_from_tle(self): + # Sample TLE from Wikipedia: + # https://en.wikipedia.org/wiki/Two-line_element_set + # ISS (ZARYA) + LINE1 = "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927" + LINE2 = "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537" + orbit = KeplerianElements.from_tle(LINE1, LINE2, body=earth) + self.assertAlmostEqual(orbit.t, 0) + self.assertTimeEqual(orbit.ref_epoch, time.Time("2008-09-20 12:25:40.0")) + + # These values for r and v were computed by SGP4. Internally, from_tle + # uses these state vectors as an intermediate step, so verify that the + # resulting r and v values match. + R = Position(4083902.4635207, -993631.9996058, 5243603.6653708) + V = Velocity(2512.8372952, 7259.888525, -583.7785365) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + + # NOTE: These expected values have just been set to match the output + # from this function. Note that the e, i, raan, arg_pe and M0 are + # given directly in the TLE data, and a can be computed from the mean + # motion in the TLE data as shown here. + # EXPECTED_N_REV_PER_DAY = 15.72125391 + # EXPECTED_N = EXPECTED_N_REV_PER_DAY * tau / 86400 # [rad/s] + # EXPECTED_A = (earth.mu / EXPECTED_N ** 2) ** (1 / 3) + # + # Yet the output values do not match. This is likely due to rounding + # errors converting to state vectors and back. The largest discrepancy + # is a, which is off by 5 km. See the note about arg_pe and M0. + # + # The comment after each line gives the actual expected value based on + # the TLE data. + self.assertAlmostEqual(orbit.a, 6725547.816501163) # 6730960.68?? + self.assertAlmostEqual(orbit.e, 0.0008330) # 0.0006703?? + self.assertAlmostEqual(orbit.i, radians(51.621653)) # 51.6416?? + self.assertAlmostEqual(orbit.raan, radians(247.45773)) # 247.4627?? + # Note: The following two values are off by a huge amount (about 18°) + # but the errors cancel out. Given the eccentricity is so low, the sum + # of these two angles is all that really matters. + self.assertAlmostEqual(orbit.arg_pe, radians(112.50348)) # 130.5360?? + self.assertAlmostEqual(orbit.M0, radians(343.056983)) # 325.0288?? + + def test_set_epoch(self): + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + self.assertTimeEqual(orbit.epoch, time.Time("2000-01-01 12:00:00")) + orbit.epoch = time.Time("2000-01-01 12:03:00") + self.assertTimeEqual(orbit.epoch, time.Time("2000-01-01 12:03:00")) + self.assertAlmostEqual(orbit.t, 180.0) + self.assertAlmostEqual(orbit.M, (180 * orbit.n) % tau) + + def test_set_M(self): + # Circular trajectory. + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.M = radians(495) + self.assertAlmostEqual(orbit.M, radians(135)) + self.assertAlmostEqual(orbit.E, radians(135)) + self.assertAlmostEqual(orbit.f, radians(135)) + # t does not get set. + self.assertAlmostEqual(orbit.t, 0.0) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.M = radians(-45) + self.assertAlmostEqual(orbit.M, radians(315)) + self.assertAlmostEqual(orbit.E, radians(315)) + self.assertAlmostEqual(orbit.f, radians(315)) + + def test_set_E(self): + # Elliptical trajectory. + A = 10000000.0 + orbit = KeplerianElements( + a=A, e=0.75, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + # Set E such that M is 90 degrees. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.E = radians(125.140095) + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(125.140095)) + self.assertAlmostEqual(orbit.f, radians(157.802569)) + # t does not get set. + self.assertAlmostEqual(orbit.t, 0.0) + + # Test values outside the range (0, tau]. + # Regression test for https://github.com/RazerM/orbital/issues/37. + # XXX Unlike setting M or f, setting E does not mod by tau. Attempting + # to read the value of E back results in non-convergence because the + # value of M is outside of the expected range. + # Commented-out asserts are failing. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.E = radians(485.140095) + # self.assertAlmostEqual(orbit.M, radians(90)) + # self.assertAlmostEqual(orbit.E, radians(125.140095)) + # self.assertAlmostEqual(orbit.f, radians(157.802569)) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.E = radians(-180) + # self.assertAlmostEqual(orbit.M, radians(180)) + # self.assertAlmostEqual(orbit.E, radians(180)) + # self.assertAlmostEqual(orbit.f, radians(180)) + + def test_set_f(self): + # Elliptical trajectory. + A = 10000000.0 + orbit = KeplerianElements( + a=A, e=0.75, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + # Set f such that M is 90 degrees. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.f = radians(157.802569) + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(125.140095)) + self.assertAlmostEqual(orbit.f, radians(157.802569)) + # t does not get set. + self.assertAlmostEqual(orbit.t, 0.0) + + # mod tau test. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.f = radians(517.802569) + self.assertAlmostEqual(orbit.M, radians(90)) + self.assertAlmostEqual(orbit.E, radians(125.140095)) + self.assertAlmostEqual(orbit.f, radians(157.802569)) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=OrbitalWarning) + orbit.f = radians(-180) + self.assertAlmostEqual(orbit.M, radians(180)) + self.assertAlmostEqual(orbit.E, radians(180)) + self.assertAlmostEqual(orbit.f, radians(180)) + + def test_set_a(self): + orbit = KeplerianElements( + a=10000000.0, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=radians(45), body=earth + ) + orbit.t = orbit.T * 0.125 # Another 45 degrees. + self.assertAlmostEqual(orbit.M, radians(90)) + orbit.a = 5000000.0 + self.assertAlmostEqual(orbit.a, 5000000.0) + # M should be preserved. + self.assertAlmostEqual(orbit.M, radians(90)) + + def test_set_v(self): + # Start with a circular trajectory. + RADIUS = 2500000.0 + R = Position(RADIUS, 0, 0) + V = Velocity(0, sqrt(earth.mu / RADIUS), 0) + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + + # Elliptical trajectory. + # Increase the velocity at the periapsis, which should grow the apoapsis + # and make an elliptical orbit. The r and periapsis should remain the + # same. + V = Velocity(0, 16703.901013, 0) + orbit.v = V + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, V) + self.assertAlmostEqual(orbit.pericenter_radius, RADIUS) + self.assertAlmostEqual(orbit.a, 10000000, places=2) + self.assertAlmostEqual(orbit.e, 0.75) + self.assertAlmostEqual(orbit.i, 0.0) + self.assertAlmostEqual(orbit.raan, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.M0, 0.0) + + def test_set_v_arg_pe_gt_180(self): + # Test the same concept as + # test_from_state_vector_elliptical_arg_pe_gt_180 but with the v setter + # instead of from_state_vector. + # Regression test for https://github.com/RazerM/orbital/issues/39. + # Elliptical orbit, at periapsis, with arg_pe = 90°. + orbit = KeplerianElements( + a=10000000.0, + e=0.75, + i=0.0, + raan=0.0, + arg_pe=radians(90), + M0=0.0, + body=earth, + ) + R = Position(0, 2500000, 0) + numpy.testing.assert_almost_equal(orbit.r, R) + numpy.testing.assert_almost_equal(orbit.v, Velocity(-16703.9010129, 0, 0)) + + # Reduce velocity to switch around the apoapsis and periapsis, so + # arg_pe = 270°. + V = Velocity(-10000, 0, 0) + + def set_v(value): + orbit.v = value + + # XXX The 'r and v changed' detection logic is triggered in this case, + # causing a RuntimeError to be raised. If this was not raised, the + # following asserts would be wildly off. + self.assertRaises(RuntimeError, set_v, V) + # numpy.testing.assert_almost_equal(orbit.r, R) + # numpy.testing.assert_almost_equal(orbit.v, V) + # arg_pe should have rotated around 180°, and M0 to match (so r is in + # the same spot as it was before). + # XXX This is incorrectly calculated as 90°. + # self.assertAlmostEqual(orbit.arg_pe, radians(270.0)) + self.assertAlmostEqual(orbit.M0, radians(180.0)) + + def test_set_n(self): + # Circular trajectory. + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + NEW_N = sqrt(earth.mu / 5000000.0**3) + orbit.n = NEW_N + self.assertAlmostEqual(orbit.n, NEW_N) + self.assertAlmostEqual(orbit.T, tau / NEW_N) + self.assertAlmostEqual(orbit.a, 5000000.0) + + def test_set_T(self): + # Circular trajectory. + RADIUS = 10000000.0 + orbit = KeplerianElements( + a=RADIUS, e=0.0, i=0.0, raan=0.0, arg_pe=0.0, M0=0.0, body=earth + ) + NEW_T = tau * sqrt(5000000.0**3 / earth.mu) + orbit.T = NEW_T + self.assertAlmostEqual(orbit.n, tau / NEW_T) + self.assertAlmostEqual(orbit.T, NEW_T) + self.assertAlmostEqual(orbit.a, 5000000.0) + + def assertUVWMatches(self, orbit): + """Check that orbit's UVW matches U, V and W. + + This should always be true. + """ + numpy.testing.assert_almost_equal(orbit.UVW[0], orbit.U) + numpy.testing.assert_almost_equal(orbit.UVW[1], orbit.V) + numpy.testing.assert_almost_equal(orbit.UVW[2], orbit.W) + + def assertTimeEqual(self, t1, t2): + """Assert that two Time values match to the nearest second.""" + self.assertEqual( + t1.strftime("%Y-%b-%d %H:%M:%S"), t2.strftime("%Y-%b-%d %H:%M:%S") + ) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 6a0773c..1514fcc 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -1,6 +1,8 @@ import math import unittest +from math import fmod, radians, tau +from orbital import utilities from orbital.utilities import mod @@ -25,3 +27,28 @@ def test_mod(self): self.assertEqual(mod(-3.0, ninf), -3.0) self.assertEqual(mod(0.0, 3.0), 0.0) self.assertEqual(mod(0.0, ninf), 0.0) + + def test_eccentric_anomaly_from_mean(self): + # Pairs of (e, E) to test round-trips. + CASES = [ + (0.0, 0.0), + (0.0, radians(90)), + (0.0, radians(270)), + (0.75, 0.0), + (0.75, radians(90)), + (0.75, radians(270)), + (0.75, radians(450)), + (0.75, radians(-90)), + (0.75, 0.000001), + (1.0, radians(90)), + (1.0, radians(270)), + ] + # These functions should round-trip, no matter what value is given (as + # long as e <= 1.0). + for i, (e, E) in enumerate(CASES): + M = utilities.mean_anomaly_from_eccentric(e, E) + self.assertAlmostEqual( + utilities.eccentric_anomaly_from_mean(e, M), + fmod(E, tau), + msg=f"Case #{i:d}: ({e:f}, {E:f})", + ) From df971d0d0c802ef29ae8f6437fdeed4f3bc32043 Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 12 Aug 2025 16:33:33 +1000 Subject: [PATCH 2/2] Explain that test_from_tle is a regression test for #43. --- tests/test_orbital.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/tests/test_orbital.py b/tests/test_orbital.py index ab76868..083d6ab 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -890,6 +890,8 @@ def test_from_tle(self): numpy.testing.assert_almost_equal(orbit.r, R) numpy.testing.assert_almost_equal(orbit.v, V) + # Regression test for https://github.com/RazerM/orbital/issues/43. + # NOTE: These expected values have just been set to match the output # from this function. Note that the e, i, raan, arg_pe and M0 are # given directly in the TLE data, and a can be computed from the mean @@ -898,21 +900,21 @@ def test_from_tle(self): # EXPECTED_N = EXPECTED_N_REV_PER_DAY * tau / 86400 # [rad/s] # EXPECTED_A = (earth.mu / EXPECTED_N ** 2) ** (1 / 3) # - # Yet the output values do not match. This is likely due to rounding - # errors converting to state vectors and back. The largest discrepancy - # is a, which is off by 5 km. See the note about arg_pe and M0. + # Yet the output values do not match. This is due to rounding errors + # converting to state vectors and back. The largest discrepancy is a, + # which is off by 5 km. See the note about arg_pe and M0. # # The comment after each line gives the actual expected value based on # the TLE data. - self.assertAlmostEqual(orbit.a, 6725547.816501163) # 6730960.68?? - self.assertAlmostEqual(orbit.e, 0.0008330) # 0.0006703?? - self.assertAlmostEqual(orbit.i, radians(51.621653)) # 51.6416?? - self.assertAlmostEqual(orbit.raan, radians(247.45773)) # 247.4627?? + self.assertAlmostEqual(orbit.a, 6725547.816501163) # 6730960.68 + self.assertAlmostEqual(orbit.e, 0.0008330) # 0.0006703 + self.assertAlmostEqual(orbit.i, radians(51.621653)) # 51.6416 + self.assertAlmostEqual(orbit.raan, radians(247.45773)) # 247.4627 # Note: The following two values are off by a huge amount (about 18°) # but the errors cancel out. Given the eccentricity is so low, the sum # of these two angles is all that really matters. - self.assertAlmostEqual(orbit.arg_pe, radians(112.50348)) # 130.5360?? - self.assertAlmostEqual(orbit.M0, radians(343.056983)) # 325.0288?? + self.assertAlmostEqual(orbit.arg_pe, radians(112.50348)) # 130.5360 + self.assertAlmostEqual(orbit.M0, radians(343.056983)) # 325.0288 def test_set_epoch(self): RADIUS = 10000000.0