Merge pull request #80 from Juanlu001/osculating-elements

Compute osculating orbital elements

CHANGELOG.txt

@@ -1,5 +1,13 @@
 # orbit-predictor changelog
 
+## 1.11.0 (2019-12-27)
+
+* Add `osculating_elements` property to Position that computes
+  osculating orbital elements
+* Add `is_sun_synchronous` function that checks if a Predictor
+  represents a Sun-synchronous orbit
+* Add Earth flattening and light speed to constants
+
 ## 1.10.0 (2019-11-27)
 
 * Add get_shadow method that computes if a satellite is in shadow

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2017-2019 Satellogic SA
+Copyright (c) 2017-2020 Satellogic SA
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

orbit_predictor/constants.py

@@ -5,8 +5,11 @@
 
 AU = 149597870.700  # km
 
+LIGHT_SPEED_KMS = 299792.458  # km / s
+
 OMEGA = 2 * pi / (86400 * 365.2421897)  # rad / s
 MU_E = wgs84.mu  # km3 / s2
 R_E_KM = wgs84.radiusearthkm  # km
+F_E = 1 / 298.257223560
 J2 = wgs84.j2
 OMEGA_E = 7.292115e-5  # rad / s

orbit_predictor/locations.py

@@ -25,11 +25,10 @@
 from math import asin, cos, degrees, radians, sin, sqrt
 import os
 
+from orbit_predictor.constants import LIGHT_SPEED_KMS
 from orbit_predictor import coordinate_systems
 from orbit_predictor.utils import reify, sun_azimuth_elevation
 
-LIGHT_SPEED_KMS = 299792.458
-
 
 class Location:
     def __init__(self, name, latitude_deg, longitude_deg, elevation_m):

orbit_predictor/predictors/accurate.py

@@ -55,7 +55,7 @@
 from orbit_predictor import coordinate_systems
 
 from ..utils import reify, timetuple_from_dt
-from .base import CartesianPredictor, logger
+from .base import CartesianPredictor
 
 # Hack Zone be warned
 

orbit_predictor/predictors/base.py

@@ -26,9 +26,12 @@
 from collections import namedtuple
 from math import pi, acos, degrees, radians
 
-from orbit_predictor.exceptions import NotReachable, PropagationError
+import numpy as np
 
+from orbit_predictor.constants import MU_E
+from orbit_predictor.exceptions import NotReachable, PropagationError
 from orbit_predictor import coordinate_systems
+from orbit_predictor.keplerian import rv2coe
 from orbit_predictor.utils import (
     cross_product,
     dot_product,
@@ -55,6 +58,28 @@ def position_llh(self):
         """Latitude (deg), longitude (deg), altitude (km)."""
         return coordinate_systems.ecef_to_llh(self.position_ecef)
 
+    @reify
+    def osculating_elements(self):
+        """Osculating Keplerian orbital elements.
+
+        Semimajor axis (km), eccentricity, inclination (deg),
+        right ascension of the ascending node or RAAN (deg),
+        argument of perigee (deg), true anomaly (deg).
+
+        """
+        gmst = gstime_from_datetime(self.when_utc)
+        position_eci = coordinate_systems.ecef_to_eci(self.position_ecef, gmst)
+        velocity_eci = coordinate_systems.ecef_to_eci(self.velocity_ecef, gmst)
+
+        # Convert position to Keplerian osculating elements
+        p, ecc, inc, raan, argp, ta = rv2coe(
+            MU_E, np.array(position_eci), np.array(velocity_eci)
+        )
+        # Transform to more familiar semimajor axis
+        sma = p / (1 - ecc ** 2)
+
+        return sma, ecc, degrees(inc), degrees(raan), degrees(argp), degrees(ta)
+
 
 class PredictedPass:
     def __init__(self, location, sate_id,

orbit_predictor/predictors/keplerian.py

@@ -21,27 +21,22 @@
 # SOFTWARE.
 
 import datetime as dt
-from math import degrees, radians, sqrt
+from math import radians, sqrt
 
-import numpy as np
-from sgp4.earth_gravity import wgs84
-
-from orbit_predictor import coordinate_systems
 from orbit_predictor.angles import ta_to_M, M_to_ta
-from orbit_predictor.keplerian import rv2coe, coe2rv
+from orbit_predictor.constants import MU_E
+from orbit_predictor.keplerian import coe2rv
 from orbit_predictor.predictors import TLEPredictor
 from orbit_predictor.predictors.base import CartesianPredictor
-from orbit_predictor.utils import gstime_from_datetime, mean_motion
-
-MU_E = wgs84.mu
+from orbit_predictor.utils import mean_motion
 
 
 def kepler(argp, delta_t_sec, ecc, inc, p, raan, sma, ta):
     # Initial mean anomaly
     M_0 = ta_to_M(ta, ecc)
 
     # Mean motion
-    n = sqrt(wgs84.mu / sma ** 3)
+    n = sqrt(MU_E / sma ** 3)
 
     # Propagation
     M = M_0 + n * delta_t_sec
@@ -113,17 +108,7 @@ def from_tle(cls, sate_id, source, date=None):
         # Retrieve TLE position at given date as starting point
         pos = TLEPredictor(sate_id, source).get_position(date)
 
-        # Convert position from ECEF to ECI
-        gmst = gstime_from_datetime(date)
-        position_eci = coordinate_systems.ecef_to_eci(pos.position_ecef, gmst)
-        velocity_eci = coordinate_systems.ecef_to_eci(pos.velocity_ecef, gmst)
-
-        # Convert position to Keplerian osculating elements
-        p, ecc, inc, raan, argp, ta = rv2coe(
-            wgs84.mu, np.array(position_eci), np.array(velocity_eci))
-        sma = p / (1 - ecc ** 2)
-
-        return cls(sma, ecc, degrees(inc), degrees(raan), degrees(argp), degrees(ta), date)
+        return cls(*pos.osculating_elements, epoch=date)
 
     def propagate_eci(self, when_utc=None):
         """Return position and velocity in the given date using ECI coordinate system.

orbit_predictor/predictors/numerical.py

@@ -23,13 +23,35 @@
 from math import degrees, radians, sqrt, cos, sin
 import datetime as dt
 
+try:
+    from math import isclose
+except ImportError:
+    def isclose(a, b, *, rel_tol=1e-09, abs_tol=0.0):
+        return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
+
 import numpy as np
 
 from orbit_predictor.constants import OMEGA, MU_E, R_E_KM, J2, OMEGA_E
 from orbit_predictor.predictors.keplerian import KeplerianPredictor
 from orbit_predictor.angles import ta_to_M, M_to_ta
 from orbit_predictor.keplerian import coe2rv
-from orbit_predictor.utils import njit, raan_from_ltan, float_to_hms
+from orbit_predictor.utils import njit, raan_from_ltan, float_to_hms, mean_motion
+
+
+def is_sun_synchronous(predictor, rtol=1e-3, epoch=None):
+    """Check if predictor corresponds to Sun-synchronous orbit within tolerance.
+
+    """
+    if epoch is None:
+        epoch = dt.datetime.now()
+
+    sma_km, ecc, inc_deg, *_ = predictor.get_position(epoch).osculating_elements
+    p = sma_km * (1 - ecc ** 2)
+    n = mean_motion(sma_km)
+
+    raan_dot_sec = - 3 * n * R_E_KM ** 2 * J2 / (2 * p ** 2) * cos(radians(inc_deg))
+
+    return isclose(raan_dot_sec, OMEGA, rel_tol=rtol)
 
 
 def sun_sync_plane_constellation(num_satellites, *,
@@ -222,7 +244,9 @@ def sun_synchronous(cls, *, alt_km=None, ecc=None, inc_deg=None, ltan_h=12, date
                     )
 
         except FloatingPointError as e:
-            raise InvalidOrbitError("Cannot find Sun-synchronous orbit with given parameters") from e
+            raise InvalidOrbitError(
+                "Cannot find Sun-synchronous orbit with given parameters"
+            ) from e
 
         # TODO: Allow change in time or location
         # Right the epoch is fixed given the LTAN, as well as the sub-satellite point

orbit_predictor/sources.py

@@ -22,7 +22,6 @@
 
 import logging
 from collections import defaultdict, namedtuple
-import warnings
 
 import requests
 from urllib import parse as urlparse

orbit_predictor/utils.py

@@ -26,11 +26,10 @@
 from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt, tan, modf
 
 import numpy as np
-from sgp4.earth_gravity import wgs84
 from sgp4.ext import jday
 from sgp4.propagation import _gstime
 
-from .constants import AU
+from .constants import AU, R_E_KM, MU_E
 from .coordinate_systems import eci_to_radec, ecef_to_eci
 
 # Inspired in https://github.com/poliastro/poliastro/blob/88edda8/src/poliastro/jit.py
@@ -322,7 +321,7 @@ def get_shadow(r, when_utc):
     return shadow(r_sun, ecef_to_eci(r, gmst))
 
 
-def shadow(r_sun, r, r_p=wgs84.radiusearthkm):
+def shadow(r_sun, r, r_p=R_E_KM):
     """
     Gives illumination of Earth satellite (2 for illuminated, 1 for penumbra, 0 for umbra).
 
@@ -405,7 +404,7 @@ def timetuple_from_dt(when_utc):
 
 
 def mean_motion(sma_km):
-    return sqrt(wgs84.mu / sma_km ** 3)  # rad / s
+    return sqrt(MU_E / sma_km ** 3)  # rad / s
 
 
 class reify:

orbit_predictor/version.py

@@ -1,2 +1,2 @@
 # https://www.python.org/dev/peps/pep-0440/
-__version__ = '2.0.dev0'
+__version__ = '1.11.0'

setup.py

@@ -2,7 +2,7 @@
 import os.path
 from setuptools import setup, find_packages
 
-# Copyright 2017-2019 Satellogic SA.
+# Copyright 2017-2020 Satellogic SA.
 
 
 # https://packaging.python.org/guides/single-sourcing-package-version/

tests/test_keplerian.py

@@ -7,8 +7,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-from sgp4.earth_gravity import wgs84
-
+from orbit_predictor.constants import MU_E
 from orbit_predictor.keplerian import coe2rv, rv2coe
 
 
@@ -25,7 +24,7 @@ def test_convert_coe_to_rv(self):
         expected_position = [6525.344, 6861.535, 6449.125]
         expected_velocity = [4.902276, 5.533124, -1.975709]
 
-        position, velocity = coe2rv(wgs84.mu, p, ecc, inc, raan, argp, ta)
+        position, velocity = coe2rv(MU_E, p, ecc, inc, raan, argp, ta)
 
         assert_allclose(position, expected_position, rtol=1e-5)
         assert_allclose(velocity, expected_velocity, rtol=1e-5)
@@ -44,7 +43,7 @@ def test_convert_rv_to_coe(self):
         expected_argp = radians(53.38)
         expected_ta = radians(92.335)
 
-        p, ecc, inc, raan, argp, ta = rv2coe(wgs84.mu, position, velocity)
+        p, ecc, inc, raan, argp, ta = rv2coe(MU_E, position, velocity)
 
         self.assertAlmostEqual(p, expected_p, places=0)
         self.assertAlmostEqual(ecc, expected_ecc, places=4)

tests/test_numerical_predictor.py

@@ -6,7 +6,9 @@
 import pytest
 
 from orbit_predictor.locations import ARG
-from orbit_predictor.predictors.numerical import J2Predictor, InvalidOrbitError, R_E_KM
+from orbit_predictor.predictors.numerical import (
+    J2Predictor, InvalidOrbitError, R_E_KM, is_sun_synchronous
+)
 
 
 class J2PredictorTests(TestCase):
@@ -53,47 +55,47 @@ def test_sun_sync_from_altitude_and_eccentricity(self):
         expected_inc = 98.6
 
         pred = J2Predictor.sun_synchronous(alt_km=800, ecc=0)
-        self.assertAlmostEqual(pred._inc, expected_inc, places=2)
+        self.assertAlmostEqual(pred.get_position().osculating_elements[2], expected_inc, places=2)
 
     def test_sun_sync_from_altitude_and_inclination(self):
         # Hardcoded from our implementation
         expected_ecc = 0.14546153131334466
 
         pred = J2Predictor.sun_synchronous(alt_km=475, inc_deg=97)
-        self.assertAlmostEqual(pred._ecc, expected_ecc, places=16)
+        self.assertAlmostEqual(pred.get_position().osculating_elements[1], expected_ecc, places=15)
 
     def test_sun_sync_from_eccentricity_and_inclination(self):
         # Vallado 3rd edition, example 11-2
         expected_sma = 7346.846
 
         pred = J2Predictor.sun_synchronous(ecc=0.2, inc_deg=98.6)
-        self.assertAlmostEqual(pred._sma, expected_sma, places=1)
+        self.assertAlmostEqual(pred.get_position().osculating_elements[0], expected_sma, places=1)
 
     def test_sun_sync_delta_true_anomaly_has_expected_anomaly_and_epoch(self):
         date = dt.datetime.today().date()
         ltan_h = 12
         expected_ref_epoch = dt.datetime(date.year, date.month, date.day, 12)
 
-        for ta_deg in [-30, 0, 30]:
+        for expected_ta_deg in [-30, 0, 30]:
             pred = J2Predictor.sun_synchronous(
-                alt_km=800, ecc=0, date=date, ltan_h=ltan_h, ta_deg=ta_deg
+                alt_km=800, ecc=0, date=date, ltan_h=ltan_h, ta_deg=expected_ta_deg
             )
 
-            self.assertEqual(pred._ta, ta_deg)
-            self.assertEqual(pred._epoch, expected_ref_epoch)
+            ta_deg = pred.get_position(expected_ref_epoch).osculating_elements[5]
+            self.assertAlmostEqual(ta_deg, expected_ta_deg % 360, places=12)
 
     def test_sun_sync_delta_true_anomaly_non_circular(self):
         date = dt.datetime.today().date()
         ltan_h = 12
         expected_ref_epoch = dt.datetime(date.year, date.month, date.day, 12)
 
-        for ta_deg in [-30, 30]:
+        for expected_ta_deg in [-30, 30]:
             pred = J2Predictor.sun_synchronous(
-                alt_km=475, ecc=0.1455, date=date, ltan_h=ltan_h, ta_deg=ta_deg
+                alt_km=475, ecc=0.1455, date=date, ltan_h=ltan_h, ta_deg=expected_ta_deg
             )
 
-            self.assertEqual(pred._ta, ta_deg)
-            self.assertEqual(pred._epoch, expected_ref_epoch)
+            ta_deg = pred.get_position(expected_ref_epoch).osculating_elements[5]
+            self.assertAlmostEqual(ta_deg, expected_ta_deg % 360, places=12)
 
 
 # Test data from Wertz et al. "Space Mission Engineering: The New SMAD" (2011), table 9-13
@@ -108,4 +110,27 @@ def test_sun_sync_delta_true_anomaly_non_circular(self):
 def test_repeated_groundtrack_sma(orbits, days, inc_deg, expected_h):
     pred = J2Predictor.repeating_ground_track(orbits=orbits, days=days, ecc=0.0, inc_deg=inc_deg)
 
-    assert_almost_equal(pred._sma - R_E_KM, expected_h, decimal=0)
+    assert_almost_equal(pred.get_position().osculating_elements[0] - R_E_KM, expected_h, decimal=0)
+
+
+def test_is_sun_sync_returns_false_for_non_sun_sync_orbit():
+    pred1 = J2Predictor(7000, 0, 0, 0, 0, 0, dt.datetime.now())
+
+    assert not is_sun_synchronous(pred1)
+
+
+def test_is_sun_sync_detects_almost_sun_sync_orbit():
+    pred2 = J2Predictor(R_E_KM + 460, 0.001, 97.4, 0, 0, 0, dt.datetime.now())
+
+    assert not is_sun_synchronous(pred2)
+    assert is_sun_synchronous(pred2, rtol=1e-1)
+
+
+def test_is_sun_sync_returns_true_for_sun_sync_orbit():
+    pred1 = J2Predictor.sun_synchronous(alt_km=500, ecc=0)
+    pred2 = J2Predictor.sun_synchronous(alt_km=500, inc_deg=97)
+    pred3 = J2Predictor.sun_synchronous(ecc=0, inc_deg=97)
+
+    assert is_sun_synchronous(pred1)
+    assert is_sun_synchronous(pred2)
+    assert is_sun_synchronous(pred3)