resolve pr#336 (except for one bug)

Author: kyle-su1

gs/backend/positioning/sgp4.py

@@ -0,0 +1,115 @@
+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)
+    )
+    """
+    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/sgp4_handler.py

@@ -63,6 +63,8 @@ def to_tle(self) -> str:
         )
         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:
@@ -177,6 +179,7 @@ def get_tle_data(object_id: int) -> TLEData:
     return parse_tle_data(output)
 
 
+
 # Testing section
 
 

gs/backend/positioning/tle.py

@@ -0,0 +1,52 @@
+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"
+    )
+
+    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))
+    # ensure it's between 0 and 1
+    #assert 0.0 <= tle.eccentricity < 1.0, "parsed eccentricity out of range!"
+    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

python_test/positioning/test_sgp4_handler.py

@@ -9,7 +9,12 @@
 @pytest.mark.parametrize(
     "tle_line, expected_checksum",
     [
-        ("1 25544U 98067A 23058.54791667 .00016717 00000+0 10270-3 0 9991", 8),
+        ("1 25544U 98067A 23058.54791667 .00016717 00000+0 10270-3 0 9991", 8), #sample value
+        ("", 0), #empty line
+        ("ABCDEFG-+", 1), #check addition on minus
+        ("1 25544U---", (1+2+5+5+4+4 + 3) % 10), #combination of digits and minus
+
+
     ],
 )
 def test_calculate_checksum(tle_line, expected_checksum):
@@ -21,9 +26,17 @@ def test_calculate_checksum(tle_line, expected_checksum):
     [
         ("123456-3", 0.00123456),  # 1.23456e-3
         ("678901+2", 678.901),  # 6.78901e2
-        # ("-456789-5", -5.6789e-6),       # -0.56789e-5
+        ("000001+0", 0.00001),  
+        ("999999-1", 0.999999),  # 9.99999e-1   
+        ("123456+0", 1.23456),  # 1.23456e0
+
+        # ("-456789-5", -5.6789e-6),       # -0.56789e-5 (negative value, fails)
     ],
 )
+#I did some research and it seems that the function won't be called for negative values for real TLEs, so failing on negative values should be fine
+
+
+
 def test_convert_decimal_point_assumed(input_value, expected_result):
     assert convert_decimal_point_assumed(input_value) == pytest.approx(expected_result, rel=1e-6)
 
@@ -38,5 +51,22 @@ def test_parse_tle_data():
     parsed_data = parse_tle_data(tle_str)
 
     assert parsed_data.satellite_number == 25544
+    assert parsed_data.first_derivative_mean_motion == pytest.approx(.00016717, rel=1e-4)
     assert parsed_data.inclination == pytest.approx(51.6435, rel=1e-4)
     assert parsed_data.right_ascension == pytest.approx(143.0464, rel=1e-4)
+
+
+
+def test_parse_tle_data_invalid(): #test for corrupted TLEs with mismatched checksums
+    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"
+    )
+
+    bad_tle1 = tle_str.replace("9997", "9998")
+    bad_tle2 = tle_str.replace("15.49815308274053", "15.49815308274054")
+    with pytest.raises(ValueError, match="Invalid checksum for line 1"):
+        parse_tle_data(bad_tle1)
+    with pytest.raises(ValueError, match="Invalid checksum for line 2"):
+        parse_tle_data(bad_tle2)