Setup SGP4 Propagator with unit testing

gs/backend/positioning/sgp4_handler.py

@@ -0,0 +1,119 @@
+from dataclasses import dataclass
+from datetime import datetime
+
+from sgp4.api import SGP4_ERRORS, Satrec, jday
+from sgp4.api import WGS72 as GRAVITY_MODEL
+from sgp4.api import accelerated
+from loguru import logger
+from math import radians
+
+from gs.backend.positioning.tle import AFSPC_MODE_IMPROVED, TLEData
+
+#log whether accelerated mode is used
+if not accelerated:
+    logger.warning(
+        "SGP4 accelerated mode is not available. "
+        "Falling back to pure Python implementation (slower)."
+    )
+else:
+    logger.info("SGP4 accelerated mode is enabled.")
+
+
+@dataclass
+class SGP4Data:
+    """Data structure representing the satellite's position and velocity."""
+
+    position_km: tuple[float, float, float]
+    velocity_km_sec: tuple[float, float, float]
+
+
+def setup_sgp4(tle: TLEData) -> Satrec:
+    """
+    Initialize the SGP4 satellite model using TLE data.  Formatting and SGP4 initialization pulled from link below
+        https://pypi.org/project/sgp4/
+    """
+
+
+    sat = Satrec()
+
+
+    sat.sgp4init( #causes error when tle.eccentricity is low while tle.drag_term is high)
+        GRAVITY_MODEL,  # gravity model
+        AFSPC_MODE_IMPROVED,  # propagation mode
+        tle.satellite_number,  # satellite number
+        tle.convert_epoch_values_to_jd(),  # epoch (Julian date)
+        tle.drag_term,  # BSTAR drag term
+        #6.2485e-05,
+        tle.first_derivative_mean_motion,  # first time derivative of mean motion
+        tle.second_derivative_mean_motion,  # second time derivative of mean motion
+        tle.eccentricity,  # eccentricity
+        radians(tle.argument_of_perigee),  # argument of perigee (radians)
+        radians(tle.inclination),  # inclination (radians) 
+        radians(tle.mean_anomaly),  # mean anomaly (radians)
+        tle.mean_motion * (2 * 3.141592653589793 / 1440.0),  # mean motion (radians/min)
+        radians(tle.right_ascension),  # RA of ascending node (radians)
+    )
+    """  These values where used for debugging to find source of eccentricity value, can use to test if issue is fixed
+    sat.sgp4init(
+        GRAVITY_MODEL,                # gravity model
+        'i',                  # 'a' = old AFSPC mode, 'i' = improved mode
+        25544,                # satnum: Satellite number
+        25545.69339541,       # epoch: days since 1949 December 31 00:00 UT
+        3.8792e-05,           # bstar: drag coefficient (1/earth radii)
+        0.0,                  # ndot: first time derivative of mean motion (radians/min^2)
+        0.0,                  # nddot: second derivative of mean motion (radians/min^3)
+        0.0007417,            # ecco: eccentricity (0..1)
+        0.3083420829620822,   # argpo: argument of perigee (radians)
+        0.9013560935706996,   # inclo: inclination (radians)
+        1.4946964807494398,   # mo: mean anomaly (radians)
+        0.06763602333248933,  # no_kozai: mean motion (radians/min)
+        3.686137125541276,    # nodeo: R.A. of ascending node (radians)
+    )
+    """
+
+    """
+    print("GRAVITY_MODEL:", GRAVITY_MODEL)
+    print("Propagation mode:", AFSPC_MODE_IMPROVED)
+    print("Satellite number:", tle.satellite_number)
+    print("Epoch (Julian date):", tle.convert_epoch_values_to_jd())
+    print("BSTAR drag term:", tle.drag_term)
+    print("First derivative of mean motion:", tle.first_derivative_mean_motion)
+    print("Second derivative of mean motion:", tle.second_derivative_mean_motion)
+    print("Eccentricity:", tle.eccentricity)
+    print("Argument of perigee (rad):", radians(tle.argument_of_perigee))
+    print("Inclination (rad):", radians(tle.inclination))
+    print("Mean anomaly (rad):", radians(tle.mean_anomaly))
+    print("Mean motion (rad/min):", tle.mean_motion * (2 * 3.141592653589793 / 1440.0))
+    print("RA of ascending node (rad):", radians(tle.right_ascension))
+    """
+
+    return sat
+
+
+def get_sat_position(tle: TLEData, dt: datetime) -> SGP4Data:
+    """
+    Compute the satellite's position and velocity at a given time.
+
+        Arguments are
+        tle(TLEData): Two-line element set representing the satellite.
+        dt(datetime): The timestamp for which to calculate the position.
+
+    """
+    sat = setup_sgp4(tle)
+    jd, fr = jday(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)
+    error_code, position, velocity = sat.sgp4(jd, fr)
+
+    #print("eccentricity (parsed):", tle.eccentricity)
+
+
+    if 0.0 <= tle.eccentricity and tle.eccentricity <= 1.0:
+        print("tle.eccentricity within expected bounds")
+
+    error_message = SGP4_ERRORS.get(error_code, None)
+    if error_message is not None:
+        raise RuntimeError(error_message)
+
+    return SGP4Data(
+        position_km=tuple(position),
+        velocity_km_sec=tuple(velocity),
+    )

gs/backend/positioning/tle.py

@@ -0,0 +1,209 @@
+from dataclasses import dataclass
+from datetime import datetime, timedelta
+from typing import Final
+
+import requests
+
+from gs.backend.sun.ephemeris import convert_date_to_jd
+
+TLE_FAILURE_RETURN_MESSAGE: Final[str] = "No GP data found"
+AFSPC_MODE_IMPROVED: Final[str] = "i"
+
+
+@dataclass
+class TLEData:
+    """Class for TLE data."""
+
+    name: str
+    satellite_number: int
+    classification: str
+    launch_year: int
+    launch_number: int
+    launch_piece: str
+    epoch_year: int  # Last two digits of the year
+    epoch_day: float  # Day of the year and fractional portion of the day
+    first_derivative_mean_motion: float
+    second_derivative_mean_motion: float
+    drag_term: float
+    ephemeris_type: int
+    element_number: int
+    inclination: float  # In degrees
+    right_ascension: float  # In degrees
+    eccentricity: float
+    argument_of_perigee: float  # In degrees
+    mean_anomaly: float  # In degrees
+    mean_motion: float  # In revolutions per day
+    revolution_number: int
+
+    def to_tle(self) -> str:
+        """
+        Convert the TLE data into a valid 2-line TLE string format
+        Formatting for this is pulled from this resource: https://www.researchgate.net/figure/Two-line-Element-Set-Format-An-example-TLE-is-shown-with-descriptions-and-units-of-each_fig3_289774073
+        """
+        first_line = (
+            f"1 {self.satellite_number:05d}{self.classification} "
+            f"{self.launch_year:02d}{self.launch_number:03d}{self.launch_piece:<3} "
+            f"{self.epoch_year:02d}{self.epoch_day:012.8f} "
+            f".{self.first_derivative_mean_motion:.8f}".split(".")[1]
+            + " "
+            + f"{self.second_derivative_mean_motion: 8.5e} "
+            f"{self.drag_term: 8.5e} "
+            f"{self.ephemeris_type} "
+            f"{self.element_number:4d}"
+        )
+        second_line = (
+            f"2 {self.satellite_number:05d} "
+            f"{self.inclination:8.4f} "
+            f"{self.right_ascension:8.4f} "
+            f"{self.eccentricity:.7f}".split(".")[1]
+            + " "
+            + f"{self.argument_of_perigee:8.4f} "
+            f"{self.mean_anomaly:8.4f} "
+            f"{self.mean_motion:11.8f}{self.revolution_number:5d}"
+        )
+        first_line += str(calculate_checksum(first_line))
+        second_line += str(calculate_checksum(second_line))
+        #print("hi!")
+        #print(f"{first_line}\n{second_line}")
+        return f"{first_line}\n{second_line}"
+
+    def format_epoch_to_date(self, epoch_year: int, epoch_day: float) -> str:
+        """
+            Convert epoch year and day-of-year to a calendar date (YYYY-MM-DD).
+
+        Arguments:
+            epoch_year (int): Last two digits of the epoch year.
+            epoch_day (float): Day of the year, including fraction.
+
+        Returns:
+            str: Formatted date string in YYYY-MM-DD format.
+        """
+        full_year = 2000 + epoch_year if epoch_year < 57 else 1900 + epoch_year
+        base_date = datetime(full_year, 1, 1) + timedelta(days=epoch_day - 1)
+        return base_date.strftime("%Y-%m-%d")
+
+    def convert_epoch_values_to_jd(self) -> float:
+        """converting epoch values to julian date time"""
+        date = self.format_epoch_to_date(self.epoch_year, self.epoch_day)
+        # if(is_valid_date(date)):
+        return convert_date_to_jd(date)
+
+
+def calculate_checksum(line: str) -> int:
+    """Calculate the checksum for a line of TLE data."""
+    output = 0
+    for i in line[:68]:
+        if i.isdigit():
+            output += int(i)
+        if i == "-":
+            output += 1
+    return output % 10
+
+
+def convert_decimal_point_assumed(value: str) -> float:
+    """Convert a string to a float, assuming a decimal point before the last digit."""
+    return float(f"{value[0]}.{value[1:6]}e{value[6:]}")
+
+
+def parse_tle_data(tle_data: str) -> TLEData:
+    """
+    @brief Parse TLE data from the argument into an object. Validates the data.
+
+    @throw ValueError if the TLE data is invalid.
+
+    @param tle_data: The TLE data to parse.
+    @return A TLEData object.
+    """
+    # Check for failure
+    if tle_data == TLE_FAILURE_RETURN_MESSAGE:
+        raise ValueError("No GP data found")
+
+    lines = tle_data.splitlines()
+    name, line1, line2 = lines
+
+    # Validate the lines
+    if (
+        len(lines) != 3
+        or not line1.startswith("1")
+        or not line2.startswith("2")
+        or len(line1) != 69
+        or len(line2) != 69
+    ):
+        raise ValueError("Invalid TLE data")
+
+    # Get the checksums
+    checksum_1 = int(line1[68])
+    checksum_2 = int(line2[68])
+
+    # Parse the data
+    output = TLEData(
+        name=name.strip(),
+        satellite_number=int(line1[2:7]),
+        classification=line1[7],
+        launch_year=int(line1[9:11]),
+        launch_number=int(line1[11:14]),
+        launch_piece=line1[14:17].strip(),
+        epoch_year=int(line1[18:20]),
+        epoch_day=float(line1[20:32]),
+        first_derivative_mean_motion=float(line1[33:43]),
+        second_derivative_mean_motion=convert_decimal_point_assumed(line1[44:52]),
+        drag_term=convert_decimal_point_assumed(line1[53:61]),
+        ephemeris_type=int(line1[62]),
+        element_number=int(line1[64:68]),
+        inclination=float(line2[8:16]),
+        right_ascension=float(line2[17:25]),
+        eccentricity=float(f"0.{line2[26:33]}"),
+        argument_of_perigee=float(line2[34:42]),
+        mean_anomaly=float(line2[43:51]),
+        mean_motion=float(line2[52:63]),
+        revolution_number=int(line2[63:68]),
+    )
+
+    # Validate the checksums
+    if checksum_1 != calculate_checksum(line1):
+        raise ValueError("Invalid checksum for line 1")
+    if checksum_2 != calculate_checksum(line2):
+        raise ValueError("Invalid checksum for line 2")
+    return output
+
+
+def get_tle_data(object_id: int) -> TLEData:
+    """
+    @brief Get TLE data for a satellite from Celestrak based on the object ID.
+
+    @throw ValueError if the object ID is invalid or the TLE data is invalid.
+    """
+    url = f"https://celestrak.org/NORAD/elements/gp.php?CATNR={object_id}&FORMAT=tle"
+    response = requests.get(url)
+    output = response.text
+    return parse_tle_data(output)
+
+
+
+
+def id_from_user() -> int:
+    """
+    Get the object ID from the user.
+    This function will keep asking for input until a valid number is entered.
+    Used for testing purposes.
+    """
+    while True:
+        try:
+            object_id = int(input("Enter the satellite ID: "))
+            return object_id
+        except ValueError:
+            print("Invalid input. Please enter a valid number.")
+
+
+def main() -> None:
+    """Example usage of the get_tle_data function."""
+    object_id = id_from_user()
+    try:
+        tle_data = get_tle_data(object_id)
+        print(tle_data)
+    except ValueError as e:
+        print(e)
+
+
+if __name__ == "__main__":
+    main()

python_test/positioning/test_sgp4_handler.py

@@ -0,0 +1,57 @@
+import pytest
+from datetime import datetime
+from gs.backend.positioning.sgp4_handler import (
+    setup_sgp4,
+    get_sat_position,
+    SGP4Data
+)
+from gs.backend.positioning.tle import parse_tle_data
+from sgp4.api import Satrec
+from sgp4.conveniences import dump_satrec
+from sgp4.earth_gravity import wgs72
+from sgp4.io import twoline2rv
+from sys import stdout
+
+
+def test_setup_sgp4(): #test setup_sgp4 returns a valid Satrec object
+    tle_str = (
+        "ISS (ZARYA)\n"
+        "1 25544U 98067A   23058.54791667  .00016717  00000+0  10270-3 0  9997\n"
+        "2 25544  51.6435 143.0464 0004767 278.8055  81.2436 15.49815308274053"
+    )
+    tle = parse_tle_data(tle_str)
+    sat = setup_sgp4(tle)
+    assert isinstance(sat, Satrec)
+
+
+def test_get_sat_position(): #test that output has expected structure and data types
+    tle_str = (
+        "CSS (MENGTIAN)\n"
+        "1 54216U 22143A   25308.16010129  .00032649  00000-0  38953-3 0  9993\n"
+        "2 54216  41.4668 251.0039 0006377 319.7576  40.2790 15.60341380158647"
+    )
+    #Note that the above tle causes eccentricity out of range bug, while tle_str2 does not. 
+    #Current implementation causes error when tle.eccentricity is low while tle.drag_term is high
+    tle_str2 = (
+        "ISS (ZARYA)\n"
+        "1 25544U 98067A   25308.35786713  .00010709  00000-0  19707-3 0  9992\n"  
+        "2 25544  51.6336 332.4903 0005031  16.0382 344.0765 15.49743270536903"
+
+    )
+    tle = parse_tle_data(tle_str)
+
+    """
+    sat = setup_sgp4(tle)
+    stdout.writelines(dump_satrec(sat))
+    assert 0.0 <= tle.eccentricity < 1.0, "parsed eccentricity out of range!" #attempt to test via direct check rather than built in error handling
+    """
+
+
+    dt = datetime(2025, 11, 12, 12, 0, 0)
+    data = get_sat_position(tle, dt)
+    assert isinstance(data, SGP4Data) 
+
+    assert len(data.position_km) == 3  #test that position and velocity have x y z components (3D, length 3)
+    assert len(data.velocity_km_sec) == 3
+    assert all(isinstance(x, float) for x in data.position_km) #loop through position_km and velocity_km,
+    assert all(isinstance(x, float) for x in data.velocity_km_sec) #   ensuring data types are floats