Add TLE.from_orbit and tle.to_lines functions

README.md

@@ -59,6 +59,11 @@ n=15.50437522, rev_num=18780)
 
 and you can then access its attributes like `t.argp`, `t.epoch`...
 
+TLE can also be written out to a string based on the orbital value:
+```python
+print(tle.to_lines())
+```
+
 ### TLE format specification
 
 Some more or less complete TLE format specifications can be found on the following websites:

tests/conftest.py

@@ -9,21 +9,53 @@ def tle_string():
 1 25544U 98067A   19249.04864348  .00001909  00000-0  40858-4 0  9990
 2 25544  51.6464 320.1755 0007999  10.9066  53.2893 15.50437522187805"""
 
+@pytest.fixture
+def tle_from_lines_to_lines_expected_string():
+    # expected result of tle = TLE.from_lines(*tle_lines).to_lines()
+    return """ISS (ZARYA)
+1 25544U 98067A   19249.04864348 +.00001909 +00000-0 +40858-4 0  9990
+2 25544  51.6464 320.1755 0007999  10.9066  53.2893 15.50437522187805"""
+
+@pytest.fixture
+def tle_from_orbit_to_lines_expected_string():
+    # expected result of tle = TLE.from_lines(*tle_lines).to_lines()
+    return """ISS (ZARYA)
+1 25544U 98067A   19249.04864348 +.00001909 +00000-0 +40858-4 0  9990
+2 25544  51.6464 320.1755 0007999  10.9066  53.2893 15.50437522187805"""
+
 @pytest.fixture
 def tle_string2():
     # This TLE tests high mean anomaly value.
     return """NOAA 18
 1 28654U 05018A   20098.54037539  .00000075  00000-0  65128-4 0  9992
 2 28654  99.0522 154.2797 0015184  73.2195 287.0641 14.12501077766909"""
 
+@pytest.fixture
+def tle_strings2_from_orbit_expected():
+    return """NOAA 18
+1 28654U 05018A   20098.54037539 +.00000075 +00000-0 +65128-4 0  9992
+2 28654  99.0522 154.2797 0015184  73.2195 287.0641 14.12501077766909"""
+
 @pytest.fixture
 def tle_lines(tle_string):
     return tle_string.splitlines()
 
+@pytest.fixture
+def tle_from_lines_to_lines_expected(tle_from_lines_to_lines_expected_string):
+    return tle_from_lines_to_lines_expected_string.splitlines()
+
+@pytest.fixture
+def tle_from_orbit_to_lines_expected(tle_from_orbit_to_lines_expected_string):
+    return tle_from_orbit_to_lines_expected_string.splitlines()
+
 @pytest.fixture
 def tle_lines2(tle_string2):
     return tle_string2.splitlines()
 
+@pytest.fixture
+def tle_lines2_from_orbit_expected(tle_strings2_from_orbit_expected):
+    return tle_strings2_from_orbit_expected.splitlines()
+
 @pytest.fixture
 def tle():
     return TLE('ISS (ZARYA)', '25544', 'U', '98067A',

tests/test_tle.py

@@ -3,14 +3,17 @@
 
 import astropy.units as u
 
+
 def test_from_lines(tle_lines):
     t = TLE.from_lines(*tle_lines)
     assert isinstance(t, TLE)
 
+
 def test_from_lines_high_M(tle_lines2):
     t = TLE.from_lines(*tle_lines2)
     assert isinstance(t, TLE)
 
+
 def test_from_lines_with_units(tle_lines):
     t = TLEu.from_lines(*tle_lines)
     assert isinstance(t, TLEu)
@@ -27,3 +30,47 @@ def test_asdict(tle):
 
 def test_astuple(tle):
     assert type(tle)(*tle.astuple()) == tle
+
+
+def test_to_lines(tle_lines, tle_from_lines_to_lines_expected):
+    t = TLE.from_lines(*tle_lines)
+    lines = t.to_lines()
+    tle_joined_lines = '\n'.join(tle_lines)
+    assert (lines == tle_from_lines_to_lines_expected)
+
+def test_orbit_to_lines(tle_lines, tle_from_orbit_to_lines_expected):
+    tle = TLE.from_lines(*tle_lines)
+    orbit = tle.to_orbit()
+    # lines = orbit.to_lines()
+
+    lines_from_orbit = TLE.from_orbit(
+        orbit,
+        name=tle.name,
+        int_desig=tle.int_desig,
+        classification=tle.classification,
+        norad=tle.norad,
+        dn_o2=tle.dn_o2,
+        ddn_o6=tle.ddn_o6,
+        bstar=tle.bstar,
+        rev_num=tle.rev_num,
+    ).to_lines()
+    assert lines_from_orbit == tle_from_orbit_to_lines_expected
+
+
+def test_from_orbit_to_lines2(tle_lines2, tle_lines2_from_orbit_expected):
+    tle = TLE.from_lines(*tle_lines2)
+    orbit = tle.to_orbit()
+    # lines = orbit.to_lines()
+
+    lines_from_orbit = TLE.from_orbit(
+        orbit,
+        name=tle.name,
+        int_desig=tle.int_desig,
+        classification=tle.classification,
+        norad=tle.norad,
+        dn_o2=tle.dn_o2,
+        ddn_o6=tle.ddn_o6,
+        bstar=tle.bstar,
+        rev_num=tle.rev_num,
+    ).to_lines()
+    assert tle_lines2_from_orbit_expected == lines_from_orbit

tletools/tle.py

@@ -31,9 +31,9 @@
 # Maybe remove them from here?
 from poliastro.twobody import Orbit as _Orbit
 from poliastro.bodies import Earth as _Earth
+from poliastro.core import angles
 
-from .utils import partition, rev as u_rev, M_to_nu as _M_to_nu
-
+from .utils import partition, rev as u_rev, M_to_nu as _M_to_nu, nu_to_M as _nu_to_M
 
 DEG2RAD = np.pi / 180.
 RAD2DEG = 180. / np.pi
@@ -69,6 +69,42 @@ def _parse_float(s):
     return float(s[0] + '.' + s[1:6] + 'e' + s[6:8])
 
 
+def _float_to_string(f: float, digits: int = 8) -> str:
+    """Convert a float to a string with implicit dot and exponential notation.
+
+    >>> _float_to_string(0.00012345, digits=5)
+    '+12345-3'
+    >>> _float_to_string(-0.00012345, digits=5)
+    '-12345-3'
+    >>> _float_to_string(0.0, digits=5)
+    '+00000-5'
+    """
+    if f == 0:
+        # zero gets a special string
+        return "+" + "0" * digits + "-0"
+    format_string = '{:+.' + str(digits) + 'E}'
+    s = format_string.format(f)
+    # skip the first zero in the exponent, and if it is +0, make it -0 to confirm to the TLE convention
+    exponent = int(s[digits + 4: digits + 7])
+    exponent = exponent + 1
+
+    # skip the decimal point, and E
+    return s[0:2] + s[3:digits + 2] + str(exponent)
+
+
+def _calculate_check_sum_on_tle_line(line: str) -> int:
+    """Calculate the checksum of a TLE line.
+
+    The checksum is calculated by taking the sum of all the digits in the line, ignoring spaces, and then taking the
+    modulo 10 of that sum.
+
+    >>> _calculate_check_sum_on_tle_line('1 25544U 98067A   19249.04864348  .00001909  00000-0  40858-4 0  9990')
+    0
+    """
+    sum_of_digits = sum(int(c) for c in line if c.isdigit()) + sum(1 for c in line if c == '-')
+    return sum_of_digits % 10
+
+
 @attr.s
 class TLE:
     """Data class representing a single TLE.
@@ -145,15 +181,15 @@ def epoch(self):
         """Epoch of the TLE, as an :class:`astropy.time.Time` object."""
         if self._epoch is None:
             year = np.datetime64(self.epoch_year - 1970, 'Y')
-            day = np.timedelta64(int((self.epoch_day - 1) * 86400 * 10**6), 'us')
+            day = np.timedelta64(int((self.epoch_day - 1) * 86400 * 10 ** 6), 'us')
             self._epoch = Time(year + day, format='datetime64', scale='utc')
         return self._epoch
 
     @property
     def a(self):
         """Semi-major axis."""
         if self._a is None:
-            self._a = (_Earth.k.value / (self.n * np.pi / 43200) ** 2) ** (1/3) / 1000
+            self._a = (_Earth.k.value / (self.n * np.pi / 43200) ** 2) ** (1 / 3) / 1000
         return self._a
 
     @property
@@ -226,6 +262,76 @@ def to_orbit(self, attractor=_Earth):
             nu=u.Quantity(self.nu, u.deg),
             epoch=self.epoch)
 
+    @classmethod
+    def from_orbit(cls, orbit: _Orbit, name: str = "UNASSIGNED", norad: str = "00000", classification: str = "U",
+                   int_desig: str = "00000A",
+                   dn_o2: float = 0.0, ddn_o6: float = 0.0, bstar: float = 0.0, set_num: int = 999, rev_num: int = 999):
+        '''Convert from a :class:`poliastro.twobody.orbit.Orbit` around the attractor into a TLE.
+        Additional information, such as the name, NORAD ID, classification, int_desig, dn_o2, ddn_o6, bstar, set_num, and rev_num needs
+        to be manually provided, or copied from a valid TLE.
+
+        >>> from poliastro.twobody import Orbit
+        >>> from astropy import units as u
+        >>> from poliastro.bodies import Earth
+        >>> from astropy.time import Time
+        >>> orbit_epoch_str = "2024-02-01T00:00:00.000000Z"
+        >>> orbit_epoch_time = Time.strptime(orbit_epoch_str, "%Y-%m-%dT%H:%M:%S.%fZ")
+        >>> orbit = Orbit.from_classical(
+        ...     attractor=Earth,
+        ...     a = 6788 << u.km,
+        ...     ecc = 0.0007999 << u.one,
+        ...     inc = 51.6464 << u.deg,
+        ...     raan= 320.1755 << u.deg,
+        ...     argp= 10.9066 << u.deg,
+        ...     nu= 53.2893 << u.deg,
+        ...     epoch = orbit_epoch_time
+        ... )
+        >>> tle = TLE.from_orbit(
+        ...     orbit,
+        ...     name="UNASSIGNED",
+        ...     norad="00000",
+        ...     classification="U",
+        ...     int_desig="00000A",
+        ...     dn_o2=0.0,
+        ...     ddn_o6=0.0,
+        ...     bstar=0.0,
+        ...     set_num=999,
+        ...     rev_num=999
+        ... )
+        >>> tle_lines = tle.to_lines()
+        >>> tle_string = "\n".join(tle_lines)
+        >>> tle_string = """UNASSIGNED
+        ... 1 00000U 00000A   24032.00000000 +.00001909 +00000-0 +40858-4 0  9999
+        ... 2 00000  51.6464 320.1755 0007999  10.9066  53.2158 15.52351307009990"""
+        '''
+
+        mean_motion = orbit.n * (24 * 60 * 60 * u.s)
+        mean_motion = mean_motion / ((2 * np.pi) << u.rad)
+        mean_motion = mean_motion.value
+        M = _nu_to_M(orbit.nu.to(u.rad).value, orbit.ecc.value) * RAD2DEG,  # mean anomaly
+        M = (M[0] + 360.0) % 360.0  # ensures the M's range is from 0 to 360, or [0, 360)
+        return cls(
+            name=name,
+            norad=norad,
+            classification=classification,
+            int_desig=int_desig,
+            epoch_year=orbit.epoch.to_datetime().year,
+            epoch_day=orbit.epoch.to_datetime().timetuple().tm_yday + orbit.epoch.to_datetime().hour / 24 + orbit.epoch.to_datetime().minute / (
+                    24 * 60) + orbit.epoch.to_datetime().second / (
+                              24 * 60 * 60) + orbit.epoch.to_datetime().microsecond / (24 * 60 * 60 * 1000000),
+            dn_o2=dn_o2,
+            ddn_o6=ddn_o6,
+            bstar=bstar,
+            set_num=set_num,
+            inc=orbit.inc.to(u.deg).value,
+            raan=orbit.raan.to(u.deg).value,
+            ecc=orbit.ecc.value,
+            argp=orbit.argp.to(u.deg).value,
+            M=M,
+            n=mean_motion,
+            rev_num=rev_num
+        )
+
     def astuple(self):
         """Return a tuple of the attributes."""
         return attr.astuple(self)
@@ -239,6 +345,30 @@ def asdict(self, computed=False, epoch=False):
             d.update(epoch=self.epoch)
         return d
 
+    def to_lines(self):
+        templates = [
+            "{name}",
+            "1 {norad}{classification} {int_desig}   {epoch_year_last_digits:02d}{epoch_day:012.8f} {dn_o2_wo_leading_zero} {ddn_o6_wo_e} {bstar_wo_e} 0  {set_num:3d}",
+            "2 {norad} {inc:8.4f} {raan:8.4f} {ecc_wo_leading_zero} {argp:8.4f} {M:8.4f} {n:11.8f}{rev_num:05d}",
+        ]
+        additional_dict = {
+            'epoch_year_last_digits': self.epoch_year % 100,
+            'dn_o2_wo_leading_zero': "{dn_o2:+.8f}".format(dn_o2=self.dn_o2).replace('0.', '.'),
+            # dn_o2 without leading zero
+            'ddn_o6_wo_e': _float_to_string(self.ddn_o6, digits=5),
+            'bstar_wo_e': _float_to_string(self.bstar, digits=5),
+            'ecc_wo_leading_zero': "{ecc:.7f}".format(ecc=self.ecc).lstrip('0').lstrip('.'),  # ecc without leading zero
+        }
+        lines = [template.format(**{**self.asdict(), **additional_dict}) for template in templates]
+        # line 1 checksum
+        line_1_mod = _calculate_check_sum_on_tle_line(lines[1])
+        # line 2 checksum
+        line_2_mod = _calculate_check_sum_on_tle_line(lines[2])
+        lines[1] = lines[1] + str(line_1_mod)
+        lines[2] = lines[2] + str(line_2_mod)
+
+        return lines
+
 
 @attr.s
 class TLEu(TLE):
@@ -256,7 +386,7 @@ class TLEu(TLE):
     def a(self):
         """Semi-major axis."""
         if self._a is None:
-            self._a = (_Earth.k.value / self.n.to_value(u.rad/u.s) ** 2) ** (1/3) * u.m
+            self._a = (_Earth.k.value / self.n.to_value(u.rad / u.s) ** 2) ** (1 / 3) * u.m
         return self._a
 
     @property
@@ -280,8 +410,8 @@ def from_lines(cls, name, line1, line2):
             int_desig=line1[9:17],
             epoch_year=line1[18:20],
             epoch_day=float(line1[20:32]),
-            dn_o2=u.Quantity(float(line1[33:43]), u_rev / u.day**2),
-            ddn_o6=u.Quantity(_parse_float(line1[44:52]), u_rev / u.day**3),
+            dn_o2=u.Quantity(float(line1[33:43]), u_rev / u.day ** 2),
+            ddn_o6=u.Quantity(_parse_float(line1[44:52]), u_rev / u.day ** 3),
             bstar=u.Quantity(_parse_float(line1[53:61]), 1 / u.earthRad),
             set_num=line1[64:68],
             inc=u.Quantity(float(line2[8:16]), u.deg),

tletools/utils.py

@@ -1,6 +1,6 @@
 import numpy as np
 import astropy.units as u
-from poliastro.core.angles import M_to_E as _M_to_E, E_to_nu as _E_to_nu
+from poliastro.core.angles import M_to_E as _M_to_E, E_to_nu as _E_to_nu, nu_to_E as _nu_to_E, E_to_M as _E_to_M
 
 #: :class:`numpy.dtype` for a date expressed as a year.
 dt_dt64_Y = np.dtype('datetime64[Y]')
@@ -17,6 +17,18 @@
 
 u.add_enabled_units(rev)
 
+def nu_to_M(nu, ecc):
+    """Mean anomaly from true anomaly.
+    :param float nu: True anomaly in radians.
+    :param float ecc: Eccentricity.
+    :returns: `M`, the mean anomaly, between -π and π radians.
+
+    **Warning**
+
+    The mean anomaly must be between -π and π radians.
+    The eccentricity must be less than 1.
+    """
+    return _E_to_M(_nu_to_E(nu, ecc), ecc)
 
 def M_to_nu(M, ecc):
     """True anomaly from mean anomaly.