diff --git a/src/orbital/utilities.py b/src/orbital/utilities.py index d6a9469..524c5f5 100644 --- a/src/orbital/utilities.py +++ b/src/orbital/utilities.py @@ -284,16 +284,23 @@ 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 = 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. - 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: @@ -301,23 +308,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 083d6ab..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) @@ -761,15 +755,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) @@ -834,19 +826,31 @@ 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°. - # (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 +1081,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):