First commit

dcajacob · Dec 19, 2022 · 88d92dc · 88d92dc
commit 88d92dc
Show file tree

Hide file tree

Showing 33 changed files with 166,785 additions and 0 deletions.
diff --git a/README.md b/README.md
@@ -0,0 +1,34 @@
+TLE Tailor
+
+TLE Tailor is a set of pure-python tools to generate Two Line Element Sets (TLEs) for Earth-orbiting satellites.
+
+Generating TLEs is more difficult than it seems. TLEs are mean to be used with the SGP4 orbit propagation model. SGP4 uses mean orbital elements, which are not the same as classical / keplerian / osculating orbital elements. This causes a lot of confusion for folks. There is no direct conversion from a state vector to a TLE. Rather, we must perform trajectory sampling and find a set of TLE elements that best fit the trajectory. To do this, we use non-linear leastsquares in a process called differential correction. Sounds scary right? Well, that's what this projectis meant to solve.
+
+Currently supported methods:
+
+    * Generate a TLE from ephemeris (e.g. GPS data from a satellite)
+    * Generate a TLE from a pre-launch or post-launch state vector
+    * Shift the epoch time of a TLE
+    * Match a TLE to a SpaceTrack TLE
+
+Why do I need this?
+
+    * Confidently identify and contact your satellite after launch
+
+How do I use this?
+    * Before launch or shortly after, if you're given a predicted or measured state vector by your launch provider, you can make a TLE to simplify your LEOP tracking
+    * If you make contact with your satellite and download some GPS data, you can use it to generate a TLE. This will allow you to confidently identify your satellite for 18 SDS
+    * You can continue to make supplemental TLEs using your GPS data
+
+Caveats:
+    * This is a brand new project and although it works, there may be significant changes as time goes on
+
+Future Work:
+    * Additional ephemeris type supported, e.g. Range, Rate, Az, El, etc.
+      * This could potentially use measurements from operators' ground stations / strf https://github.com/cbassa/strf to generate a TLE
+
+Notes:
+    * Differential correction usually uses finite differencing to calculate the Jacobian of complex or black-box models. We might call SGP4 a gre-box, since we have the code (although the AFSPC version may be different). Finite differencing is a numerical method to dodge this complication and works great. But I have also implemented a neat way to automatically derive the analytical Jacobian using JAX's autograd capability. This can significantly speed up the calculations, although it incurs an upfront JIT compiling penalty, so it's best for cases where you're processing a lot of data or TLEs.
+
+References:
+    * [Mention Vallado and his previous work]
diff --git a/TLEFit - COE - FD.ipynb b/TLEFit - COE - FD.ipynb
diff --git a/TLEFit - COE - JAX.ipynb b/TLEFit - COE - JAX.ipynb
diff --git a/TLEFit - EQN - FD.ipynb b/TLEFit - EQN - FD.ipynb
diff --git a/TLEFit - EQN - GPS FD - GPS Sat.ipynb b/TLEFit - EQN - GPS FD - GPS Sat.ipynb
diff --git a/TLEFit - EQN - GPS FD - Icesat.ipynb b/TLEFit - EQN - GPS FD - Icesat.ipynb
diff --git a/TLEFit - EQN - GPS JAX - GPS Sat.ipynb b/TLEFit - EQN - GPS JAX - GPS Sat.ipynb
diff --git a/TLEFit - EQN - GPS JAX - Icesat.ipynb b/TLEFit - EQN - GPS JAX - Icesat.ipynb
diff --git a/TLEFit - EQN - JAX.ipynb b/TLEFit - EQN - JAX.ipynb
diff --git a/coarse_fit.py b/coarse_fit.py
@@ -0,0 +1,151 @@
+import numpy as np
+
+from skyfield.api import EarthSatellite, load
+from sgp4.api import Satrec, WGS72
+from sgp4.model import wgs72, wgs84
+from sgp4.ext import rv2coe
+
+ts = load.timescale()
+
+
+def coarse_fit(satellite, rms_epsilon=0.002, debug=False):
+    """A hack to get a first-cut estimate of TLE mean elements from keplerian elements
+
+    Based on some code I found on the internet. FIXME: Find it again and cite properly.
+    """
+
+    # Form initial state estimate
+    _, r, v = satellite.model.sgp4_tsince(0)
+
+    for ix in range(15):
+        coe_nom = rv2coe(r, v, wgs72.mu)
+        p, a, ecc, incl, omega, argp, nu, m, arglat, truelon, lonper = coe_nom
+        n = np.sqrt(wgs72.mu / a**3) * 60  # radians / min
+        bstar = 0
+
+        elements = [n, a, ecc, incl, omega, argp, m, bstar]
+
+        a_nom = np.cbrt(wgs72.mu / (satellite.model.no_kozai / 60) ** 2)
+        orig_elements = [
+            satellite.model.no_kozai,
+            a_nom,
+            satellite.model.ecco,
+            satellite.model.inclo,
+            satellite.model.nodeo,
+            satellite.model.argpo,
+            satellite.model.mo,
+            satellite.model.bstar,
+        ]
+
+        if debug and ix == 0:
+            print(
+                "Guess ",
+                elements[0],
+                elements[1],
+                elements[2],
+                np.degrees(elements[3]),
+                np.degrees(elements[4]),
+                np.degrees(elements[5]),
+                np.degrees(elements[6]),
+                elements[7],
+            )
+
+        # Print the initial difference
+        satrec = Satrec()
+        satrec.sgp4init(
+            WGS72,
+            "i",
+            satellite.model.satnum,
+            round(
+                satellite.model.jdsatepoch + satellite.model.jdsatepochF - 2433281.5, 8
+            ),
+            bstar,
+            0.0,
+            0.0,
+            ecc,
+            argp,
+            incl,
+            m,
+            n,
+            omega,
+        )
+
+        pert_sat = EarthSatellite.from_satrec(satrec, ts)
+        pert_sat.model.jdsatepochF = satellite.model.jdsatepochF
+
+        res = np.array(pert_sat.model.sgp4_tsince(0)[1:]) - np.array(
+            satellite.model.sgp4_tsince(0)[1:]
+        )
+        b = np.concatenate((res[0], res[1]))
+        b_epoch = b
+
+        if debug:
+            print(np.linalg.norm(b[0:3]), np.linalg.norm(b[3:6]))
+
+        sigma_new = np.sqrt(b.T @ b)
+
+        #         convergence_test = np.abs((sigma_old - sigma_new) / sigma_old)
+
+        #         if convergence_test <= rms_epsilon and  sigmasigma_new <= rms_epsilon:
+        if sigma_new <= rms_epsilon:
+            break
+
+        r -= b[0:3]
+        v -= b[3:6]
+        sigma_old = sigma_new
+
+    coe_nom = rv2coe(r, v, wgs72.mu)
+    p, a, ecc, incl, omega, argp, nu, m, arglat, truelon, lonper = coe_nom
+    n = np.sqrt(wgs72.mu / a**3) * 60  # radians / min
+    elements = [n, a, ecc, incl, omega, argp, m, bstar]
+
+    if debug:
+        print(
+            "Solution ",
+            elements[0],
+            elements[1],
+            elements[2],
+            np.degrees(elements[3]),
+            np.degrees(elements[4]),
+            np.degrees(elements[5]),
+            np.degrees(elements[6]),
+            elements[7],
+        )
+        print(
+            "Original ",
+            orig_elements[0],
+            orig_elements[1],
+            orig_elements[2],
+            np.degrees(orig_elements[3]),
+            np.degrees(orig_elements[4]),
+            np.degrees(orig_elements[5]),
+            np.degrees(orig_elements[6]),
+            orig_elements[7],
+        )
+
+    elements = elements[1:]
+
+    a, ecc, incl, omega, argp, m, bstar = elements
+
+    satrec = Satrec()
+    satrec.sgp4init(
+        WGS72,
+        "i",
+        satellite.model.satnum,
+        round(satellite.model.jdsatepoch + satellite.model.jdsatepochF - 2433281.5, 8),
+        bstar,
+        0.0,
+        0.0,
+        ecc,
+        argp,
+        incl,
+        m,
+        n,
+        omega,
+    )
+
+    close_sat = EarthSatellite.from_satrec(satrec, ts)
+    close_sat.model.jdsatepochF = satellite.model.jdsatepochF
+
+    # return elements
+    return close_sat