diff --git a/gs/backend/positioning/__init__.py b/gs/backend/positioning/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/gs/backend/positioning/sgp4_handler.py b/gs/backend/positioning/sgp4_handler.py new file mode 100644 index 000000000..6a179b44c --- /dev/null +++ b/gs/backend/positioning/sgp4_handler.py @@ -0,0 +1,128 @@ +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. + :param position_km: tuple representing the position of the satellite in km + :param velocity_km_sec: tuple representing velocity of the satellite in km/second + + + """ + + 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/ + + :warning: currently broken for TLEs where eccentricity is low relative to drag term + :param tle(TLEData): The TLE string used to initialize Satrec + :return: Returns initialized Satrec + """ + + + 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. + + :warning: currently broken for TLEs where eccentricity is low relative to drag term + + :param tle(TLEData): Two-line element set representing the satellite. + :param dt(datetime): The timestamp for which to calculate the position. + :return: Returns location data with custom SGP4Data object + """ + + 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) + + 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), + ) diff --git a/gs/backend/positioning/tle.py b/gs/backend/positioning/tle.py new file mode 100644 index 000000000..791a89350 --- /dev/null +++ b/gs/backend/positioning/tle.py @@ -0,0 +1,216 @@ +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)) + 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). + + :param epoch_year (int): Last two digits of the epoch year. + :param epoch_day (float): Day of the year, including fraction. + + :return: Returns formatted date string in YYYY-MM-DD format as a string. + """ + 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. + + :param line(str): The single line of TLE data to calculate the checksum for + :return: Returns the checksum integer + + """ + 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. + + :param value(str): the string to convert + :returns: data from string converted into float + """ + 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() diff --git a/python_test/positioning/test_sgp4_handler.py b/python_test/positioning/test_sgp4_handler.py new file mode 100644 index 000000000..96d071f84 --- /dev/null +++ b/python_test/positioning/test_sgp4_handler.py @@ -0,0 +1,51 @@ +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) + 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 + +""" \ No newline at end of file diff --git a/python_test/positioning/test_tle.py b/python_test/positioning/test_tle.py new file mode 100644 index 000000000..9d231e9a7 --- /dev/null +++ b/python_test/positioning/test_tle.py @@ -0,0 +1,72 @@ +import pytest +from gs.backend.positioning.tle import ( + calculate_checksum, + convert_decimal_point_assumed, + parse_tle_data, +) + + +@pytest.mark.parametrize( + "tle_line, expected_checksum", + [ + ("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): + assert calculate_checksum(tle_line) == expected_checksum + + +@pytest.mark.parametrize( # The convert_decimal_point_assumed fails on negative values. Is that supposed to happen? + "input_value, expected_result", + [ + ("123456-3", 0.00123456), # 1.23456e-3 + ("678901+2", 678.901), # 6.78901e2 + ("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) + + +def test_parse_tle_data(): + 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" + ) + + 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)