From 2eaa4442cc07ee6f73f2633f0b8e27db2d5e67df Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 16 Sep 2025 17:49:18 +1000 Subject: [PATCH 1/4] =?UTF-8?q?from=5Fstate=5Fvector:=20Fixed=20arg=5Fpe?= =?UTF-8?q?=20incorrect=20when=20it=20should=20be=20>=20180=C2=B0.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a missing check for the case where the eccentricity vector is in the third or fourth quadrant, requiring arg_pe to be negated. Note that this causes an edge case where arg_pe should be 0° but it comes out as 360° because we are not using modulo, but rather adding 2*pi. I didn't want to fix that in this commit because all the calculations in from_state_vector behave the same way. Fixes #39. --- src/orbital/utilities.py | 2 ++ tests/test_orbital.py | 29 +++++++++++------------------ 2 files changed, 13 insertions(+), 18 deletions(-) diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py index d6a9469..0d17802 100644 --- a/src/orbital/utilities.py +++ b/src/orbital/utilities.py @@ -284,6 +284,8 @@ def elements_from_state_vector(r, v, mu): # Argument of periapsis is the angle between # eccentricity vector and its x component. arg_pe = acos(ev.x / norm(ev)) + if ev.y < 0: + arg_pe = 2 * pi - arg_pe else: # Right ascension of ascending node is the angle # between the node vector and its x component. diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 083d6ab..3e935d0 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -749,7 +749,9 @@ def test_from_state_vector_elliptical(self): 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) + # XXX This should really be 0, but from_state_vector does not properly + # mod 2*pi, so it sometimes returns exactly 2*pi. + self.assertAlmostEqual(orbit.arg_pe, radians(360)) self.assertAlmostEqual(orbit.M0, radians(90)) self.assertAlmostEqual(orbit.t, 0.0) self.assertAlmostEqual(orbit.M, radians(90)) @@ -761,15 +763,13 @@ def test_from_state_vector_elliptical_arg_pe_gt_180(self): 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) + 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.arg_pe, radians(270.0)) self.assertAlmostEqual(orbit.M0, 0.0) self.assertAlmostEqual(orbit.ref_epoch, J2000) @@ -846,7 +846,7 @@ def test_from_state_vector_roundtrips(self): # 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)), + (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°. @@ -1077,19 +1077,12 @@ def test_set_v_arg_pe_gt_180(self): # 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) + orbit.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.arg_pe, radians(270.0)) self.assertAlmostEqual(orbit.M0, radians(180.0)) def test_set_n(self): From b52ac6cca23cee462c5dd60567e42a337ee0d36a Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 16 Sep 2025 18:04:08 +1000 Subject: [PATCH 2/4] from_state_vector: Fix some angular values returning 2*pi instead of 0. When negating an angle, consistently use mod 2*pi instead of adding 2*pi. This causes 0 to remain as 0 instead of becoming 2*pi. --- src/orbital/utilities.py | 12 ++++++------ tests/test_orbital.py | 4 +--- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py index 0d17802..2f24710 100644 --- a/src/orbital/utilities.py +++ b/src/orbital/utilities.py @@ -285,13 +285,13 @@ def elements_from_state_vector(r, v, mu): # eccentricity vector and its x component. arg_pe = acos(ev.x / norm(ev)) if ev.y < 0: - arg_pe = 2 * pi - arg_pe + arg_pe = mod(-arg_pe, 2 * pi) else: # Right ascension of ascending node is the angle # between the node vector and its x component. raan = acos(n.x / norm(n)) if n.y < 0: - raan = 2 * pi - raan + raan = mod(-raan, 2 * pi) # Argument of periapsis is angle between # node and eccentricity vectors. @@ -303,23 +303,23 @@ def elements_from_state_vector(r, v, mu): # vector and its x component. f = acos(r.x / norm(r)) if v.x > 0: - f = 2 * pi - f + f = mod(-f, 2 * pi) else: # True anomaly is angle between node # vector and position vector. f = acos(dot(n, r) / (norm(n) * norm(r))) if dot(n, v) > 0: - f = 2 * pi - f + f = mod(-f, 2 * pi) else: if ev.z < 0: - arg_pe = 2 * pi - arg_pe + arg_pe = mod(-arg_pe, 2 * pi) # True anomaly is angle between eccentricity # vector and position vector. f = acos(dot(ev, r) / (norm(ev) * norm(r))) if dot(r, v) < 0: - f = 2 * pi - f + f = mod(-f, 2 * pi) return OrbitalElements(a=a, e=e, i=i, raan=raan, arg_pe=arg_pe, f=f) diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 3e935d0..659d44a 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -749,9 +749,7 @@ def test_from_state_vector_elliptical(self): self.assertAlmostEqual(orbit.e, 0.75) self.assertAlmostEqual(orbit.i, 0.0) self.assertAlmostEqual(orbit.raan, 0.0) - # XXX This should really be 0, but from_state_vector does not properly - # mod 2*pi, so it sometimes returns exactly 2*pi. - self.assertAlmostEqual(orbit.arg_pe, radians(360)) + 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)) From d33869071b60de8831e9aee49c2d48da2f3e9288 Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 16 Sep 2025 18:15:38 +1000 Subject: [PATCH 3/4] from_state_vector: Fixed circular inclined orbits having arg_pe = nan. This was caused by missing check for the case where e ~= 0 when orbits are inclined, thus dividing by the eccentricity resulted in a divide-by-zero. Fixes #38. --- src/orbital/utilities.py | 11 ++++++++--- tests/test_orbital.py | 36 +++++++++++++++++++++--------------- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py index 2f24710..524c5f5 100644 --- a/src/orbital/utilities.py +++ b/src/orbital/utilities.py @@ -293,9 +293,14 @@ def elements_from_state_vector(r, v, mu): if n.y < 0: raan = mod(-raan, 2 * pi) - # Argument of periapsis is angle between - # node and eccentricity vectors. - arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev))) + if abs(e) < SMALL_NUMBER: + # For circular orbits, place periapsis + # at ascending node by convention + arg_pe = 0 + else: + # Argument of periapsis is angle between + # node and eccentricity vectors. + arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev))) if abs(e) < SMALL_NUMBER: if abs(i) < SMALL_NUMBER: diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 659d44a..9cfbd6f 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -700,20 +700,14 @@ def test_from_state_vector_inclined(self): 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) + 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, radians(45)) self.assertAlmostEqual(orbit.raan, radians(90)) - # self.assertAlmostEqual(orbit.arg_pe, 0.0) + self.assertAlmostEqual(orbit.arg_pe, 0.0) self.assertAlmostEqual(orbit.M0, radians(90)) self.assertAlmostEqual(orbit.ref_epoch, J2000) @@ -832,15 +826,27 @@ def test_from_state_vector_roundtrips(self): # 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)), + ( + 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)), + ( + 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)), + ( + 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)), + ( + Position(0, -7071067.811865476, -7071067.811865476), + Velocity(6313.4811435530555, 0, 0), + ), # Circular polar. - # (Position(10000000, 0, 0), Velocity(0, 0, 6313.4811435530555)), + (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°. From 0b374733a7c4e6b77e6d03d7b2ebb005782a454b Mon Sep 17 00:00:00 2001 From: Matt Giuca Date: Tue, 16 Sep 2025 18:21:22 +1000 Subject: [PATCH 4/4] from_state_vector: Fixed assertion error for retrograde non-inclined orbits. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This was calculating both raan and f as nan, since i=180° was going into the inclined orbit code path, when it is actually non-inclined. When invoked via KeplerianElements.from_state_vector, the nan value would fail an assertion. This also adds an important case needed to get the correct value of arg_pe for retrograde non-inclined orbits: if i=180°, arg_pe needs to be inverted since the orbit itself is upside down. Added two regression test cases for this case. Fixes #18 and #25. --- src/orbital/utilities.py | 11 ++++++++--- tests/test_orbital.py | 39 ++++++++++++++++++--------------------- 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py index 524c5f5..baaac86 100644 --- a/src/orbital/utilities.py +++ b/src/orbital/utilities.py @@ -271,8 +271,9 @@ def elements_from_state_vector(r, v, mu): # Inclination is the angle between the angular # momentum vector and its z component. i = acos(h.z / norm(h)) + i_small = abs(i) < SMALL_NUMBER or abs(i - pi) < SMALL_NUMBER - if abs(i) < SMALL_NUMBER: + if i_small: # For non-inclined orbits, raan is undefined; # set to zero by convention raan = 0 @@ -285,7 +286,11 @@ def elements_from_state_vector(r, v, mu): # eccentricity vector and its x component. arg_pe = acos(ev.x / norm(ev)) if ev.y < 0: - arg_pe = mod(-arg_pe, 2 * pi) + arg_pe = -arg_pe + # If plane is upside down, arg_pe needs to be inverted. + if abs(i) > pi * 0.5: + arg_pe = -arg_pe + arg_pe = mod(arg_pe, 2 * pi) else: # Right ascension of ascending node is the angle # between the node vector and its x component. @@ -303,7 +308,7 @@ def elements_from_state_vector(r, v, mu): arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev))) if abs(e) < SMALL_NUMBER: - if abs(i) < SMALL_NUMBER: + if i_small: # True anomaly is angle between position # vector and its x component. f = acos(r.x / norm(r)) diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 9cfbd6f..3881e95 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -672,26 +672,19 @@ def test_from_state_vector_circular_retrograde(self): 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) + 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, 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) + 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. @@ -824,7 +817,7 @@ def test_from_state_vector_roundtrips(self): # 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)), + (Position(10000000, 0, 0), Velocity(0, -6313.4811435530555, 0)), # Circular inclined, prograde. ( Position(10000000, 0, 0), @@ -852,7 +845,11 @@ def test_from_state_vector_roundtrips(self): # 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)), + (Position(2500000, 0, 0), Velocity(0, -16703.901013, 0)), + # Elliptical flat (e=0.75), retrograde, arg_pe = 90°, at periapsis. + (Position(0, -2500000, 0), Velocity(-16703.901013, 0, 0)), + # Elliptical flat, retrograde, arg_pe = 90°, at apoapsis. + (Position(0, -2500000, 0), Velocity(-10000.0, 0, 0)), # Elliptical flat (e=0.75), f > 180°. ( Position(-13255776.4031414, -5408888.899183, 0),