From 2eaa4442cc07ee6f73f2633f0b8e27db2d5e67df Mon Sep 17 00:00:00 2001
From: Matt Giuca <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 17:49:18 +1000
Subject: [PATCH 1/7] =?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 <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 18:04:08 +1000
Subject: [PATCH 2/7] 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 <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 18:15:38 +1000
Subject: [PATCH 3/7] 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 <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 18:21:22 +1000
Subject: [PATCH 4/7] 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),

From cbd925eaba14e1c73d12b7688735874472971a44 Mon Sep 17 00:00:00 2001
From: Matt Giuca <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 18:58:20 +1000
Subject: [PATCH 5/7] utilities.elements_from_state_vector: Abstract
 vector_angle to a function.

This abstracts the repeated code that finds the unsigned minimum angle between
two vectors using acos to a standalone function.
---
 src/orbital/utilities.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py
index baaac86..b64ada0 100644
--- a/src/orbital/utilities.py
+++ b/src/orbital/utilities.py
@@ -305,7 +305,7 @@ def elements_from_state_vector(r, v, mu):
         else:
             # Argument of periapsis is angle between
             # node and eccentricity vectors.
-            arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev)))
+            arg_pe = vector_angle(n, ev)
 
     if abs(e) < SMALL_NUMBER:
         if i_small:
@@ -317,7 +317,7 @@ def elements_from_state_vector(r, v, mu):
         else:
             # True anomaly is angle between node
             # vector and position vector.
-            f = acos(dot(n, r) / (norm(n) * norm(r)))
+            f = vector_angle(n, r)
             if dot(n, v) > 0:
                 f = mod(-f, 2 * pi)
     else:
@@ -326,8 +326,7 @@ def elements_from_state_vector(r, v, mu):
 
         # True anomaly is angle between eccentricity
         # vector and position vector.
-        f = acos(dot(ev, r) / (norm(ev) * norm(r)))
-
+        f = vector_angle(ev, r)
         if dot(r, v) < 0:
             f = mod(-f, 2 * pi)
 
@@ -385,6 +384,14 @@ def divmod(x, y):
     return (floor(x / y), mod(x, y))
 
 
+def vector_angle(v1, v2):
+    """Return the unsigned minimum angle between two XyzVectors [rad].
+
+    The result will always be between 0 and pi inclusive.
+    """
+    return acos(dot(v1, v2) / (norm(v1) * norm(v2)))
+
+
 # Objects for package
 
 

From 9b811671700389ee0f2acb11e9b259a92d59c2df Mon Sep 17 00:00:00 2001
From: Matt Giuca <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 18:59:31 +1000
Subject: [PATCH 6/7] utilities.elements_from_state_vector: fix nan values for
 f.

This was caused by the argument to acos being slightly higher than 1.0 due to
rounding errors, resulting in nan (when it should just be 0).

Fixes #40.
---
 src/orbital/utilities.py |  5 ++++-
 tests/test_orbital.py    | 36 +++++++++++++++++-------------------
 2 files changed, 21 insertions(+), 20 deletions(-)

diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py
index b64ada0..e8f28b3 100644
--- a/src/orbital/utilities.py
+++ b/src/orbital/utilities.py
@@ -389,7 +389,10 @@ def vector_angle(v1, v2):
 
     The result will always be between 0 and pi inclusive.
     """
-    return acos(dot(v1, v2) / (norm(v1) * norm(v2)))
+    # The parameter to acos is clamped to the range [-1.0, 1.0]. If this is not
+    # done, rounding errors can make the value slightly below -1.0 or above 1.0,
+    # resulting in nan. These should simply result in pi and 0.0, respectively.
+    return acos(np.clip(dot(v1, v2) / (norm(v1) * norm(v2)), -1.0, 1.0))
 
 
 # Objects for package
diff --git a/tests/test_orbital.py b/tests/test_orbital.py
index 3881e95..74df48c 100644
--- a/tests/test_orbital.py
+++ b/tests/test_orbital.py
@@ -769,24 +769,19 @@ def test_from_state_vector_f_at_periapsis(self):
         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)
+        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=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
@@ -861,7 +856,10 @@ def test_from_state_vector_roundtrips(self):
                 Velocity(0, 11811.441678561141, 11811.441678561141),
             ),
             # Elliptical inclined, arg_pe > 180°.
-            # (Position(0, -1767766.952966369, -1767766.952966369), Velocity(16703.901013, 0, 0)),
+            (
+                Position(0, -1767766.952966369, -1767766.952966369),
+                Velocity(16703.901013, 0, 0),
+            ),
         ]
 
         for i, (r, v) in enumerate(CASES):

From 7a6f48adf9d7b364df6de42b63ef5f0a5b3eae94 Mon Sep 17 00:00:00 2001
From: Matt Giuca <matt.giuca@gmail.com>
Date: Tue, 16 Sep 2025 22:29:45 +1000
Subject: [PATCH 7/7] from_state_vector: Fixed arg_pe bug where orbit is very
 slightly inclined.
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

If the Z coordinate is negative, but the inclination is extremely small, it
triggers the non-inclined case, which hard-codes raan to 0, yet still runs the
code path which inverts arg_pe if ev.z < 0. This code path only makes sense if
raan is also rotated 180°, which it is not.

Moved the ev.z < 0 code path into the i-not-small condition. This also makes the
code more readable, as now all the code that sets arg_pe is in one place.

Added regression test cases for this issue.

Fixes #44.
---
 src/orbital/utilities.py |  5 ++---
 tests/test_orbital.py    | 30 ++++++++++++++++++++++++++++++
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py
index e8f28b3..a9a8b48 100644
--- a/src/orbital/utilities.py
+++ b/src/orbital/utilities.py
@@ -306,6 +306,8 @@ def elements_from_state_vector(r, v, mu):
             # Argument of periapsis is angle between
             # node and eccentricity vectors.
             arg_pe = vector_angle(n, ev)
+            if ev.z < 0:
+                arg_pe = mod(-arg_pe, 2 * pi)
 
     if abs(e) < SMALL_NUMBER:
         if i_small:
@@ -321,9 +323,6 @@ def elements_from_state_vector(r, v, mu):
             if dot(n, v) > 0:
                 f = mod(-f, 2 * pi)
     else:
-        if ev.z < 0:
-            arg_pe = mod(-arg_pe, 2 * pi)
-
         # True anomaly is angle between eccentricity
         # vector and position vector.
         f = vector_angle(ev, r)
diff --git a/tests/test_orbital.py b/tests/test_orbital.py
index 74df48c..182e805 100644
--- a/tests/test_orbital.py
+++ b/tests/test_orbital.py
@@ -761,6 +761,30 @@ def test_from_state_vector_elliptical_arg_pe_gt_180(self):
         self.assertEqual(orbit.body, earth)
         self.assertAlmostEqual(orbit.t, 0.0)
 
+    def test_from_state_vector_elliptical_almost_flat(self):
+        # Elliptical orbit, at periapsis, arg_pe = 90°.
+        # Almost zero inclination, but Z is very slightly negative.
+        # Regression test for https://github.com/RazerM/orbital/issues/44
+        # The issue is that the inclination is so small, it is considered to be
+        # a non-inclined case, but the eccentricity vector's Z coordinate is
+        # very slightly negative, which would flip arg_pe erroneously.
+        R = Position(0, 2500000, -0.01)
+        V = Velocity(-16703.901013, 0, 0)
+
+        orbit = KeplerianElements.from_state_vector(R, V, body=earth)
+        numpy.testing.assert_almost_equal(orbit.r, R, decimal=2)
+        numpy.testing.assert_almost_equal(orbit.v, V, decimal=2)
+        self.assertAlmostEqual(orbit.a, 10000000.0, places=2)
+        self.assertAlmostEqual(orbit.e, 0.75)
+        self.assertAlmostEqual(orbit.i, 0.0)
+        # The absolute value of raan and arg_pe doesn't matter, because i is
+        # small. All that matters is that they are 90° apart.
+        # The current implementation always sets raan = 0 when i is small.
+        self.assertAlmostEqual(orbit.raan, 0.0)
+        self.assertAlmostEqual(orbit.arg_pe, radians(90.0))
+        self.assertAlmostEqual(orbit.M0, 0.0)
+        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.
@@ -837,6 +861,12 @@ def test_from_state_vector_roundtrips(self):
             (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), prograde, arg_pe = 90°, at periapsis.
+            (Position(0, 2500000, 0), Velocity(-16703.901013, 0, 0)),
+            # Elliptical almost-flat (e=0.75, i ~= 0 but Z slightly positive).
+            (Position(0, 2500000, 0.000000001), Velocity(-16703.901013, 0, 0)),
+            # Elliptical almost-flat (e=0.75, i ~= 0 but Z slightly negative).
+            (Position(0, 2500000, -0.000000001), Velocity(-16703.901013, 0, 0)),
             # Elliptical flat (e=0.75), arg_pe > 180°.
             (Position(0, -2500000, 0), Velocity(16703.901013, 0, 0)),
             # Elliptical flat (e=0.75), retrograde.