Merge e19fdd5 into f47ae93

orbit_predictor/groundtrack.py

@@ -0,0 +1,24 @@
+class BaseElevationAPI:
+    def get_elevation(self, *, longitude, latitude):
+        raise NotImplementedError
+
+    def get_ground_point(self, position):
+        latitude, longitude, _ = position.position_llh
+        return (
+            longitude,
+            latitude,
+            self.get_elevation(longitude=longitude, latitude=latitude),
+        )
+
+    def get_groundtrack(self, positions):
+        return [self.get_ground_point(position) for position in positions]
+
+
+class ZeroElevation(BaseElevationAPI):
+    def get_elevation(self, *, longitude, latitude):
+        return 0.0
+
+
+def compute_groundtrack(predictor, times, elevation_api=ZeroElevation()):
+    positions = [predictor.get_position(time) for time in times]
+    return elevation_api.get_groundtrack(positions)

orbit_predictor/predictors/accurate.py

@@ -101,24 +101,15 @@ def _propagate_only_position_ecef(self, timetuple):
         gmst = _gstime(jday(*timetuple))
         return coordinate_systems.eci_to_ecef(position_eci, gmst)
 
-    def _propagate_eci(self, when_utc=None):
-        """Return position and velocity in the given date using ECI coordinate system."""
-        tle = self.source.get_tle(self.sate_id, when_utc)
-        logger.debug("Propagating using ECI. sate_id: %s, when_utc: %s, tle: %s",
-                     self.sate_id, when_utc, tle)
-        tle_line_1, tle_line_2 = tle.lines
-        sgp4_sate = twoline2rv(tle_line_1, tle_line_2, wgs84)
-        timetuple = when_utc.timetuple()[:6]
-        timetuple[5] = timetuple[5] + when_utc.microsecond * 1e-6
-        position_eci, velocity_eci = sgp4_sate.propagate(*timetuple)
-        return position_eci, velocity_eci
-
     def _propagate_ecef(self, when_utc):
         """Return position and velocity in the given date using ECEF coordinate system."""
         timetuple = (when_utc.year, when_utc.month, when_utc.day,
                      when_utc.hour, when_utc.minute, when_utc.second + when_utc.microsecond * 1e-6)
 
         position_eci, velocity_eci = self.propagator.propagate(*timetuple)
+        if self.propagator.error != 0:
+            raise RuntimeError(self.propagator.error_message)
+
         gmst = _gstime(jday(*timetuple))
         position_ecef = coordinate_systems.eci_to_ecef(position_eci, gmst)
         velocity_ecef = coordinate_systems.eci_to_ecef(velocity_eci, gmst)

orbit_predictor/predictors/keplerian.py

@@ -25,7 +25,6 @@
 
 import numpy as np
 from sgp4.earth_gravity import wgs84
-from sgp4.io import twoline2rv
 
 from orbit_predictor import coordinate_systems
 from orbit_predictor.angles import ta_to_M, M_to_ta
@@ -111,14 +110,11 @@ def from_tle(cls, sate_id, source, date=None):
         if date is None:
             date = dt.datetime.utcnow()
 
-        tle = source.get_tle(sate_id, date)
-
-        # Retrieve TLE epoch and corresponding position
-        epoch = twoline2rv(tle.lines[0], tle.lines[1], wgs84).epoch
-        pos = TLEPredictor(sate_id, source).get_position(epoch)
+        # 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(epoch)
+        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)
 
@@ -127,7 +123,7 @@ def from_tle(cls, sate_id, source, date=None):
             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), epoch)
+        return cls(sma, ecc, degrees(inc), degrees(raan), degrees(argp), degrees(ta), date)
 
     def _propagate_eci(self, when_utc):
         """Return position and velocity in the given date using ECI coordinate system.

orbit_predictor/predictors/numerical.py

@@ -99,13 +99,13 @@ def repeating_ground_track_sma(orbits, days=1, *, ecc, inc_deg=0, tolerance=1e-8
     while True:
         sma_new = np.cbrt(MU_E * (1 / n) ** 2)
         p = sma_new * (1 - ecc ** 2)
-        Omega_dot = - 3 * n * J2 / 2 * (R_E_KM / p) ** 2 * np.cos(np.radians(inc_deg))
-        omega_dot = 3 * n * J2 / 4 * (R_E_KM / p) ** 2 * (4 - 5 * np.sin(np.radians(inc_deg)) ** 2)
+        node_dot = - 3 * n * J2 / 2 * (R_E_KM / p) ** 2 * np.cos(np.radians(inc_deg))
+        argp_dot = 3 * n * J2 / 4 * (R_E_KM / p) ** 2 * (4 - 5 * np.sin(np.radians(inc_deg)) ** 2)
         M0_dot = (
             3 * n * J2 / 4 * (R_E_KM / p) ** 2 * np.sqrt(1 - ecc ** 2)
             * (2 - 3 * np.sin(np.radians(inc_deg)) ** 2)
         )
-        n = k * (OMEGA_E - Omega_dot) - (M0_dot + omega_dot)
+        n = k * (OMEGA_E - node_dot) - (M0_dot + argp_dot)
         sma = np.cbrt(MU_E * (1 / n) ** 2)
         if np.isclose(sma, sma_new, rtol=tolerance):
             break

tests/test_keplerian_predictor.py

@@ -122,6 +122,35 @@ def test_from_tle_returns_same_initial_conditions_on_epoch(self):
         assert_allclose(pos_keplerian.position_ecef, pos_tle.position_ecef, rtol=1e-9)
         assert_allclose(pos_keplerian.velocity_ecef, pos_tle.velocity_ecef, rtol=1e-13)
 
+    def test_from_tle_returns_same_initial_conditions_tle_date(self):
+        start = dt.datetime(2017, 3, 6, 7, 51)
+        db = MemoryTLESource()
+        db.add_tle(self.SATE_ID, self.LINES, start)
+
+        keplerian_predictor = KeplerianPredictor.from_tle(self.SATE_ID, db, start)
+        tle_predictor = TLEPredictor(self.SATE_ID, db)
+
+        pos_keplerian = keplerian_predictor.get_position(start)
+        pos_tle = tle_predictor.get_position(start)
+
+        assert_allclose(pos_keplerian.position_ecef, pos_tle.position_ecef, rtol=1e-9)
+        assert_allclose(pos_keplerian.velocity_ecef, pos_tle.velocity_ecef, rtol=1e-13)
+
+    def test_from_tle_returns_same_initial_conditions_future_date(self):
+        start = dt.datetime(2017, 3, 6, 7, 51)
+        future = dt.datetime(2017, 9, 6, 7, 51)
+        db = MemoryTLESource()
+        db.add_tle(self.SATE_ID, self.LINES, start)
+
+        keplerian_predictor = KeplerianPredictor.from_tle(self.SATE_ID, db, future)
+        tle_predictor = TLEPredictor(self.SATE_ID, db)
+
+        pos_keplerian = keplerian_predictor.get_position(future)
+        pos_tle = tle_predictor.get_position(future)
+
+        assert_allclose(pos_keplerian.position_ecef, pos_tle.position_ecef, rtol=1e-9)
+        assert_allclose(pos_keplerian.velocity_ecef, pos_tle.velocity_ecef, rtol=1e-13)
+
     def test_from_tle_with_no_date_works(self):
         # https://github.com/satellogic/orbit-predictor/issues/52
         start = dt.datetime(2017, 3, 6, 7, 51)