From 85662c64507c065fc2b2f6175960f6a63055e0aa Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Fri, 10 Apr 2026 22:28:27 -0700 Subject: [PATCH 01/19] setup for tle, sgp4 adcs position determination --- src/lib/tle.py | 148 +++++++++++++++++++++++++++++++++++++++++ src/tasks/adcs/sgp4.py | 12 ++++ 2 files changed, 160 insertions(+) create mode 100644 src/lib/tle.py create mode 100644 src/tasks/adcs/sgp4.py diff --git a/src/lib/tle.py b/src/lib/tle.py new file mode 100644 index 0000000..e38aa19 --- /dev/null +++ b/src/lib/tle.py @@ -0,0 +1,148 @@ +""" +Functions and Variables used by ADCS to update and use TLE (two-line element data) + +Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, +""" + +def _conv_year(s): + """Interpret a two-digit year string.""" + if isinstance(s, int): + return s + y = int(s) + return y + (1900 if y >= 57 else 2000) + +def _parse_decimal(s): + """Parse a floating point with implicit leading dot. + + >>> _parse_decimal('378') + 0.378 + """ + return float('.' + s) + +def _parse_float(s): + """Parse a floating point with implicit dot and exponential notation. + + >>> _parse_float(' 12345-3') + 0.00012345 + >>> _parse_float('+12345-3') + 0.00012345 + >>> _parse_float('-12345-3') + -0.00012345 + """ + return float(s[0] + '.' + s[1:6] + 'e' + s[6:8]) + +class TLE: + """ + Two line-elements (TLEs) are unpacked from both given and propagated data. + This implementation uses Keplerian orbital parameters + + A two-line element set (TLE) is a data format encoding a list of orbital + elements of an Earth-orbiting object for a given point in time, the epoch. + + All the attributes parsed from the TLE are expressed in the same units that + are used in the TLE format. + + :str name: + Name of the satellite. + :str norad: + NORAD catalog number (https://en.wikipedia.org/wiki/Satellite_Catalog_Number). + :str classification: + 'U', 'C', 'S' for unclassified, classified, secret. + :str int_desig: + International designator (https://en.wikipedia.org/wiki/International_Designator), + :int epoch_year: + Year of the epoch. + :float epoch_day: + Day of the year plus fraction of the day. + :float dn_o2: + First time derivative of the mean motion divided by 2. + :float ddn_o6: + Second time derivative of the mean motion divided by 6. + :float bstar: + BSTAR coefficient (https://en.wikipedia.org/wiki/BSTAR). + :int set_num: + Element set number. + :float inc: + Inclination. + :float raan: + Right ascension of the ascending node. + :float ecc: + Eccentricity. + :float argp: + Argument of perigee. + :float M: + Mean anomaly. + :float n: + Mean motion. + :int rev_num: + Revolution number. + """ + + def __init__(self, + # ID parameters, Line 1 + name:str, norad:str, classification:str, int_desig:str, + # time (derivative) parameters, line 1 + epoch_year:int, epoch_day:float, dn_o2:float, ddn_o6:float, bstar:float, set_num:int, + # keplerian parameters, line 2 + inc:float, raan:float, ecc:float, argp:float, M:float, n:float, rev_num:int): + # Oh my lord prepare for absolute misery on earth + + self.name = str.strip(name) + self.norad = str.strip(norad) + self.classification = classification + self.int_desig = str.strip(int_desig) + + self.epoch_year = _conv_year(epoch_year) + self.epoch_day = epoch_day + self.dn_o2 = dn_o2 + self.ddn_o6 = ddn_o6 + self.bstar = bstar + self.set_num = int(set_num) + + self.inc = inc + self.raan = raan + self.ecc = ecc + self.argp = argp + self.M = M + self.n = n + self.rev_num = int(rev_num) + + @classmethod + def from_lines(cls, name, line1, line2): + """Parse a TLE from its constituent lines. + + All the attributes parsed from the TLE are expressed in the same units that + are used in the TLE format. + """ + return cls( + name=name, + norad=line1[2:7], + classification=line1[7], + int_desig=line1[9:17], + epoch_year=line1[18:20], + epoch_day=float(line1[20:32]), + dn_o2=float(line1[33:43]), + ddn_o6=_parse_float(line1[44:52]), + bstar=_parse_float(line1[53:61]), + set_num=line1[64:68], + inc=float(line2[8:16]), + raan=float(line2[17:25]), + ecc=_parse_decimal(line2[26:33]), + argp=float(line2[34:42]), + M=float(line2[43:51]), + n=float(line2[52:63]), + rev_num=line2[63:68]) + + @classmethod + def from_file(cls, filename): + """Load TLE from a file.""" + if isinstance(filename, str): + with open(filename) as fp: + return [cls.from_lines(*fp.readlines[:2])] + + @classmethod + def from_str(cls, string): + """Load TLE from a string.""" + return [cls.from_lines(*string.split('\n')[:2])] + + diff --git a/src/tasks/adcs/sgp4.py b/src/tasks/adcs/sgp4.py new file mode 100644 index 0000000..54ecb3d --- /dev/null +++ b/src/tasks/adcs/sgp4.py @@ -0,0 +1,12 @@ +""" + +""" + +import math + +try: + import ulab.numpy as np # For CircuitPython +except ImportError: + import numpy as np # For GitHub Actions / PC testing + + From 407d45fd1687ff04f03400c4181f17a9ccc38c3c Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Sun, 19 Apr 2026 18:33:28 -0700 Subject: [PATCH 02/19] code for propagation function, beginning some reorganisation --- src/lib/datastores/adcs.py | 40 +++++-- src/lib/tle.py | 35 +++++- src/tasks/adcs/sgp4.py | 229 ++++++++++++++++++++++++++++++++++++- src/tasks/adcs/triad.py | 1 - 4 files changed, 294 insertions(+), 11 deletions(-) diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 9c9042c..6d216f7 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -4,6 +4,8 @@ are set to `None` throughout this module. """ +import tle + class Datastore: """ Datastore class for adcs processes. Holds time, sensor, and attitude data to be used system-wide @@ -33,7 +35,7 @@ def __init__(self): None # Quaternion representing attitude from body frame to inertial frame ) self.mode = self.DETUMBLE - self.tle: TLE = TLE() + self.tle: tle.TLE = tle.TLE() class AdcsTime: """ @@ -54,11 +56,35 @@ def __init__(self): self.magnetometer = None self.gyroscope = None -class TLE: +class Satrec: """ - Attitude helper class + Parameters, constants that are commonly used across the sgp4 logical flow + + Usually built from TLE """ - def __init__(self): - # reference vectors in inertial frame - self.ref_vec1 = 0.0 # more accurate vector - self.ref_vec2 = 0.0 # less accurate vector + + """ + t - time since + mo, mdot - mean anomaly + argpo, argpdot - argument of perigee + nodeo, nodedot, nodecf - RAAN value, drift, and correction respectively + bstar - atmospheric drag + cc1, cc4, cc5 - drag coefficient terms + ecco, inclo - eccentricity, inclination (no derivatives) + mm, nm - mean motion before and after corrections + error, error_message - might want to leave this up to cdh but implement + this here? + """ + + def __init__(self, t, mo, # mean anomaly + argpo, # argument of perigee + nodeo, # RAAN + ): + pass + + @classmethod + def from_tle_array(cls): + + obj = cls() + + return obj diff --git a/src/lib/tle.py b/src/lib/tle.py index e38aa19..626338a 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -4,6 +4,8 @@ Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, """ +from datastores.adcs import Satrec + def _conv_year(s): """Interpret a two-digit year string.""" if isinstance(s, int): @@ -78,9 +80,9 @@ class TLE: Revolution number. """ - def __init__(self, + def __init__(self, name:str, # ID parameters, Line 1 - name:str, norad:str, classification:str, int_desig:str, + norad:str, classification:str, int_desig:str, # time (derivative) parameters, line 1 epoch_year:int, epoch_day:float, dn_o2:float, ddn_o6:float, bstar:float, set_num:int, # keplerian parameters, line 2 @@ -88,6 +90,7 @@ def __init__(self, # Oh my lord prepare for absolute misery on earth self.name = str.strip(name) + self.norad = str.strip(norad) self.classification = classification self.int_desig = str.strip(int_desig) @@ -145,4 +148,32 @@ def from_str(cls, string): """Load TLE from a string.""" return [cls.from_lines(*string.split('\n')[:2])] + def to_array(self): + """ + Return 2D array of TLE values + + Indexed as + [line, col] + + n is mean motion, d suggests time derivative + + name: [0,0] + + norad: [1,0] classification: [1,1] int_desig: [1,2] epoch_year: [1,3] day: [1,4] + dn/2: [1,5] ddn/6: [1,6] bstar: [1,7] set_num: [1,8] + + inclination: [2,0] RAAN: [2,1] eccentricity: [2,2] arg_perigee: [2,3] Mean Anomaly: [2,4] n: [2,5] rev_num: [2,6] + """ + TLE() + return [ + [self.name], # Line 0 + [self.norad, self.classification, self.int_desig, # line 1 ID + self.epoch_year, self.epoch_day, self.dn_o2, self.ddn_o6, self.bstar, self.set_num], # line 1 time-derivative + [self.inc, self.raan, self.ecc, self.argp, self.M, self.n, self.rev_num] # line 2 orbital params + ] + + def to_sgp4_params(self): + """Return a formatted Satrec object that can immediately used in SGP4""" + return Satrec.from_tle_array(self.to_array()) + \ No newline at end of file diff --git a/src/tasks/adcs/sgp4.py b/src/tasks/adcs/sgp4.py index 54ecb3d..d5ac928 100644 --- a/src/tasks/adcs/sgp4.py +++ b/src/tasks/adcs/sgp4.py @@ -1,12 +1,239 @@ """ +Simplified General Perturbations Model 4 Implementation for Orbit Propagation +As outlined in https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf """ -import math +from tle import Satrec + +# Error codes +ECCENTRICITY = 1 # eccentricity is not within 0-1 +MOTION = 2 # error in propagating mean motion +SEMIRECT = 4 # apoapsis, periapsis characteristics error +DECAY = 6 # orbit has decayed try: import ulab.numpy as np # For CircuitPython except ImportError: import numpy as np # For GitHub Actions / PC testing +def sgp4_update(satrec: Satrec, tsince): + """ + Transforms to cartesian (r, v) from a Satrec (TLE 7+) input, then performs matrix operations + to propagate the system + + Current implementation uses an analytical approach to compute the Jacobian A matrix + Possible alternative is numerical stepping (for lower accuracy) of associated partials + """ + # -- Set mathematical constants + x2o3 = 2.0 / 3.0; + tau = 2.0 * np.pi + vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; + + satrec.t = tsince + + # -- Update for secular gravity and atmospheric drag + xmdf = satrec.mo + satrec.mdot * satrec.t; + argpdf = satrec.argpo + satrec.argpdot * satrec.t; + nodedf = satrec.nodeo + satrec.nodedot * satrec.t; + argpm = argpdf + mm = xmdf + t2 = satrec.t * satrec.t; + nodem = nodedf + satrec.nodecf * t2; + tempa = 1.0 - satrec.cc1 * satrec.t; + tempe = satrec.bstar * satrec.cc4 * satrec.t; + templ = satrec.t2cof * t2; + + # -- Extra mean quantities + + # -- Add lunar-solar periodics + + # -- Long period periodics + + # -- Solve Kepler's + + # -- Short period + # Preliminary quantities + + # - Update for short period periodics + + # - Orientation vectors + + # - Compute r and v + if satrec.isimp != 1: + + delomg = satrec.omgcof * satrec.t; + # sgp4fix use mutliply for speed instead of pow + delmtemp = 1.0 + satrec.eta * np.cos(xmdf); + delm = satrec.xmcof * \ + (delmtemp * delmtemp * delmtemp - + satrec.delmo); + temp = delomg + delm; + mm = xmdf + temp; + argpm = argpdf - temp; + t3 = t2 * satrec.t; + t4 = t3 * satrec.t; + tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ + satrec.d4 * t4; + tempe = tempe + satrec.bstar * satrec.cc5 * (np.sin(mm) - + satrec.sinmao); + templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + + satrec.t * satrec.t5cof); + + nm = satrec.no_unkozai; + em = satrec.ecco; + inclm = satrec.inclo; + + if nm <= 0.0: + + satrec.error_message = ('mean motion {0:f} is less than zero' + .format(nm)) + satrec.error = MOTION + # sgp4fix add return + return False, False; + + am = pow((satrec.xke / nm),x2o3) * tempa * tempa; + nm = satrec.xke / pow(am, 1.5); + em = em - tempe; + + # fix tolerance for error recognition + # sgp4fix am is fixed from the previous nm check + if em >= 1.0 or em < -0.001: # || (am < 0.95) + + satrec.error_message = ('mean eccentricity {0:f} not within' + ' range 0.0 <= e < 1.0'.format(em)) + satrec.error = ECCENTRICITY + # sgp4fix to return if there is an error in eccentricity + return False, False; + + # sgp4fix fix tolerance to avoid a divide by zero + if em < 1.0e-6: + em = 1.0e-6; + mm = mm + satrec.no_unkozai * templ; + xlm = mm + argpm + nodem; + emsq = em * em; + temp = 1.0 - emsq; + + nodem = nodem % tau if nodem >= 0.0 else -(-nodem % tau) + argpm = argpm % tau + xlm = xlm % tau + mm = (xlm - argpm - nodem) % tau + + # sgp4fix recover singly averaged mean elements + satrec.am = am; + satrec.em = em; + satrec.im = inclm; + satrec.Om = nodem; + satrec.om = argpm; + satrec.mm = mm; + satrec.nm = nm; + + # ----------------- compute extra mean quantities ------------- + sinim = np.sin(inclm); + cosim = np.cos(inclm); + + # -------------------- add lunar-solar periodics -------------- + ep = em; + xincp = inclm; + argpp = argpm; + nodep = nodem; + mp = mm; + sinip = sinim; + cosip = cosim; + + # -------------------- long period periodics ------------------ + axnl = ep * np.cos(argpp); + temp = 1.0 / (am * (1.0 - ep * ep)); + aynl = ep* np.sin(argpp) + temp * satrec.aycof; + xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; + + # --------------------- solve kepler's equation --------------- + u = (xl - nodep) % tau + eo1 = u; + tem5 = 9999.9; + ktr = 1; + # sgp4fix for kepler iteration + # the following iteration needs better limits on corrections + while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: + + sineo1 = np.sin(eo1); + coseo1 = np.cos(eo1); + tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; + tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; + if np.fabs(tem5) >= 0.95: + tem5 = 0.95 if tem5 > 0.0 else -0.95; + eo1 = eo1 + tem5; + ktr = ktr + 1; + + # ------------- short period preliminary quantities ----------- + ecose = axnl*coseo1 + aynl*sineo1; + esine = axnl*sineo1 - aynl*coseo1; + el2 = axnl*axnl + aynl*aynl; + pl = am*(1.0-el2); + if pl < 0.0: + + satrec.error_message = ('semilatus rectum {0:f} is less than zero' + .format(pl)) + satrec.error = SEMIRECT + # sgp4fix add return + return False, False; + + else: + + rl = am * (1.0 - ecose); + rdotl = np.sqrt(am) * esine/rl; + rvdotl = np.sqrt(pl) / rl; + betal = np.sqrt(1.0 - el2); + temp = esine / (1.0 + betal); + sinu = am / rl * (sineo1 - aynl - axnl * temp); + cosu = am / rl * (coseo1 - axnl + aynl * temp); + su = np.atan2(sinu, cosu); + sin2u = (cosu + cosu) * sinu; + cos2u = 1.0 - 2.0 * sinu * sinu; + temp = 1.0 / pl; + temp1 = 0.5 * satrec.j2 * temp; + temp2 = temp1 * temp; + + # -------------- update for short period periodics ------------ + mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ + 0.5 * temp1 * satrec.x1mth2 * cos2u; + su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; + xnode = nodep + 1.5 * temp2 * cosip * sin2u; + xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; + mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke; + rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + + 1.5 * satrec.con41) / satrec.xke; + + # --------------------- orientation vectors ------------------- + sinsu = np.sin(su); + cossu = np.cos(su); + snod = np.sin(xnode); + cnod = np.cos(xnode); + sini = np.sin(xinc); + cosi = np.cos(xinc); + xmx = -snod * cosi; + xmy = cnod * cosi; + ux = xmx * sinsu + cnod * cossu; + uy = xmy * sinsu + snod * cossu; + uz = sini * sinsu; + vx = xmx * cossu - cnod * sinsu; + vy = xmy * cossu - snod * sinsu; + vz = sini * cossu; + + # --------- position and velocity (in km and km/sec) ---------- + _mr = mrt * satrec.radiusearthkm + r = (_mr * ux, _mr * uy, _mr * uz) + v = ((mvt * ux + rvdot * vx) * vkmpersec, + (mvt * uy + rvdot * vy) * vkmpersec, + (mvt * uz + rvdot * vz) * vkmpersec) + + # sgp4fix for decaying satellites + if mrt < 1.0: + + satrec.error_message = ('mrt {0:f} is less than 1.0 indicating' + ' the satellite has decayed'.format(mrt)) + satrec.error = DECAY + + return r, v + diff --git a/src/tasks/adcs/triad.py b/src/tasks/adcs/triad.py index bf6b3d1..6e8728e 100644 --- a/src/tasks/adcs/triad.py +++ b/src/tasks/adcs/triad.py @@ -12,7 +12,6 @@ SINGULAR = 3 # Singular: failure from insufficient information to estimate attitude NORM_ERR = 4 # Normalization error: failure from prevented division by zero - # pylint: disable=too-many-locals def triad_algorithm( r1: np.ndarray, r2: np.ndarray, b1: np.ndarray, b2: np.ndarray From c1c3b90569d150a95a2d44bf34b44917e04a2769 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Fri, 24 Apr 2026 16:06:25 -0700 Subject: [PATCH 03/19] sgp4 prelim implementation --- src/lib/datastores/adcs.py | 33 ----- src/lib/tle.py | 270 ++++++++++++++++++++++++++++++++++++- src/tasks/adcs/sgp4.py | 85 ++++-------- 3 files changed, 294 insertions(+), 94 deletions(-) diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 6d216f7..7ae8d3c 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -55,36 +55,3 @@ def __init__(self): self.sun = None self.magnetometer = None self.gyroscope = None - -class Satrec: - """ - Parameters, constants that are commonly used across the sgp4 logical flow - - Usually built from TLE - """ - - """ - t - time since - mo, mdot - mean anomaly - argpo, argpdot - argument of perigee - nodeo, nodedot, nodecf - RAAN value, drift, and correction respectively - bstar - atmospheric drag - cc1, cc4, cc5 - drag coefficient terms - ecco, inclo - eccentricity, inclination (no derivatives) - mm, nm - mean motion before and after corrections - error, error_message - might want to leave this up to cdh but implement - this here? - """ - - def __init__(self, t, mo, # mean anomaly - argpo, # argument of perigee - nodeo, # RAAN - ): - pass - - @classmethod - def from_tle_array(cls): - - obj = cls() - - return obj diff --git a/src/lib/tle.py b/src/lib/tle.py index 626338a..cbcad74 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -4,7 +4,275 @@ Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, """ -from datastores.adcs import Satrec +try: + import ulab.numpy as np # For CircuitPython +except ImportError: + import numpy as np # For GitHub Actions / PC testing + +def _gstime(jdut1): + deg2rad = np.pi / 180.0 + tau = np.pi*2 + tut1 = (jdut1 - 2451545.0) / 36525.0 + temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \ + (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841 # sec + temp = (temp * deg2rad / 240.0) % tau # 360/86400 = 1/240, to deg, to rad + + # ------------------------ check quadrants --------------------- + if temp < 0.0: + temp += tau + + return temp + +def _initl(xke, j2, + ecco, epoch, inclo, no): + # ----------------------- earth constants ---------------------- + # sgp4fix identify constants and allow alternate values + # only xke and j2 are used here so pass them in directly + # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst + x2o3 = 2.0 / 3.0; + + # ------------- calculate auxillary epoch quantities ---------- + eccsq = ecco * ecco; + omeosq = 1.0 - eccsq; + rteosq = np.sqrt(omeosq); + cosio = np.cos(inclo); + cosio2 = cosio * cosio; + + # ------------------ un-kozai the mean motion ----------------- + ak = pow(xke / no, x2o3); + d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); + del_ = d1 / (ak * ak); + adel = ak * (1.0 - del_ * del_ - del_ * + (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); + del_ = d1/(adel * adel) + no = no / (1.0 + del_) + + ao = pow(xke / no, x2o3) + sinio = np.sin(inclo) + po = ao * omeosq + con42 = 1.0 - 5.0 * cosio2 + con41 = -con42-cosio2-cosio2 + posq = po * po + rp = ao * (1.0 - ecco) + + gsto = _gstime(epoch + 2433281.5) + + return ( + no, + ao, con41, con42, cosio, + cosio2, omeosq, posq, + rp, rteosq,sinio , gsto + ) + +class Satrec: + """ + Parameters, constants that are commonly used across the sgp4 logical flow + + Usually built from TLE + """ + + # Error codes + ECCENTRICITY = 1 # eccentricity is not within 0-1 + MOTION = 2 # error in propagating mean motion + SEMIRECT = 4 # apoapsis, periapsis characteristics error + DECAY = 6 # orbit has decayed + + """ + t - time since + mo, mdot - mean anomaly + argpo, argpdot - argument of perigee + nodeo, nodedot, nodecf - RAAN value, drift, and correction respectively + bstar - atmospheric drag + cc1, cc4, cc5 - drag coefficient terms + ecco, inclo - eccentricity, inclination (no derivatives) + mm, nm - mean motion before and after corrections + error, error_message - might want to leave this up to cdh but implement + this here? + isimp - using simplified model? + + """ + + """ Error Messages + MOTION + satrec.error_message = ('mean motion {0:f} is less than zero' + .format(nm)) + + ECCENTRICITY + satrec.error_message = ('mean eccentricity {0:f} not within' + ' range 0.0 <= e < 1.0'.format(em)) + + SEMIRECT + satrec.error_message = ('semilatus rectum {0:f} is less than zero' + .format(pl)) + + DECAY + satrec.error_message = ('mrt {0:f} is less than 1.0 indicating' + ' the satellite has decayed'.format(mrt)) + """ + + def __init__(self, satn, epoch, + bstar, ndot, nddot, ecco, argpo, + inclo, mo, no_kozai, + nodeo + ): + + temp4 = 1.5e-12 + + # Near Earth Variables + self.isimp = 0; self.aycof = 0.0 + self.con41 = 0.0; self.cc1 = 0.0; self.cc4 = 0.0 + self.cc5 = 0.0; self.d2 = 0.0; self.d3 = 0.0 + self.d4 = 0.0; self.delmo = 0.0; self.eta = 0.0 + self.argpdot = 0.0; self.omgcof = 0.0; self.sinmao = 0.0 + self.t = 0.0; self.t2cof = 0.0; self.t3cof = 0.0 + self.t4cof = 0.0; self.t5cof = 0.0; self.x1mth2 = 0.0 + self.x7thm1 = 0.0; self.mdot = 0.0; self.nodedot = 0.0 + self.xlcof = 0.0; self.xmcof = 0.0; self.nodecf = 0.0 + + # Earth Constants + self.mu = 398600.5; # in km3 / s2 + self.radiusearthkm = 6378.137 # km + self.xke = 0.07436685317 + self.tumin = 13.44685108 + self.j2 = 0.00108262998905 + self.j3 = -0.00000253215306 + self.j4 = -0.00000161098761 + self.j3oj2 = -0.002338890559 + + ss = 1.012229276 + qzms2ttemp = 0.00658499496 + qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; + x2o3 = 2.0 / 3.0 + + # -- initialisation markers + self.init = 'y' + self.t = 0.0 + + # -- + self.satnum_str = satn + self.classification = 'U' + + # -- + self.bstar = bstar + self.ndot = ndot + self.nddot = nddot + self.ecco = ecco + self.argpo = argpo + self.inclo = inclo + self.mo = mo + self.nodeo = nodeo + self.no_kozai = no_kozai + + # single averaged mean elements + self.am = 0.0 + self.em = 0.0 + self.im = 0.0 + self.Om = 0.0 + self.mm = 0.0 + self.nm = 0.0 + + self.error = 0 + + # -- + ( + self.no_unkozai, + ao, self.con41, con42, cosio, + cosio2, omeosq, posq, + rp, rteosq,sinio , self.gsto, + ) = _initl( + self.xke, self.j2, self.ecco, epoch, self.inclo, self.no_kozai + ) + self.a = pow( self.no_unkozai*self.tumin , (-2.0/3.0) ); + self.alta = self.a*(1.0 + self.ecco) - 1.0; + self.altp = self.a*(1.0 - self.ecco) - 1.0; + + if omeosq >= 0.0 or self.no_unkozai >= 0.0: + self.isimp = 0 + if rp < 220.0 / self.radiusearthkm + 1.0: + self.isimp = 1 + sfour = ss + qzms24 = qzms2t + pinvsq = 1.0 / posq; + + tsi = 1.0 / (ao - sfour); + self.eta = ao * self.ecco * tsi; + etasq = self.eta * self.eta; + eeta = self.ecco * self.eta; + psisq = np.fabs(1.0 - etasq); + coef = qzms24 * pow(tsi, 4.0); + coef1 = coef / pow(psisq, 3.5); + cc2 = coef1 * self.no_unkozai * (ao * (1.0 + 1.5 * etasq + eeta * + (4.0 + etasq)) + 0.375 * self.j2 * tsi / psisq * self.con41 * + (8.0 + 3.0 * etasq * (8.0 + etasq))); + self.cc1 = self.bstar * cc2; + cc3 = 0.0; + if self.ecco > 1.0e-4: + cc3 = -2.0 * coef * tsi * self.j3oj2 * self.no_unkozai * sinio / self.ecco; + self.x1mth2 = 1.0 - cosio2; + self.cc4 = 2.0* self.no_unkozai * coef1 * ao * omeosq * \ + (self.eta * (2.0 + 0.5 * etasq) + self.ecco * + (0.5 + 2.0 * etasq) - self.j2 * tsi / (ao * psisq) * + (-3.0 * self.con41 * (1.0 - 2.0 * eeta + etasq * + (1.5 - 0.5 * eeta)) + 0.75 * self.x1mth2 * + (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * self.argpo))); + self.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * + (etasq + eeta) + eeta * etasq); + cosio4 = cosio2 * cosio2; + temp1 = 1.5 * self.j2 * pinvsq * self.no_unkozai; + temp2 = 0.5 * temp1 * self.j2 * pinvsq; + temp3 = -0.46875 * self.j4 * pinvsq * pinvsq * self.no_unkozai; + self.mdot = self.no_unkozai + 0.5 * temp1 * rteosq * self.con41 + 0.0625 * \ + temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); + self.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * + (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + + temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)); + xhdot1 = -temp1 * cosio; + self.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; + self.omgcof = self.bstar * cc3 * np.cos(self.argpo); + self.xmcof = 0.0; + if self.ecco > 1.0e-4: + self.xmcof = -x2o3 * coef * self.bstar / eeta; + self.nodecf = 3.5 * omeosq * xhdot1 * self.cc1; + self.t2cof = 1.5 * self.cc1; + + if np.fabs(cosio+1.0) > 1.5e-12: + self.xlcof = -0.25 * self.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); + else: + self.xlcof = -0.25 * self.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; + self.aycof = -0.5 * self.j3oj2 * sinio; + + delmotemp = 1.0 + self.eta * np.cos(self.mo); + self.delmo = delmotemp * delmotemp * delmotemp; + self.sinmao = np.sin(self.mo); + self.x7thm1 = 7.0 * cosio2 - 1.0; + + if self.isimp != 1: + cc1sq = self.cc1 * self.cc1; + self.d2 = 4.0 * ao * tsi * cc1sq; + temp = self.d2 * tsi * self.cc1 / 3.0; + self.d3 = (17.0 * ao + sfour) * temp; + self.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ + self.cc1; + self.t3cof = self.d2 + 2.0 * cc1sq; + self.t4cof = 0.25 * (3.0 * self.d3 + self.cc1 * + (12.0 * self.d2 + 10.0 * cc1sq)); + self.t5cof = 0.2 * (3.0 * self.d4 + + 12.0 * self.cc1 * self.d3 + + 6.0 * self.d2 * self.d2 + + 15.0 * cc1sq * (2.0 * self.d2 + cc1sq)); + + # need to propagate to epoch 0.0 to really instantiate everything else before + # init flag is set to 'n' + + return True + + @classmethod + def from_tle_array(cls): + + obj = cls() + + return obj def _conv_year(s): """Interpret a two-digit year string.""" diff --git a/src/tasks/adcs/sgp4.py b/src/tasks/adcs/sgp4.py index d5ac928..606f19e 100644 --- a/src/tasks/adcs/sgp4.py +++ b/src/tasks/adcs/sgp4.py @@ -4,7 +4,7 @@ As outlined in https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf """ -from tle import Satrec +from datastores.adcs import Satrec # Error codes ECCENTRICITY = 1 # eccentricity is not within 0-1 @@ -26,7 +26,7 @@ def sgp4_update(satrec: Satrec, tsince): Possible alternative is numerical stepping (for lower accuracy) of associated partials """ # -- Set mathematical constants - x2o3 = 2.0 / 3.0; + x2o3 = 2.0 / 3.0 tau = 2.0 * np.pi vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; @@ -44,22 +44,6 @@ def sgp4_update(satrec: Satrec, tsince): tempe = satrec.bstar * satrec.cc4 * satrec.t; templ = satrec.t2cof * t2; - # -- Extra mean quantities - - # -- Add lunar-solar periodics - - # -- Long period periodics - - # -- Solve Kepler's - - # -- Short period - # Preliminary quantities - - # - Update for short period periodics - - # - Orientation vectors - - # - Compute r and v if satrec.isimp != 1: delomg = satrec.omgcof * satrec.t; @@ -85,28 +69,18 @@ def sgp4_update(satrec: Satrec, tsince): inclm = satrec.inclo; if nm <= 0.0: - - satrec.error_message = ('mean motion {0:f} is less than zero' - .format(nm)) - satrec.error = MOTION - # sgp4fix add return + satrec.error = satrec.MOTION return False, False; am = pow((satrec.xke / nm),x2o3) * tempa * tempa; nm = satrec.xke / pow(am, 1.5); em = em - tempe; - # fix tolerance for error recognition - # sgp4fix am is fixed from the previous nm check - if em >= 1.0 or em < -0.001: # || (am < 0.95) - - satrec.error_message = ('mean eccentricity {0:f} not within' - ' range 0.0 <= e < 1.0'.format(em)) - satrec.error = ECCENTRICITY - # sgp4fix to return if there is an error in eccentricity + if em >= 1.0 or em < -0.001: + satrec.error = satrec.ECCENTRICITY + return False, False; - # sgp4fix fix tolerance to avoid a divide by zero if em < 1.0e-6: em = 1.0e-6; mm = mm + satrec.no_unkozai * templ; @@ -119,7 +93,6 @@ def sgp4_update(satrec: Satrec, tsince): xlm = xlm % tau mm = (xlm - argpm - nodem) % tau - # sgp4fix recover singly averaged mean elements satrec.am = am; satrec.em = em; satrec.im = inclm; @@ -129,8 +102,8 @@ def sgp4_update(satrec: Satrec, tsince): satrec.nm = nm; # ----------------- compute extra mean quantities ------------- - sinim = np.sin(inclm); - cosim = np.cos(inclm); + sinim = np.sin(inclm) + cosim = np.cos(inclm) # -------------------- add lunar-solar periodics -------------- ep = em; @@ -152,8 +125,7 @@ def sgp4_update(satrec: Satrec, tsince): eo1 = u; tem5 = 9999.9; ktr = 1; - # sgp4fix for kepler iteration - # the following iteration needs better limits on corrections + while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: sineo1 = np.sin(eo1); @@ -171,28 +143,25 @@ def sgp4_update(satrec: Satrec, tsince): el2 = axnl*axnl + aynl*aynl; pl = am*(1.0-el2); if pl < 0.0: - - satrec.error_message = ('semilatus rectum {0:f} is less than zero' - .format(pl)) - satrec.error = SEMIRECT - # sgp4fix add return + satrec.error = satrec.SEMIRECT + return False, False; else: - rl = am * (1.0 - ecose); - rdotl = np.sqrt(am) * esine/rl; - rvdotl = np.sqrt(pl) / rl; - betal = np.sqrt(1.0 - el2); - temp = esine / (1.0 + betal); - sinu = am / rl * (sineo1 - aynl - axnl * temp); - cosu = am / rl * (coseo1 - axnl + aynl * temp); - su = np.atan2(sinu, cosu); - sin2u = (cosu + cosu) * sinu; - cos2u = 1.0 - 2.0 * sinu * sinu; - temp = 1.0 / pl; - temp1 = 0.5 * satrec.j2 * temp; - temp2 = temp1 * temp; + rl = am * (1.0 - ecose) + rdotl = np.sqrt(am) * esine/rl + rvdotl = np.sqrt(pl) / rl + betal = np.sqrt(1.0 - el2) + temp = esine / (1.0 + betal) + sinu = am / rl * (sineo1 - aynl - axnl * temp) + cosu = am / rl * (coseo1 - axnl + aynl * temp) + su = np.atan2(sinu, cosu) + sin2u = (cosu + cosu) * sinu + cos2u = 1.0 - 2.0 * sinu * sinu + temp = 1.0 / pl + temp1 = 0.5 * satrec.j2 * temp + temp2 = temp1 * temp # -------------- update for short period periodics ------------ mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ @@ -227,12 +196,8 @@ def sgp4_update(satrec: Satrec, tsince): (mvt * uy + rvdot * vy) * vkmpersec, (mvt * uz + rvdot * vz) * vkmpersec) - # sgp4fix for decaying satellites if mrt < 1.0: - - satrec.error_message = ('mrt {0:f} is less than 1.0 indicating' - ' the satellite has decayed'.format(mrt)) - satrec.error = DECAY + satrec.error = satrec.DECAY return r, v From 232262161f4509000857484c8672c2f2f2c070e8 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 27 Apr 2026 18:08:09 -0700 Subject: [PATCH 04/19] continuing to implement satrec and associated logical flow --- src/lib/sgp4.py | 410 +++++++++++++++++++++++++++++++++++++++++ src/lib/tle.py | 232 ++--------------------- src/tasks/adcs/sgp4.py | 204 -------------------- 3 files changed, 425 insertions(+), 421 deletions(-) create mode 100644 src/lib/sgp4.py delete mode 100644 src/tasks/adcs/sgp4.py diff --git a/src/lib/sgp4.py b/src/lib/sgp4.py new file mode 100644 index 0000000..a0f88a0 --- /dev/null +++ b/src/lib/sgp4.py @@ -0,0 +1,410 @@ +""" +Simplified General Perturbations Model 4 Implementation for Orbit Propagation + +As outlined in https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf +""" + +# Error codes +ECCENTRICITY = 1 # eccentricity is not within 0-1 +MOTION = 2 # error in propagating mean motion +SEMIRECT = 4 # apoapsis, periapsis characteristics error +DECAY = 6 # orbit has decayed + +try: + import ulab.numpy as np # For CircuitPython +except ImportError: + import numpy as np # For GitHub Actions / PC testing + +def _gstime(jdut1): + deg2rad = np.pi / 180.0 + tau = np.pi*2 + tut1 = (jdut1 - 2451545.0) / 36525.0 + temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \ + (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841 # sec + temp = (temp * deg2rad / 240.0) % tau # 360/86400 = 1/240, to deg, to rad + + # ------------------------ check quadrants --------------------- + if temp < 0.0: + temp += tau + + return temp + +def _initl(xke, j2, + ecco, epoch, inclo, no): + # ----------------------- earth constants ---------------------- + # sgp4fix identify constants and allow alternate values + # only xke and j2 are used here so pass them in directly + # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst + x2o3 = 2.0 / 3.0; + + # ------------- calculate auxillary epoch quantities ---------- + eccsq = ecco * ecco; + omeosq = 1.0 - eccsq; + rteosq = np.sqrt(omeosq); + cosio = np.cos(inclo); + cosio2 = cosio * cosio; + + # ------------------ un-kozai the mean motion ----------------- + ak = pow(xke / no, x2o3); + d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); + del_ = d1 / (ak * ak); + adel = ak * (1.0 - del_ * del_ - del_ * + (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); + del_ = d1/(adel * adel) + no = no / (1.0 + del_) + + ao = pow(xke / no, x2o3) + sinio = np.sin(inclo) + po = ao * omeosq + con42 = 1.0 - 5.0 * cosio2 + con41 = -con42-cosio2-cosio2 + posq = po * po + rp = ao * (1.0 - ecco) + + gsto = _gstime(epoch + 2433281.5) + + return (no, + ao, con41, con42, cosio, + cosio2, omeosq, posq, + rp, rteosq,sinio , gsto + ) + +def sgp4_update(satrec, tsince): + """ + Transforms to cartesian (r, v) from a Satrec (TLE 7+) input, then performs matrix operations + to propagate the system + + Current implementation uses an analytical approach to compute the Jacobian A matrix + Possible alternative is numerical stepping (for lower accuracy) of associated partials + """ + # -- Set mathematical constants + x2o3 = 2.0 / 3.0 + tau = 2.0 * np.pi + vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; + + satrec.t = tsince + + # -- Update for secular gravity and atmospheric drag + xmdf = satrec.mo + satrec.mdot * satrec.t; + argpdf = satrec.argpo + satrec.argpdot * satrec.t; + nodedf = satrec.nodeo + satrec.nodedot * satrec.t; + argpm = argpdf + mm = xmdf + t2 = satrec.t * satrec.t; + nodem = nodedf + satrec.nodecf * t2; + tempa = 1.0 - satrec.cc1 * satrec.t; + tempe = satrec.bstar * satrec.cc4 * satrec.t; + templ = satrec.t2cof * t2; + + if satrec.isimp != 1: + + delomg = satrec.omgcof * satrec.t; + # sgp4fix use mutliply for speed instead of pow + delmtemp = 1.0 + satrec.eta * np.cos(xmdf); + delm = satrec.xmcof * \ + (delmtemp * delmtemp * delmtemp - + satrec.delmo); + temp = delomg + delm; + mm = xmdf + temp; + argpm = argpdf - temp; + t3 = t2 * satrec.t; + t4 = t3 * satrec.t; + tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ + satrec.d4 * t4; + tempe = tempe + satrec.bstar * satrec.cc5 * (np.sin(mm) - + satrec.sinmao); + templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + + satrec.t * satrec.t5cof); + + nm = satrec.no_unkozai; + em = satrec.ecco; + inclm = satrec.inclo; + + if nm <= 0.0: + satrec.error = satrec.MOTION + return False, False; + + am = pow((satrec.xke / nm),x2o3) * tempa * tempa; + nm = satrec.xke / pow(am, 1.5); + em = em - tempe; + + if em >= 1.0 or em < -0.001: + satrec.error = satrec.ECCENTRICITY + + return False, False; + + if em < 1.0e-6: + em = 1.0e-6; + mm = mm + satrec.no_unkozai * templ; + xlm = mm + argpm + nodem; + emsq = em * em; + temp = 1.0 - emsq; + + nodem = nodem % tau if nodem >= 0.0 else -(-nodem % tau) + argpm = argpm % tau + xlm = xlm % tau + mm = (xlm - argpm - nodem) % tau + + satrec.am = am + satrec.em = em + satrec.im = inclm + satrec.Om = nodem + satrec.om = argpm + satrec.mm = mm + satrec.nm = nm + + # ----------------- compute extra mean quantities ------------- + sinim = np.sin(inclm) + cosim = np.cos(inclm) + + # -------------------- add lunar-solar periodics -------------- + ep = em; + xincp = inclm; + argpp = argpm; + nodep = nodem; + mp = mm; + sinip = sinim; + cosip = cosim; + + # -------------------- long period periodics ------------------ + axnl = ep * np.cos(argpp); + temp = 1.0 / (am * (1.0 - ep * ep)); + aynl = ep* np.sin(argpp) + temp * satrec.aycof; + xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; + + # --------------------- solve kepler's equation --------------- + u = (xl - nodep) % tau + eo1 = u; + tem5 = 9999.9; + ktr = 1; + + while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: + + sineo1 = np.sin(eo1); + coseo1 = np.cos(eo1); + tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; + tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; + if np.fabs(tem5) >= 0.95: + tem5 = 0.95 if tem5 > 0.0 else -0.95; + eo1 = eo1 + tem5; + ktr = ktr + 1; + + # ------------- short period preliminary quantities ----------- + ecose = axnl*coseo1 + aynl*sineo1; + esine = axnl*sineo1 - aynl*coseo1; + el2 = axnl*axnl + aynl*aynl; + pl = am*(1.0-el2); + if pl < 0.0: + satrec.error = satrec.SEMIRECT + + return False, False; + + else: + + rl = am * (1.0 - ecose) + rdotl = np.sqrt(am) * esine/rl + rvdotl = np.sqrt(pl) / rl + betal = np.sqrt(1.0 - el2) + temp = esine / (1.0 + betal) + sinu = am / rl * (sineo1 - aynl - axnl * temp) + cosu = am / rl * (coseo1 - axnl + aynl * temp) + su = np.atan2(sinu, cosu) + sin2u = (cosu + cosu) * sinu + cos2u = 1.0 - 2.0 * sinu * sinu + temp = 1.0 / pl + temp1 = 0.5 * satrec.j2 * temp + temp2 = temp1 * temp + + # -------------- update for short period periodics ------------ + mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ + 0.5 * temp1 * satrec.x1mth2 * cos2u; + su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; + xnode = nodep + 1.5 * temp2 * cosip * sin2u; + xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; + mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke; + rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + + 1.5 * satrec.con41) / satrec.xke; + + # --------------------- orientation vectors ------------------- + sinsu = np.sin(su); + cossu = np.cos(su); + snod = np.sin(xnode); + cnod = np.cos(xnode); + sini = np.sin(xinc); + cosi = np.cos(xinc); + xmx = -snod * cosi; + xmy = cnod * cosi; + ux = xmx * sinsu + cnod * cossu; + uy = xmy * sinsu + snod * cossu; + uz = sini * sinsu; + vx = xmx * cossu - cnod * sinsu; + vy = xmy * cossu - snod * sinsu; + vz = sini * cossu; + + # --------- position and velocity (in km and km/sec) ---------- + _mr = mrt * satrec.radiusearthkm + r = (_mr * ux, _mr * uy, _mr * uz) + v = ((mvt * ux + rvdot * vx) * vkmpersec, + (mvt * uy + rvdot * vy) * vkmpersec, + (mvt * uz + rvdot * vz) * vkmpersec) + + if mrt < 1.0: + satrec.error = satrec.DECAY + + return r, v + +def sgp4_init(satrec, satn, epoch, + bstar, ndot, nddot, ecco, argpo, + inclo, mo, no_kozai, + nodeo + ): + + temp4 = 1.5e-12 + + # Near Earth Variables + satrec.isimp = 0; satrec.aycof = 0.0 + satrec.con41 = 0.0; satrec.cc1 = 0.0; satrec.cc4 = 0.0 + satrec.cc5 = 0.0; satrec.d2 = 0.0; satrec.d3 = 0.0 + satrec.d4 = 0.0; satrec.delmo = 0.0; satrec.eta = 0.0 + satrec.argpdot = 0.0; satrec.omgcof = 0.0; satrec.sinmao = 0.0 + satrec.t = 0.0; satrec.t2cof = 0.0; satrec.t3cof = 0.0 + satrec.t4cof = 0.0; satrec.t5cof = 0.0; satrec.x1mth2 = 0.0 + satrec.x7thm1 = 0.0; satrec.mdot = 0.0; satrec.nodedot = 0.0 + satrec.xlcof = 0.0; satrec.xmcof = 0.0; satrec.nodecf = 0.0 + + # Earth Constants + satrec.mu = 398600.5; # in km3 / s2 + satrec.radiusearthkm = 6378.137 # km + satrec.xke = 0.07436685317 + satrec.tumin = 13.44685108 + satrec.j2 = 0.00108262998905 + satrec.j3 = -0.00000253215306 + satrec.j4 = -0.00000161098761 + satrec.j3oj2 = -0.002338890559 + + ss = 1.012229276 + qzms2ttemp = 0.00658499496 + qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; + x2o3 = 2.0 / 3.0 + + # -- initialisation markers + satrec.t = 0.0 + + # -- + satrec.satnum_str = satn + satrec.classification = 'U' + + # -- + satrec.bstar = bstar + satrec.ndot = ndot + satrec.nddot = nddot + satrec.ecco = ecco + satrec.argpo = argpo + satrec.inclo = inclo + satrec.mo = mo + satrec.nodeo = nodeo + satrec.no_kozai = no_kozai + + # single averaged mean elements + satrec.am = 0.0 + satrec.em = 0.0 + satrec.im = 0.0 + satrec.Om = 0.0 + satrec.mm = 0.0 + satrec.nm = 0.0 + + satrec.error = 0 + + # -- + ( + satrec.no_unkozai, + ao, satrec.con41, con42, cosio, + cosio2, omeosq, posq, + rp, rteosq,sinio , satrec.gsto, + ) = _initl( + satrec.xke, satrec.j2, satrec.ecco, epoch, satrec.inclo, satrec.no_kozai + ) + satrec.a = pow( satrec.no_unkozai*satrec.tumin , (-2.0/3.0) ); + satrec.alta = satrec.a*(1.0 + satrec.ecco) - 1.0; + satrec.altp = satrec.a*(1.0 - satrec.ecco) - 1.0; + + if omeosq >= 0.0 or satrec.no_unkozai >= 0.0: + satrec.isimp = 0 + if rp < 220.0 / satrec.radiusearthkm + 1.0: + satrec.isimp = 1 + sfour = ss + qzms24 = qzms2t + pinvsq = 1.0 / posq; + + tsi = 1.0 / (ao - sfour); + satrec.eta = ao * satrec.ecco * tsi; + etasq = satrec.eta * satrec.eta; + eeta = satrec.ecco * satrec.eta; + psisq = np.fabs(1.0 - etasq); + coef = qzms24 * pow(tsi, 4.0); + coef1 = coef / pow(psisq, 3.5); + cc2 = coef1 * satrec.no_unkozai * (ao * (1.0 + 1.5 * etasq + eeta * + (4.0 + etasq)) + 0.375 * satrec.j2 * tsi / psisq * satrec.con41 * + (8.0 + 3.0 * etasq * (8.0 + etasq))); + satrec.cc1 = satrec.bstar * cc2; + cc3 = 0.0; + if satrec.ecco > 1.0e-4: + cc3 = -2.0 * coef * tsi * satrec.j3oj2 * satrec.no_unkozai * sinio / satrec.ecco; + satrec.x1mth2 = 1.0 - cosio2; + satrec.cc4 = 2.0* satrec.no_unkozai * coef1 * ao * omeosq * \ + (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco * + (0.5 + 2.0 * etasq) - satrec.j2 * tsi / (ao * psisq) * + (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * + (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * + (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * satrec.argpo))); + satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * + (etasq + eeta) + eeta * etasq); + cosio4 = cosio2 * cosio2; + temp1 = 1.5 * satrec.j2 * pinvsq * satrec.no_unkozai; + temp2 = 0.5 * temp1 * satrec.j2 * pinvsq; + temp3 = -0.46875 * satrec.j4 * pinvsq * pinvsq * satrec.no_unkozai; + satrec.mdot = satrec.no_unkozai + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * \ + temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); + satrec.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * + (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + + temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)); + xhdot1 = -temp1 * cosio; + satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; + satrec.omgcof = satrec.bstar * cc3 * np.cos(satrec.argpo); + satrec.xmcof = 0.0; + if satrec.ecco > 1.0e-4: + satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta; + satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1; + satrec.t2cof = 1.5 * satrec.cc1; + + if np.fabs(cosio+1.0) > 1.5e-12: + satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); + else: + satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; + satrec.aycof = -0.5 * satrec.j3oj2 * sinio; + + delmotemp = 1.0 + satrec.eta * np.cos(satrec.mo); + satrec.delmo = delmotemp * delmotemp * delmotemp; + satrec.sinmao = np.sin(satrec.mo); + satrec.x7thm1 = 7.0 * cosio2 - 1.0; + + if satrec.isimp != 1: + cc1sq = satrec.cc1 * satrec.cc1; + satrec.d2 = 4.0 * ao * tsi * cc1sq; + temp = satrec.d2 * tsi * satrec.cc1 / 3.0; + satrec.d3 = (17.0 * ao + sfour) * temp; + satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ + satrec.cc1; + satrec.t3cof = satrec.d2 + 2.0 * cc1sq; + satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * + (12.0 * satrec.d2 + 10.0 * cc1sq)); + satrec.t5cof = 0.2 * (3.0 * satrec.d4 + + 12.0 * satrec.cc1 * satrec.d3 + + 6.0 * satrec.d2 * satrec.d2 + + 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)); + + # propagate to 0 + sgp4_update(satrec, 0) + + return True diff --git a/src/lib/tle.py b/src/lib/tle.py index cbcad74..ac8c0ea 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -4,71 +4,14 @@ Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, """ -try: - import ulab.numpy as np # For CircuitPython -except ImportError: - import numpy as np # For GitHub Actions / PC testing - -def _gstime(jdut1): - deg2rad = np.pi / 180.0 - tau = np.pi*2 - tut1 = (jdut1 - 2451545.0) / 36525.0 - temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \ - (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841 # sec - temp = (temp * deg2rad / 240.0) % tau # 360/86400 = 1/240, to deg, to rad - - # ------------------------ check quadrants --------------------- - if temp < 0.0: - temp += tau - - return temp - -def _initl(xke, j2, - ecco, epoch, inclo, no): - # ----------------------- earth constants ---------------------- - # sgp4fix identify constants and allow alternate values - # only xke and j2 are used here so pass them in directly - # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst - x2o3 = 2.0 / 3.0; - - # ------------- calculate auxillary epoch quantities ---------- - eccsq = ecco * ecco; - omeosq = 1.0 - eccsq; - rteosq = np.sqrt(omeosq); - cosio = np.cos(inclo); - cosio2 = cosio * cosio; - - # ------------------ un-kozai the mean motion ----------------- - ak = pow(xke / no, x2o3); - d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); - del_ = d1 / (ak * ak); - adel = ak * (1.0 - del_ * del_ - del_ * - (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); - del_ = d1/(adel * adel) - no = no / (1.0 + del_) - - ao = pow(xke / no, x2o3) - sinio = np.sin(inclo) - po = ao * omeosq - con42 = 1.0 - 5.0 * cosio2 - con41 = -con42-cosio2-cosio2 - posq = po * po - rp = ao * (1.0 - ecco) - - gsto = _gstime(epoch + 2433281.5) - - return ( - no, - ao, con41, con42, cosio, - cosio2, omeosq, posq, - rp, rteosq,sinio , gsto - ) +from sgp4 import sgp4_init, sgp4_update class Satrec: """ + Satellite record object Parameters, constants that are commonly used across the sgp4 logical flow - Usually built from TLE + In this implementation, built from TLE data """ # Error codes @@ -89,7 +32,7 @@ class Satrec: error, error_message - might want to leave this up to cdh but implement this here? isimp - using simplified model? - + revnum - """ """ Error Messages @@ -110,162 +53,9 @@ class Satrec: ' the satellite has decayed'.format(mrt)) """ - def __init__(self, satn, epoch, - bstar, ndot, nddot, ecco, argpo, - inclo, mo, no_kozai, - nodeo - ): - - temp4 = 1.5e-12 - - # Near Earth Variables - self.isimp = 0; self.aycof = 0.0 - self.con41 = 0.0; self.cc1 = 0.0; self.cc4 = 0.0 - self.cc5 = 0.0; self.d2 = 0.0; self.d3 = 0.0 - self.d4 = 0.0; self.delmo = 0.0; self.eta = 0.0 - self.argpdot = 0.0; self.omgcof = 0.0; self.sinmao = 0.0 - self.t = 0.0; self.t2cof = 0.0; self.t3cof = 0.0 - self.t4cof = 0.0; self.t5cof = 0.0; self.x1mth2 = 0.0 - self.x7thm1 = 0.0; self.mdot = 0.0; self.nodedot = 0.0 - self.xlcof = 0.0; self.xmcof = 0.0; self.nodecf = 0.0 - - # Earth Constants - self.mu = 398600.5; # in km3 / s2 - self.radiusearthkm = 6378.137 # km - self.xke = 0.07436685317 - self.tumin = 13.44685108 - self.j2 = 0.00108262998905 - self.j3 = -0.00000253215306 - self.j4 = -0.00000161098761 - self.j3oj2 = -0.002338890559 - - ss = 1.012229276 - qzms2ttemp = 0.00658499496 - qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; - x2o3 = 2.0 / 3.0 - - # -- initialisation markers - self.init = 'y' - self.t = 0.0 - - # -- - self.satnum_str = satn - self.classification = 'U' - - # -- - self.bstar = bstar - self.ndot = ndot - self.nddot = nddot - self.ecco = ecco - self.argpo = argpo - self.inclo = inclo - self.mo = mo - self.nodeo = nodeo - self.no_kozai = no_kozai - - # single averaged mean elements - self.am = 0.0 - self.em = 0.0 - self.im = 0.0 - self.Om = 0.0 - self.mm = 0.0 - self.nm = 0.0 - - self.error = 0 - - # -- - ( - self.no_unkozai, - ao, self.con41, con42, cosio, - cosio2, omeosq, posq, - rp, rteosq,sinio , self.gsto, - ) = _initl( - self.xke, self.j2, self.ecco, epoch, self.inclo, self.no_kozai - ) - self.a = pow( self.no_unkozai*self.tumin , (-2.0/3.0) ); - self.alta = self.a*(1.0 + self.ecco) - 1.0; - self.altp = self.a*(1.0 - self.ecco) - 1.0; - - if omeosq >= 0.0 or self.no_unkozai >= 0.0: - self.isimp = 0 - if rp < 220.0 / self.radiusearthkm + 1.0: - self.isimp = 1 - sfour = ss - qzms24 = qzms2t - pinvsq = 1.0 / posq; - - tsi = 1.0 / (ao - sfour); - self.eta = ao * self.ecco * tsi; - etasq = self.eta * self.eta; - eeta = self.ecco * self.eta; - psisq = np.fabs(1.0 - etasq); - coef = qzms24 * pow(tsi, 4.0); - coef1 = coef / pow(psisq, 3.5); - cc2 = coef1 * self.no_unkozai * (ao * (1.0 + 1.5 * etasq + eeta * - (4.0 + etasq)) + 0.375 * self.j2 * tsi / psisq * self.con41 * - (8.0 + 3.0 * etasq * (8.0 + etasq))); - self.cc1 = self.bstar * cc2; - cc3 = 0.0; - if self.ecco > 1.0e-4: - cc3 = -2.0 * coef * tsi * self.j3oj2 * self.no_unkozai * sinio / self.ecco; - self.x1mth2 = 1.0 - cosio2; - self.cc4 = 2.0* self.no_unkozai * coef1 * ao * omeosq * \ - (self.eta * (2.0 + 0.5 * etasq) + self.ecco * - (0.5 + 2.0 * etasq) - self.j2 * tsi / (ao * psisq) * - (-3.0 * self.con41 * (1.0 - 2.0 * eeta + etasq * - (1.5 - 0.5 * eeta)) + 0.75 * self.x1mth2 * - (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * self.argpo))); - self.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * - (etasq + eeta) + eeta * etasq); - cosio4 = cosio2 * cosio2; - temp1 = 1.5 * self.j2 * pinvsq * self.no_unkozai; - temp2 = 0.5 * temp1 * self.j2 * pinvsq; - temp3 = -0.46875 * self.j4 * pinvsq * pinvsq * self.no_unkozai; - self.mdot = self.no_unkozai + 0.5 * temp1 * rteosq * self.con41 + 0.0625 * \ - temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); - self.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * - (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + - temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)); - xhdot1 = -temp1 * cosio; - self.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + - 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; - self.omgcof = self.bstar * cc3 * np.cos(self.argpo); - self.xmcof = 0.0; - if self.ecco > 1.0e-4: - self.xmcof = -x2o3 * coef * self.bstar / eeta; - self.nodecf = 3.5 * omeosq * xhdot1 * self.cc1; - self.t2cof = 1.5 * self.cc1; - - if np.fabs(cosio+1.0) > 1.5e-12: - self.xlcof = -0.25 * self.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); - else: - self.xlcof = -0.25 * self.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; - self.aycof = -0.5 * self.j3oj2 * sinio; - - delmotemp = 1.0 + self.eta * np.cos(self.mo); - self.delmo = delmotemp * delmotemp * delmotemp; - self.sinmao = np.sin(self.mo); - self.x7thm1 = 7.0 * cosio2 - 1.0; - - if self.isimp != 1: - cc1sq = self.cc1 * self.cc1; - self.d2 = 4.0 * ao * tsi * cc1sq; - temp = self.d2 * tsi * self.cc1 / 3.0; - self.d3 = (17.0 * ao + sfour) * temp; - self.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ - self.cc1; - self.t3cof = self.d2 + 2.0 * cc1sq; - self.t4cof = 0.25 * (3.0 * self.d3 + self.cc1 * - (12.0 * self.d2 + 10.0 * cc1sq)); - self.t5cof = 0.2 * (3.0 * self.d4 + - 12.0 * self.cc1 * self.d3 + - 6.0 * self.d2 * self.d2 + - 15.0 * cc1sq * (2.0 * self.d2 + cc1sq)); - - # need to propagate to epoch 0.0 to really instantiate everything else before - # init flag is set to 'n' - - return True + def __init__(self): + self.revnum = 0 + pass @classmethod def from_tle_array(cls): @@ -273,6 +63,14 @@ def from_tle_array(cls): obj = cls() return obj + + def sgp4_update(self, jd, fr): + """ + For a julian date (jd) and its fractional representation (fr), + propagate Satrec using sgp4 + """ + + sgp4_update() def _conv_year(s): """Interpret a two-digit year string.""" diff --git a/src/tasks/adcs/sgp4.py b/src/tasks/adcs/sgp4.py deleted file mode 100644 index 606f19e..0000000 --- a/src/tasks/adcs/sgp4.py +++ /dev/null @@ -1,204 +0,0 @@ -""" -Simplified General Perturbations Model 4 Implementation for Orbit Propagation - -As outlined in https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf -""" - -from datastores.adcs import Satrec - -# Error codes -ECCENTRICITY = 1 # eccentricity is not within 0-1 -MOTION = 2 # error in propagating mean motion -SEMIRECT = 4 # apoapsis, periapsis characteristics error -DECAY = 6 # orbit has decayed - -try: - import ulab.numpy as np # For CircuitPython -except ImportError: - import numpy as np # For GitHub Actions / PC testing - -def sgp4_update(satrec: Satrec, tsince): - """ - Transforms to cartesian (r, v) from a Satrec (TLE 7+) input, then performs matrix operations - to propagate the system - - Current implementation uses an analytical approach to compute the Jacobian A matrix - Possible alternative is numerical stepping (for lower accuracy) of associated partials - """ - # -- Set mathematical constants - x2o3 = 2.0 / 3.0 - tau = 2.0 * np.pi - vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; - - satrec.t = tsince - - # -- Update for secular gravity and atmospheric drag - xmdf = satrec.mo + satrec.mdot * satrec.t; - argpdf = satrec.argpo + satrec.argpdot * satrec.t; - nodedf = satrec.nodeo + satrec.nodedot * satrec.t; - argpm = argpdf - mm = xmdf - t2 = satrec.t * satrec.t; - nodem = nodedf + satrec.nodecf * t2; - tempa = 1.0 - satrec.cc1 * satrec.t; - tempe = satrec.bstar * satrec.cc4 * satrec.t; - templ = satrec.t2cof * t2; - - if satrec.isimp != 1: - - delomg = satrec.omgcof * satrec.t; - # sgp4fix use mutliply for speed instead of pow - delmtemp = 1.0 + satrec.eta * np.cos(xmdf); - delm = satrec.xmcof * \ - (delmtemp * delmtemp * delmtemp - - satrec.delmo); - temp = delomg + delm; - mm = xmdf + temp; - argpm = argpdf - temp; - t3 = t2 * satrec.t; - t4 = t3 * satrec.t; - tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ - satrec.d4 * t4; - tempe = tempe + satrec.bstar * satrec.cc5 * (np.sin(mm) - - satrec.sinmao); - templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + - satrec.t * satrec.t5cof); - - nm = satrec.no_unkozai; - em = satrec.ecco; - inclm = satrec.inclo; - - if nm <= 0.0: - satrec.error = satrec.MOTION - return False, False; - - am = pow((satrec.xke / nm),x2o3) * tempa * tempa; - nm = satrec.xke / pow(am, 1.5); - em = em - tempe; - - if em >= 1.0 or em < -0.001: - satrec.error = satrec.ECCENTRICITY - - return False, False; - - if em < 1.0e-6: - em = 1.0e-6; - mm = mm + satrec.no_unkozai * templ; - xlm = mm + argpm + nodem; - emsq = em * em; - temp = 1.0 - emsq; - - nodem = nodem % tau if nodem >= 0.0 else -(-nodem % tau) - argpm = argpm % tau - xlm = xlm % tau - mm = (xlm - argpm - nodem) % tau - - satrec.am = am; - satrec.em = em; - satrec.im = inclm; - satrec.Om = nodem; - satrec.om = argpm; - satrec.mm = mm; - satrec.nm = nm; - - # ----------------- compute extra mean quantities ------------- - sinim = np.sin(inclm) - cosim = np.cos(inclm) - - # -------------------- add lunar-solar periodics -------------- - ep = em; - xincp = inclm; - argpp = argpm; - nodep = nodem; - mp = mm; - sinip = sinim; - cosip = cosim; - - # -------------------- long period periodics ------------------ - axnl = ep * np.cos(argpp); - temp = 1.0 / (am * (1.0 - ep * ep)); - aynl = ep* np.sin(argpp) + temp * satrec.aycof; - xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; - - # --------------------- solve kepler's equation --------------- - u = (xl - nodep) % tau - eo1 = u; - tem5 = 9999.9; - ktr = 1; - - while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: - - sineo1 = np.sin(eo1); - coseo1 = np.cos(eo1); - tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; - tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; - if np.fabs(tem5) >= 0.95: - tem5 = 0.95 if tem5 > 0.0 else -0.95; - eo1 = eo1 + tem5; - ktr = ktr + 1; - - # ------------- short period preliminary quantities ----------- - ecose = axnl*coseo1 + aynl*sineo1; - esine = axnl*sineo1 - aynl*coseo1; - el2 = axnl*axnl + aynl*aynl; - pl = am*(1.0-el2); - if pl < 0.0: - satrec.error = satrec.SEMIRECT - - return False, False; - - else: - - rl = am * (1.0 - ecose) - rdotl = np.sqrt(am) * esine/rl - rvdotl = np.sqrt(pl) / rl - betal = np.sqrt(1.0 - el2) - temp = esine / (1.0 + betal) - sinu = am / rl * (sineo1 - aynl - axnl * temp) - cosu = am / rl * (coseo1 - axnl + aynl * temp) - su = np.atan2(sinu, cosu) - sin2u = (cosu + cosu) * sinu - cos2u = 1.0 - 2.0 * sinu * sinu - temp = 1.0 / pl - temp1 = 0.5 * satrec.j2 * temp - temp2 = temp1 * temp - - # -------------- update for short period periodics ------------ - mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ - 0.5 * temp1 * satrec.x1mth2 * cos2u; - su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; - xnode = nodep + 1.5 * temp2 * cosip * sin2u; - xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; - mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke; - rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + - 1.5 * satrec.con41) / satrec.xke; - - # --------------------- orientation vectors ------------------- - sinsu = np.sin(su); - cossu = np.cos(su); - snod = np.sin(xnode); - cnod = np.cos(xnode); - sini = np.sin(xinc); - cosi = np.cos(xinc); - xmx = -snod * cosi; - xmy = cnod * cosi; - ux = xmx * sinsu + cnod * cossu; - uy = xmy * sinsu + snod * cossu; - uz = sini * sinsu; - vx = xmx * cossu - cnod * sinsu; - vy = xmy * cossu - snod * sinsu; - vz = sini * cossu; - - # --------- position and velocity (in km and km/sec) ---------- - _mr = mrt * satrec.radiusearthkm - r = (_mr * ux, _mr * uy, _mr * uz) - v = ((mvt * ux + rvdot * vx) * vkmpersec, - (mvt * uy + rvdot * vy) * vkmpersec, - (mvt * uz + rvdot * vz) * vkmpersec) - - if mrt < 1.0: - satrec.error = satrec.DECAY - - return r, v - - From 40a4bf0d7784b6d314c912566433dd036a761524 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Fri, 1 May 2026 17:46:41 -0700 Subject: [PATCH 05/19] TLE and Satrec consolidated, ready to begin unittest --- src/lib/datastores/adcs.py | 2 +- src/lib/sgp4.py | 5 +- src/lib/tle.py | 251 +++++++++++++++++++++++------------- unit_tests/lib/sgp4_test.py | 11 ++ 4 files changed, 177 insertions(+), 92 deletions(-) create mode 100644 unit_tests/lib/sgp4_test.py diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 7ae8d3c..13d70fb 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -4,7 +4,7 @@ are set to `None` throughout this module. """ -import tle +import tle as tle class Datastore: """ diff --git a/src/lib/sgp4.py b/src/lib/sgp4.py index a0f88a0..305514e 100644 --- a/src/lib/sgp4.py +++ b/src/lib/sgp4.py @@ -255,7 +255,7 @@ def sgp4_update(satrec, tsince): def sgp4_init(satrec, satn, epoch, bstar, ndot, nddot, ecco, argpo, - inclo, mo, no_kozai, + inclo, mo, nodeo ): @@ -303,7 +303,6 @@ def sgp4_init(satrec, satn, epoch, satrec.inclo = inclo satrec.mo = mo satrec.nodeo = nodeo - satrec.no_kozai = no_kozai # single averaged mean elements satrec.am = 0.0 @@ -322,7 +321,7 @@ def sgp4_init(satrec, satn, epoch, cosio2, omeosq, posq, rp, rteosq,sinio , satrec.gsto, ) = _initl( - satrec.xke, satrec.j2, satrec.ecco, epoch, satrec.inclo, satrec.no_kozai + satrec.xke, satrec.j2, satrec.ecco, epoch, satrec.inclo ) satrec.a = pow( satrec.no_unkozai*satrec.tumin , (-2.0/3.0) ); satrec.alta = satrec.a*(1.0 + satrec.ecco) - 1.0; diff --git a/src/lib/tle.py b/src/lib/tle.py index ac8c0ea..9b08a88 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -6,71 +6,78 @@ from sgp4 import sgp4_init, sgp4_update -class Satrec: +try: + import ulab.numpy as np # For CircuitPython +except ImportError: + import numpy as np # For GitHub Actions / PC testing + +minutes_per_day = 1440 +epoch0 = 2433281.5 # jan 0 1950 + +def _day_of_year_to_month_day(day_of_year, is_leap): + """Core logic for turning days into months, for easy testing.""" + february_bump = (2 - is_leap) * (day_of_year >= 60 + is_leap) + august = day_of_year >= 215 + month, day = divmod(2 * (day_of_year - 1 + 30 * august + february_bump), 61) + month += 1 - august + day //= 2 + day += 1 + return month, day + +def _jday(year, mon, day, hr, minute, sec): + """Return two floats (jd, fr) that, when added, produce the specified Julian date + jd (Julian Date) and fr (Fractional) components + + >>> jd, fr = jday(2020, 2, 11, 13, 57, 0) + >>> jd + 2458890.5 + >>> fr + 0.58125 + + Note the first float, which gives the moment of midnight that + commences the given calendar date, always ends in + ``.5`` because Julian dates begin and end at noon. This made + Julian dates more convenient for astronomers in Europe, by making + the whole night belong to a single Julian date. """ - Satellite record object - Parameters, constants that are commonly used across the sgp4 logical flow - - In this implementation, built from TLE data + jd = (367.0 * year + - 7 * (year + ((mon + 9) // 12.0)) * 0.25 // 1.0 + + 275 * mon / 9.0 // 1.0 + + day + + 1721013.5) + fr = (sec + minute * 60.0 + hr * 3600.0) / 86400.0; + return jd, fr + +def _days2mdhms(year, days, round_to_microsecond=6): + """Convert a float point number of days into the year into date and time. + + >>> days2mdhms(2000, 32.0) # February 1 + (2, 1, 0, 0, 0.0) + >>> days2mdhms(2000, 366.0) # December 31, since 2000 was a leap year + (12, 31, 0, 0, 0.0) + + The floating point seconds are rounded to an even number of + microseconds if ``round_to_microsecond`` is true. """ + second = days * 86400.0 + if round_to_microsecond: + second = round(second, round_to_microsecond) - # Error codes - ECCENTRICITY = 1 # eccentricity is not within 0-1 - MOTION = 2 # error in propagating mean motion - SEMIRECT = 4 # apoapsis, periapsis characteristics error - DECAY = 6 # orbit has decayed - - """ - t - time since - mo, mdot - mean anomaly - argpo, argpdot - argument of perigee - nodeo, nodedot, nodecf - RAAN value, drift, and correction respectively - bstar - atmospheric drag - cc1, cc4, cc5 - drag coefficient terms - ecco, inclo - eccentricity, inclination (no derivatives) - mm, nm - mean motion before and after corrections - error, error_message - might want to leave this up to cdh but implement - this here? - isimp - using simplified model? - revnum - - """ + minute, second = divmod(second, 60.0) + if round_to_microsecond: + second = round(second, round_to_microsecond) - """ Error Messages - MOTION - satrec.error_message = ('mean motion {0:f} is less than zero' - .format(nm)) - - ECCENTRICITY - satrec.error_message = ('mean eccentricity {0:f} not within' - ' range 0.0 <= e < 1.0'.format(em)) + minute = int(minute) + hour, minute = divmod(minute, 60) + day_of_year, hour = divmod(hour, 24) - SEMIRECT - satrec.error_message = ('semilatus rectum {0:f} is less than zero' - .format(pl)) - - DECAY - satrec.error_message = ('mrt {0:f} is less than 1.0 indicating' - ' the satellite has decayed'.format(mrt)) - """ + is_leap = year % 400 == 0 or (year % 4 == 0 and year % 100 != 0) + month, day = _day_of_year_to_month_day(day_of_year, is_leap) + if month == 13: # behave like the original in case of overflow + month = 12 + day += 31 - def __init__(self): - self.revnum = 0 - pass - - @classmethod - def from_tle_array(cls): - - obj = cls() - - return obj - - def sgp4_update(self, jd, fr): - """ - For a julian date (jd) and its fractional representation (fr), - propagate Satrec using sgp4 - """ - - sgp4_update() + return month, day, int(hour), int(minute), second def _conv_year(s): """Interpret a two-digit year string.""" @@ -99,14 +106,16 @@ def _parse_float(s): """ return float(s[0] + '.' + s[1:6] + 'e' + s[6:8]) -class TLE: +class Satrec: """ + Satellite record object + Includes parameters, constants that are commonly used across the sgp4 logical flow + + In this implementation, built from TLE data + Two line-elements (TLEs) are unpacked from both given and propagated data. This implementation uses Keplerian orbital parameters - A two-line element set (TLE) is a data format encoding a list of orbital - elements of an Earth-orbiting object for a given point in time, the epoch. - All the attributes parsed from the TLE are expressed in the same units that are used in the TLE format. @@ -122,10 +131,10 @@ class TLE: Year of the epoch. :float epoch_day: Day of the year plus fraction of the day. - :float dn_o2: - First time derivative of the mean motion divided by 2. - :float ddn_o6: - Second time derivative of the mean motion divided by 6. + :float dn: + First time derivative of the mean motion (divided by 2 in TLE) + :float ddn: + Second time derivative of the mean motion (divided by 6 in TLE). :float bstar: BSTAR coefficient (https://en.wikipedia.org/wiki/BSTAR). :int set_num: @@ -146,16 +155,23 @@ class TLE: Revolution number. """ + # Error codes + ECCENTRICITY = 1 # eccentricity is not within 0-1 + MOTION = 2 # error in propagating mean motion + SEMIRECT = 4 # apoapsis, periapsis characteristics error + DECAY = 6 # orbit has decayed + def __init__(self, name:str, # ID parameters, Line 1 norad:str, classification:str, int_desig:str, # time (derivative) parameters, line 1 - epoch_year:int, epoch_day:float, dn_o2:float, ddn_o6:float, bstar:float, set_num:int, + epoch_year:int, epoch_day:float, dn:float, ddn:float, bstar:float, ephtype: str, set_num:int, # keplerian parameters, line 2 - inc:float, raan:float, ecc:float, argp:float, M:float, n:float, rev_num:int): - # Oh my lord prepare for absolute misery on earth + inc:float, raan:float, ecc:float, argp:float, M:float, n:float, rev_num:int, + # for the purposes of keeping the tle around as future-proofing + tle_str:str ): - self.name = str.strip(name) + self.name = str.strip(name) self.norad = str.strip(norad) self.classification = classification @@ -163,21 +179,24 @@ def __init__(self, name:str, self.epoch_year = _conv_year(epoch_year) self.epoch_day = epoch_day - self.dn_o2 = dn_o2 - self.ddn_o6 = ddn_o6 + self.dn = dn + self.ddn = ddn self.bstar = bstar + self.ephtype = ephtype self.set_num = int(set_num) self.inc = inc self.raan = raan self.ecc = ecc self.argp = argp - self.M = M - self.n = n + self.M = M + self.n = n # mean motion self.rev_num = int(rev_num) + + self.tle_str = tle_str @classmethod - def from_lines(cls, name, line1, line2): + def from_tle_lines(cls, name, line1, line2): """Parse a TLE from its constituent lines. All the attributes parsed from the TLE are expressed in the same units that @@ -186,13 +205,14 @@ def from_lines(cls, name, line1, line2): return cls( name=name, norad=line1[2:7], - classification=line1[7], + classification=line1[7] or 'U', int_desig=line1[9:17], epoch_year=line1[18:20], epoch_day=float(line1[20:32]), - dn_o2=float(line1[33:43]), - ddn_o6=_parse_float(line1[44:52]), + dn=float(line1[33:43]), + ddn=_parse_float(line1[44:52]), bstar=_parse_float(line1[53:61]), + ephtype = line1[62], set_num=line1[64:68], inc=float(line2[8:16]), raan=float(line2[17:25]), @@ -203,17 +223,60 @@ def from_lines(cls, name, line1, line2): rev_num=line2[63:68]) @classmethod - def from_file(cls, filename): + def from_tle_file(cls, filename): """Load TLE from a file.""" if isinstance(filename, str): with open(filename) as fp: return [cls.from_lines(*fp.readlines[:2])] @classmethod - def from_str(cls, string): + def from_tle_str(cls, string): """Load TLE from a string.""" return [cls.from_lines(*string.split('\n')[:2])] + @classmethod + def sgp4_init(cls, tle: Satrec): + """ + Creates a satrec object specifically modified to be used in sgp4. + + Changes units, activates certain new parameters, etc. + """ + + self = tle + + # constants for unit change + deg2rad = np.pi / 180.0; # 0.0174532925199433 + xpdotp = 1440.0 / (2.0 *np.pi); # 229.1831180523293 + + # ---- convert to sgp4 units ---- + self.n = self.n / xpdotp + self.dn = self.dn / (xpdotp*1440.0) + self.ddn= self.ddn / (xpdotp*1440.0*1440) + + # ---- find standard orbital elements ---- + self.inc = self.inc * deg2rad + self.raan = self.raan * deg2rad + self.argp = self.argp * deg2rad + self.M = self.M * deg2rad + + yr = self.epoch_year + + # Build Julian Date + if yr < 57: + year = yr + 2000 + else: + year = yr + 1900 + + mon,day,hr,minute,sec = _days2mdhms(year, self.epoch_day) + self.jdsatepoch = _jday(year,mon,day,hr,minute,sec); + epoch0 = 2433281.5 + + sgp4_init(self, self.set_num, self.jdsatepoch - epoch0, self.bstar, + self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, + self.raan) + + return self + def to_array(self): """ Return 2D array of TLE values @@ -226,20 +289,32 @@ def to_array(self): name: [0,0] norad: [1,0] classification: [1,1] int_desig: [1,2] epoch_year: [1,3] day: [1,4] - dn/2: [1,5] ddn/6: [1,6] bstar: [1,7] set_num: [1,8] + dn: [1,5] ddn: [1,6] bstar: [1,7] set_num: [1,8] inclination: [2,0] RAAN: [2,1] eccentricity: [2,2] arg_perigee: [2,3] Mean Anomaly: [2,4] n: [2,5] rev_num: [2,6] """ - TLE() + return [ [self.name], # Line 0 [self.norad, self.classification, self.int_desig, # line 1 ID - self.epoch_year, self.epoch_day, self.dn_o2, self.ddn_o6, self.bstar, self.set_num], # line 1 time-derivative + self.epoch_year, self.epoch_day, self.dn, self.ddn, self.bstar, self.set_num], # line 1 time-derivative [self.inc, self.raan, self.ecc, self.argp, self.M, self.n, self.rev_num] # line 2 orbital params ] - def to_sgp4_params(self): - """Return a formatted Satrec object that can immediately used in SGP4""" + + def sgp4_update(self, jd, fr): + """ + For a julian date (jd) and its fractional representation (fr), + propagate Satrec using sgp4 + """ + + tsince = ((jd - self.sat_epoch) * minutes_per_day + + (fr - self.sat_epochF) * minutes_per_day) + r, v = sgp4_update(self, tsince) + + return self.error, r, v - return Satrec.from_tle_array(self.to_array()) - \ No newline at end of file + def error_message(self): + if self.error == self.MOTION: + return ('mean motion {0:f} is less than zero' + .format(self.n)) diff --git a/unit_tests/lib/sgp4_test.py b/unit_tests/lib/sgp4_test.py new file mode 100644 index 0000000..06b6a96 --- /dev/null +++ b/unit_tests/lib/sgp4_test.py @@ -0,0 +1,11 @@ +import unittest +import math +from tle import Satrec + +try: + import ulab.numpy as np # For CircuitPython +except ImportError: + import numpy as np # For GitHub Actions / PC testing + +class PropagatorTest(unittest.TestCase): + pass From 9bf2b7f7cb33c122a8b218a211914da36284c3c1 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 4 May 2026 16:06:49 -0700 Subject: [PATCH 06/19] sgp4 basic test case (propagate to may 3) implemented, static analyser fixes --- src/lib/sgp4.py | 480 ++++++++++++++++++------------------ src/lib/tle.py | 125 +++++----- unit_tests/lib/sgp4_test.py | 27 +- 3 files changed, 335 insertions(+), 297 deletions(-) diff --git a/src/lib/sgp4.py b/src/lib/sgp4.py index 305514e..71dbd1f 100644 --- a/src/lib/sgp4.py +++ b/src/lib/sgp4.py @@ -30,7 +30,7 @@ def _gstime(jdut1): return temp def _initl(xke, j2, - ecco, epoch, inclo, no): + ecc, epoch, inclo, no): # ----------------------- earth constants ---------------------- # sgp4fix identify constants and allow alternate values # only xke and j2 are used here so pass them in directly @@ -38,7 +38,7 @@ def _initl(xke, j2, x2o3 = 2.0 / 3.0; # ------------- calculate auxillary epoch quantities ---------- - eccsq = ecco * ecco; + eccsq = ecc * ecc; omeosq = 1.0 - eccsq; rteosq = np.sqrt(omeosq); cosio = np.cos(inclo); @@ -59,7 +59,7 @@ def _initl(xke, j2, con42 = 1.0 - 5.0 * cosio2 con41 = -con42-cosio2-cosio2 posq = po * po - rp = ao * (1.0 - ecco) + rp = ao * (1.0 - ecc) gsto = _gstime(epoch + 2433281.5) @@ -85,60 +85,60 @@ def sgp4_update(satrec, tsince): satrec.t = tsince # -- Update for secular gravity and atmospheric drag - xmdf = satrec.mo + satrec.mdot * satrec.t; - argpdf = satrec.argpo + satrec.argpdot * satrec.t; - nodedf = satrec.nodeo + satrec.nodedot * satrec.t; + xmdf = satrec.mo + satrec.mdot * satrec.t + argpdf = satrec.argp + satrec.argpdot * satrec.t + nodedf = satrec.raan + satrec.nodedot * satrec.t argpm = argpdf mm = xmdf - t2 = satrec.t * satrec.t; - nodem = nodedf + satrec.nodecf * t2; - tempa = 1.0 - satrec.cc1 * satrec.t; - tempe = satrec.bstar * satrec.cc4 * satrec.t; - templ = satrec.t2cof * t2; + t2 = satrec.t * satrec.t + nodem = nodedf + satrec.nodecf * t2 + tempa = 1.0 - satrec.cc1 * satrec.t + tempe = satrec.bstar * satrec.cc4 * satrec.t + templ = satrec.t2cof * t2 if satrec.isimp != 1: - delomg = satrec.omgcof * satrec.t; + delomg = satrec.omgcof * satrec.t # sgp4fix use mutliply for speed instead of pow - delmtemp = 1.0 + satrec.eta * np.cos(xmdf); + delmtemp = 1.0 + satrec.eta * np.cos(xmdf) delm = satrec.xmcof * \ (delmtemp * delmtemp * delmtemp - - satrec.delmo); - temp = delomg + delm; - mm = xmdf + temp; - argpm = argpdf - temp; - t3 = t2 * satrec.t; - t4 = t3 * satrec.t; + satrec.delmo) + temp = delomg + delm + mm = xmdf + temp + argpm = argpdf - temp + t3 = t2 * satrec.t + t4 = t3 * satrec.t tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ - satrec.d4 * t4; + satrec.d4 * t4 tempe = tempe + satrec.bstar * satrec.cc5 * (np.sin(mm) - - satrec.sinmao); + satrec.sinmao) templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + - satrec.t * satrec.t5cof); + satrec.t * satrec.t5cof) - nm = satrec.no_unkozai; - em = satrec.ecco; - inclm = satrec.inclo; + nm = satrec.n + em = satrec.ecc + inclm = satrec.inclo if nm <= 0.0: satrec.error = satrec.MOTION - return False, False; + return False, False - am = pow((satrec.xke / nm),x2o3) * tempa * tempa; - nm = satrec.xke / pow(am, 1.5); - em = em - tempe; + am = pow((satrec.xke / nm),x2o3) * tempa * tempa + nm = satrec.xke / pow(am, 1.5) + em = em - tempe - if em >= 1.0 or em < -0.001: + if em >= 1.0 or em < -0.001: satrec.error = satrec.ECCENTRICITY - - return False, False; + + return False, False if em < 1.0e-6: - em = 1.0e-6; - mm = mm + satrec.no_unkozai * templ; - xlm = mm + argpm + nodem; - emsq = em * em; - temp = 1.0 - emsq; + em = 1.0e-6 + mm = mm + satrec.n * templ + xlm = mm + argpm + nodem + emsq = em * em + temp = 1.0 - emsq nodem = nodem % tau if nodem >= 0.0 else -(-nodem % tau) argpm = argpm % tau @@ -158,46 +158,46 @@ def sgp4_update(satrec, tsince): cosim = np.cos(inclm) # -------------------- add lunar-solar periodics -------------- - ep = em; - xincp = inclm; - argpp = argpm; - nodep = nodem; - mp = mm; - sinip = sinim; - cosip = cosim; + ep = em + xincp = inclm + argpp = argpm + nodep = nodem + mp = mm + sinip = sinim + cosip = cosim # -------------------- long period periodics ------------------ - axnl = ep * np.cos(argpp); - temp = 1.0 / (am * (1.0 - ep * ep)); - aynl = ep* np.sin(argpp) + temp * satrec.aycof; - xl = mp + argpp + nodep + temp * satrec.xlcof * axnl; + axnl = ep * np.cos(argpp) + temp = 1.0 / (am * (1.0 - ep * ep)) + aynl = ep* np.sin(argpp) + temp * satrec.aycof + xl = mp + argpp + nodep + temp * satrec.xlcof * axnl # --------------------- solve kepler's equation --------------- u = (xl - nodep) % tau - eo1 = u; - tem5 = 9999.9; - ktr = 1; + eo1 = u + tem5 = 9999.9 + ktr = 1 while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: - sineo1 = np.sin(eo1); - coseo1 = np.cos(eo1); - tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl; - tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5; + sineo1 = np.sin(eo1) + coseo1 = np.cos(eo1) + tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl + tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5 if np.fabs(tem5) >= 0.95: - tem5 = 0.95 if tem5 > 0.0 else -0.95; - eo1 = eo1 + tem5; - ktr = ktr + 1; + tem5 = 0.95 if tem5 > 0.0 else -0.95 + eo1 = eo1 + tem5 + ktr = ktr + 1 # ------------- short period preliminary quantities ----------- - ecose = axnl*coseo1 + aynl*sineo1; - esine = axnl*sineo1 - aynl*coseo1; - el2 = axnl*axnl + aynl*aynl; - pl = am*(1.0-el2); + ecose = axnl*coseo1 + aynl*sineo1 + esine = axnl*sineo1 - aynl*coseo1 + el2 = axnl*axnl + aynl*aynl + pl = am*(1.0-el2) if pl < 0.0: satrec.error = satrec.SEMIRECT - - return False, False; + + return False, False else: @@ -217,29 +217,29 @@ def sgp4_update(satrec, tsince): # -------------- update for short period periodics ------------ mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ - 0.5 * temp1 * satrec.x1mth2 * cos2u; - su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u; - xnode = nodep + 1.5 * temp2 * cosip * sin2u; - xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u; - mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke; + 0.5 * temp1 * satrec.x1mth2 * cos2u + su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u + xnode = nodep + 1.5 * temp2 * cosip * sin2u + xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u + mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + - 1.5 * satrec.con41) / satrec.xke; + 1.5 * satrec.con41) / satrec.xke # --------------------- orientation vectors ------------------- - sinsu = np.sin(su); - cossu = np.cos(su); - snod = np.sin(xnode); - cnod = np.cos(xnode); - sini = np.sin(xinc); - cosi = np.cos(xinc); - xmx = -snod * cosi; - xmy = cnod * cosi; - ux = xmx * sinsu + cnod * cossu; - uy = xmy * sinsu + snod * cossu; - uz = sini * sinsu; - vx = xmx * cossu - cnod * sinsu; - vy = xmy * cossu - snod * sinsu; - vz = sini * cossu; + sinsu = np.sin(su) + cossu = np.cos(su) + snod = np.sin(xnode) + cnod = np.cos(xnode) + sini = np.sin(xinc) + cosi = np.cos(xinc) + xmx = -snod * cosi + xmy = cnod * cosi + ux = xmx * sinsu + cnod * cossu + uy = xmy * sinsu + snod * cossu + uz = sini * sinsu + vx = xmx * cossu - cnod * sinsu + vy = xmy * cossu - snod * sinsu + vz = sini * cossu # --------- position and velocity (in km and km/sec) ---------- _mr = mrt * satrec.radiusearthkm @@ -252,158 +252,162 @@ def sgp4_update(satrec, tsince): satrec.error = satrec.DECAY return r, v - + def sgp4_init(satrec, satn, epoch, - bstar, ndot, nddot, ecco, argpo, - inclo, mo, - nodeo - ): - - temp4 = 1.5e-12 - - # Near Earth Variables - satrec.isimp = 0; satrec.aycof = 0.0 - satrec.con41 = 0.0; satrec.cc1 = 0.0; satrec.cc4 = 0.0 - satrec.cc5 = 0.0; satrec.d2 = 0.0; satrec.d3 = 0.0 - satrec.d4 = 0.0; satrec.delmo = 0.0; satrec.eta = 0.0 - satrec.argpdot = 0.0; satrec.omgcof = 0.0; satrec.sinmao = 0.0 - satrec.t = 0.0; satrec.t2cof = 0.0; satrec.t3cof = 0.0 - satrec.t4cof = 0.0; satrec.t5cof = 0.0; satrec.x1mth2 = 0.0 - satrec.x7thm1 = 0.0; satrec.mdot = 0.0; satrec.nodedot = 0.0 - satrec.xlcof = 0.0; satrec.xmcof = 0.0; satrec.nodecf = 0.0 - - # Earth Constants - satrec.mu = 398600.5; # in km3 / s2 - satrec.radiusearthkm = 6378.137 # km - satrec.xke = 0.07436685317 - satrec.tumin = 13.44685108 - satrec.j2 = 0.00108262998905 - satrec.j3 = -0.00000253215306 - satrec.j4 = -0.00000161098761 - satrec.j3oj2 = -0.002338890559 - - ss = 1.012229276 - qzms2ttemp = 0.00658499496 - qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp; - x2o3 = 2.0 / 3.0 - - # -- initialisation markers - satrec.t = 0.0 - - # -- - satrec.satnum_str = satn - satrec.classification = 'U' - - # -- - satrec.bstar = bstar - satrec.ndot = ndot - satrec.nddot = nddot - satrec.ecco = ecco - satrec.argpo = argpo - satrec.inclo = inclo - satrec.mo = mo - satrec.nodeo = nodeo - - # single averaged mean elements - satrec.am = 0.0 - satrec.em = 0.0 - satrec.im = 0.0 - satrec.Om = 0.0 - satrec.mm = 0.0 - satrec.nm = 0.0 - - satrec.error = 0 - - # -- - ( - satrec.no_unkozai, - ao, satrec.con41, con42, cosio, - cosio2, omeosq, posq, - rp, rteosq,sinio , satrec.gsto, - ) = _initl( - satrec.xke, satrec.j2, satrec.ecco, epoch, satrec.inclo - ) - satrec.a = pow( satrec.no_unkozai*satrec.tumin , (-2.0/3.0) ); - satrec.alta = satrec.a*(1.0 + satrec.ecco) - 1.0; - satrec.altp = satrec.a*(1.0 - satrec.ecco) - 1.0; - - if omeosq >= 0.0 or satrec.no_unkozai >= 0.0: - satrec.isimp = 0 - if rp < 220.0 / satrec.radiusearthkm + 1.0: - satrec.isimp = 1 - sfour = ss - qzms24 = qzms2t - pinvsq = 1.0 / posq; - - tsi = 1.0 / (ao - sfour); - satrec.eta = ao * satrec.ecco * tsi; - etasq = satrec.eta * satrec.eta; - eeta = satrec.ecco * satrec.eta; - psisq = np.fabs(1.0 - etasq); - coef = qzms24 * pow(tsi, 4.0); - coef1 = coef / pow(psisq, 3.5); - cc2 = coef1 * satrec.no_unkozai * (ao * (1.0 + 1.5 * etasq + eeta * - (4.0 + etasq)) + 0.375 * satrec.j2 * tsi / psisq * satrec.con41 * - (8.0 + 3.0 * etasq * (8.0 + etasq))); - satrec.cc1 = satrec.bstar * cc2; - cc3 = 0.0; - if satrec.ecco > 1.0e-4: - cc3 = -2.0 * coef * tsi * satrec.j3oj2 * satrec.no_unkozai * sinio / satrec.ecco; - satrec.x1mth2 = 1.0 - cosio2; - satrec.cc4 = 2.0* satrec.no_unkozai * coef1 * ao * omeosq * \ - (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco * - (0.5 + 2.0 * etasq) - satrec.j2 * tsi / (ao * psisq) * - (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * - (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * - (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * satrec.argpo))); - satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * - (etasq + eeta) + eeta * etasq); - cosio4 = cosio2 * cosio2; - temp1 = 1.5 * satrec.j2 * pinvsq * satrec.no_unkozai; - temp2 = 0.5 * temp1 * satrec.j2 * pinvsq; - temp3 = -0.46875 * satrec.j4 * pinvsq * pinvsq * satrec.no_unkozai; - satrec.mdot = satrec.no_unkozai + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * \ - temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); - satrec.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * - (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + - temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)); - xhdot1 = -temp1 * cosio; - satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + - 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; - satrec.omgcof = satrec.bstar * cc3 * np.cos(satrec.argpo); - satrec.xmcof = 0.0; - if satrec.ecco > 1.0e-4: - satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta; - satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1; - satrec.t2cof = 1.5 * satrec.cc1; - - if np.fabs(cosio+1.0) > 1.5e-12: - satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); - else: - satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4; - satrec.aycof = -0.5 * satrec.j3oj2 * sinio; - - delmotemp = 1.0 + satrec.eta * np.cos(satrec.mo); - satrec.delmo = delmotemp * delmotemp * delmotemp; - satrec.sinmao = np.sin(satrec.mo); - satrec.x7thm1 = 7.0 * cosio2 - 1.0; - - if satrec.isimp != 1: - cc1sq = satrec.cc1 * satrec.cc1; - satrec.d2 = 4.0 * ao * tsi * cc1sq; - temp = satrec.d2 * tsi * satrec.cc1 / 3.0; - satrec.d3 = (17.0 * ao + sfour) * temp; - satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ - satrec.cc1; - satrec.t3cof = satrec.d2 + 2.0 * cc1sq; - satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * - (12.0 * satrec.d2 + 10.0 * cc1sq)); - satrec.t5cof = 0.2 * (3.0 * satrec.d4 + - 12.0 * satrec.cc1 * satrec.d3 + - 6.0 * satrec.d2 * satrec.d2 + - 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)); - - # propagate to 0 - sgp4_update(satrec, 0) - - return True + bstar, dn, ddn, ecc, argp, + inclo, mo, + raan + ): + """ + (Math Heavy) initialisation of sgp4 which instantiates all important + variables in the satrec object + """ + + temp4 = 1.5e-12 + + # Near Earth Variables + satrec.isimp = 0; satrec.aycof = 0.0 + satrec.con41 = 0.0; satrec.cc1 = 0.0; satrec.cc4 = 0.0 + satrec.cc5 = 0.0; satrec.d2 = 0.0; satrec.d3 = 0.0 + satrec.d4 = 0.0; satrec.delmo = 0.0; satrec.eta = 0.0 + satrec.argpdot = 0.0; satrec.omgcof = 0.0; satrec.sinmao = 0.0 + satrec.t = 0.0; satrec.t2cof = 0.0; satrec.t3cof = 0.0 + satrec.t4cof = 0.0; satrec.t5cof = 0.0; satrec.x1mth2 = 0.0 + satrec.x7thm1 = 0.0; satrec.mdot = 0.0; satrec.nodedot = 0.0 + satrec.xlcof = 0.0; satrec.xmcof = 0.0; satrec.nodecf = 0.0 + + # Earth Constants + satrec.mu = 398600.5 # in km3 / s2 + satrec.radiusearthkm = 6378.137 # km + satrec.xke = 0.07436685317 + satrec.tumin = 13.44685108 + satrec.j2 = 0.00108262998905 + satrec.j3 = -0.00000253215306 + satrec.j4 = -0.00000161098761 + satrec.j3oj2 = -0.002338890559 + + ss = 1.012229276 + qzms2ttemp = 0.00658499496 + qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp + x2o3 = 2.0 / 3.0 + + # -- initialisation markers + satrec.t = 0.0 + + # -- + satrec.norad = satn + satrec.classification = 'U' + + # -- + satrec.bstar = bstar + satrec.dn = dn + satrec.ddn = ddn + satrec.ecc = ecc + satrec.argp = argp + satrec.inclo = inclo + satrec.mo = mo + satrec.raan = raan + + # single averaged mean elements + satrec.am = 0.0 + satrec.em = 0.0 + satrec.im = 0.0 + satrec.Om = 0.0 + satrec.mm = 0.0 + satrec.nm = 0.0 + + satrec.error = 0 + + # -- + ( + satrec.n, + ao, satrec.con41, con42, cosio, + cosio2, omeosq, posq, + rp, rteosq,sinio , satrec.gsto, + ) = _initl( + satrec.xke, satrec.j2, satrec.ecc, epoch, satrec.inclo, satrec.n + ) + satrec.a = pow( satrec.n*satrec.tumin , (-2.0/3.0) ) + satrec.alta = satrec.a*(1.0 + satrec.ecc) - 1.0 + satrec.altp = satrec.a*(1.0 - satrec.ecc) - 1.0 + + if omeosq >= 0.0 or satrec.n >= 0.0: + satrec.isimp = 0 + if rp < 220.0 / satrec.radiusearthkm + 1.0: + satrec.isimp = 1 + sfour = ss + qzms24 = qzms2t + pinvsq = 1.0 / posq + + tsi = 1.0 / (ao - sfour) + satrec.eta = ao * satrec.ecc * tsi + etasq = satrec.eta * satrec.eta + eeta = satrec.ecc * satrec.eta + psisq = np.fabs(1.0 - etasq) + coef = qzms24 * pow(tsi, 4.0); + coef1 = coef / pow(psisq, 3.5) + cc2 = coef1 * satrec.n * (ao * (1.0 + 1.5 * etasq + eeta * + (4.0 + etasq)) + 0.375 * satrec.j2 * tsi / psisq * satrec.con41 * + (8.0 + 3.0 * etasq * (8.0 + etasq))) + satrec.cc1 = satrec.bstar * cc2 + cc3 = 0.0 + if satrec.ecc > 1.0e-4: + cc3 = -2.0 * coef * tsi * satrec.j3oj2 * satrec.n * sinio / satrec.ecc + satrec.x1mth2 = 1.0 - cosio2 + satrec.cc4 = 2.0* satrec.n * coef1 * ao * omeosq * \ + (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecc * + (0.5 + 2.0 * etasq) - satrec.j2 * tsi / (ao * psisq) * + (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * + (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * + (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * satrec.argp))) + satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * + (etasq + eeta) + eeta * etasq) + cosio4 = cosio2 * cosio2 + temp1 = 1.5 * satrec.j2 * pinvsq * satrec.n + temp2 = 0.5 * temp1 * satrec.j2 * pinvsq + temp3 = -0.46875 * satrec.j4 * pinvsq * pinvsq * satrec.n + satrec.mdot = satrec.n + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * \ + temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4) + satrec.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * + (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + + temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)) + xhdot1 = -temp1 * cosio + satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio + satrec.omgcof = satrec.bstar * cc3 * np.cos(satrec.argp) + satrec.xmcof = 0.0 + if satrec.ecc > 1.0e-4: + satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta + satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1 + satrec.t2cof = 1.5 * satrec.cc1 + + if np.fabs(cosio+1.0) > 1.5e-12: + satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio) + else: + satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4 + satrec.aycof = -0.5 * satrec.j3oj2 * sinio + + delmotemp = 1.0 + satrec.eta * np.cos(satrec.mo) + satrec.delmo = delmotemp * delmotemp * delmotemp + satrec.sinmao = np.sin(satrec.mo) + satrec.x7thm1 = 7.0 * cosio2 - 1.0 + + if satrec.isimp != 1: + cc1sq = satrec.cc1 * satrec.cc1 + satrec.d2 = 4.0 * ao * tsi * cc1sq + temp = satrec.d2 * tsi * satrec.cc1 / 3.0 + satrec.d3 = (17.0 * ao + sfour) * temp + satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ + satrec.cc1 + satrec.t3cof = satrec.d2 + 2.0 * cc1sq + satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * + (12.0 * satrec.d2 + 10.0 * cc1sq)) + satrec.t5cof = 0.2 * (3.0 * satrec.d4 + + 12.0 * satrec.cc1 * satrec.d3 + + 6.0 * satrec.d2 * satrec.d2 + + 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)) + + # propagate to 0 + sgp4_update(satrec, 0) + + return True diff --git a/src/lib/tle.py b/src/lib/tle.py index 9b08a88..045bb63 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -11,8 +11,8 @@ except ImportError: import numpy as np # For GitHub Actions / PC testing -minutes_per_day = 1440 -epoch0 = 2433281.5 # jan 0 1950 +MIN_PER_DAY = 1440 +EPOCH0 = 2433281.5 # jan 0 1950 def _day_of_year_to_month_day(day_of_year, is_leap): """Core logic for turning days into months, for easy testing.""" @@ -24,30 +24,6 @@ def _day_of_year_to_month_day(day_of_year, is_leap): day += 1 return month, day -def _jday(year, mon, day, hr, minute, sec): - """Return two floats (jd, fr) that, when added, produce the specified Julian date - jd (Julian Date) and fr (Fractional) components - - >>> jd, fr = jday(2020, 2, 11, 13, 57, 0) - >>> jd - 2458890.5 - >>> fr - 0.58125 - - Note the first float, which gives the moment of midnight that - commences the given calendar date, always ends in - ``.5`` because Julian dates begin and end at noon. This made - Julian dates more convenient for astronomers in Europe, by making - the whole night belong to a single Julian date. - """ - jd = (367.0 * year - - 7 * (year + ((mon + 9) // 12.0)) * 0.25 // 1.0 - + 275 * mon / 9.0 // 1.0 - + day - + 1721013.5) - fr = (sec + minute * 60.0 + hr * 3600.0) / 86400.0; - return jd, fr - def _days2mdhms(year, days, round_to_microsecond=6): """Convert a float point number of days into the year into date and time. @@ -106,6 +82,32 @@ def _parse_float(s): """ return float(s[0] + '.' + s[1:6] + 'e' + s[6:8]) +def _sgp4_jday(year, mon, day, hr, minute, sec): + """ + Converts jdsatepoch into a compatible number for sgp4 + """ + + return (367.0 * year - + 7.0 * (year + ((mon + 9.0) // 12.0)) * 0.25 // 1.0 + + 275.0 * mon // 9.0 + + day + 1721013.5 + + ((sec / 60.0 + minute) / 60.0 + hr) / 24.0 + ) + +def jday(year, mon, day, hr, minute, sec): + """ + From a date, return a Julian date in its date + fractional form + + Used to build the jd, fr parameters in the function call of sgp4_update + """ + jd = (367.0 * year + - 7 * (year + ((mon + 9) // 12.0)) * 0.25 // 1.0 + + 275 * mon / 9.0 // 1.0 + + day + + 1721013.5) + fr = (sec + minute * 60.0 + hr * 3600.0) / 86400.0 + return jd, fr + class Satrec: """ Satellite record object @@ -147,7 +149,7 @@ class Satrec: Eccentricity. :float argp: Argument of perigee. - :float M: + :float mo: Mean anomaly. :float n: Mean motion. @@ -163,15 +165,16 @@ class Satrec: def __init__(self, name:str, # ID parameters, Line 1 - norad:str, classification:str, int_desig:str, + norad:str, classification:str, int_desig:str, # time (derivative) parameters, line 1 - epoch_year:int, epoch_day:float, dn:float, ddn:float, bstar:float, ephtype: str, set_num:int, + epoch_year:int, epoch_day:float, dn:float, ddn:float, bstar:float, + ephtype: str, set_num:int, # keplerian parameters, line 2 - inc:float, raan:float, ecc:float, argp:float, M:float, n:float, rev_num:int, + inc:float, raan:float, ecc:float, argp:float, mo:float, n:float, rev_num:int, # for the purposes of keeping the tle around as future-proofing tle_str:str ): - - self.name = str.strip(name) + + self.name = str.strip(name) self.norad = str.strip(norad) self.classification = classification @@ -189,12 +192,16 @@ def __init__(self, name:str, self.raan = raan self.ecc = ecc self.argp = argp - self.M = M + self.mo = mo self.n = n # mean motion self.rev_num = int(rev_num) self.tle_str = tle_str - + + # for sgp + self.jdsatepoch = 0 + self.jdsatepoch_f = 0 + @classmethod def from_tle_lines(cls, name, line1, line2): """Parse a TLE from its constituent lines. @@ -218,21 +225,22 @@ def from_tle_lines(cls, name, line1, line2): raan=float(line2[17:25]), ecc=_parse_decimal(line2[26:33]), argp=float(line2[34:42]), - M=float(line2[43:51]), + mo=float(line2[43:51]), n=float(line2[52:63]), - rev_num=line2[63:68]) + rev_num=line2[63:68], + tle_str=name+line1+line2) @classmethod def from_tle_file(cls, filename): """Load TLE from a file.""" if isinstance(filename, str): with open(filename) as fp: - return [cls.from_lines(*fp.readlines[:2])] + return cls.from_tle_lines(*fp.readlines[:2]) @classmethod def from_tle_str(cls, string): """Load TLE from a string.""" - return [cls.from_lines(*string.split('\n')[:2])] + return cls.from_tle_lines(*string.split('\n')[:3]) @classmethod def sgp4_init(cls, tle: Satrec): @@ -245,8 +253,8 @@ def sgp4_init(cls, tle: Satrec): self = tle # constants for unit change - deg2rad = np.pi / 180.0; # 0.0174532925199433 - xpdotp = 1440.0 / (2.0 *np.pi); # 229.1831180523293 + deg2rad = np.pi / 180.0 # 0.0174532925199433 + xpdotp = 1440.0 / (2.0 *np.pi) # 229.1831180523293 # ---- convert to sgp4 units ---- self.n = self.n / xpdotp @@ -257,9 +265,11 @@ def sgp4_init(cls, tle: Satrec): self.inc = self.inc * deg2rad self.raan = self.raan * deg2rad self.argp = self.argp * deg2rad - self.M = self.M * deg2rad + self.mo = self.mo * deg2rad yr = self.epoch_year + _, fraction = divmod(self.epoch_day, 1.0) + self.jdsatepoch_f = round(fraction, 8) # exact number of digits in TLE # Build Julian Date if yr < 57: @@ -268,13 +278,12 @@ def sgp4_init(cls, tle: Satrec): year = yr + 1900 mon,day,hr,minute,sec = _days2mdhms(year, self.epoch_day) - self.jdsatepoch = _jday(year,mon,day,hr,minute,sec); - epoch0 = 2433281.5 + self.jdsatepoch = _sgp4_jday(year,mon,day,hr,minute,sec) - sgp4_init(self, self.set_num, self.jdsatepoch - epoch0, self.bstar, - self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, + sgp4_init(self, self.set_num, self.jdsatepoch-EPOCH0, self.bstar, + self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, self.raan) - + return self def to_array(self): @@ -291,30 +300,34 @@ def to_array(self): norad: [1,0] classification: [1,1] int_desig: [1,2] epoch_year: [1,3] day: [1,4] dn: [1,5] ddn: [1,6] bstar: [1,7] set_num: [1,8] - inclination: [2,0] RAAN: [2,1] eccentricity: [2,2] arg_perigee: [2,3] Mean Anomaly: [2,4] n: [2,5] rev_num: [2,6] + inclination: [2,0] RAAN: [2,1] eccentricity: [2,2] arg_perigee: [2,3] Mean Anomaly: [2,4] + n: [2,5] rev_num: [2,6] """ - + return [ [self.name], # Line 0 [self.norad, self.classification, self.int_desig, # line 1 ID - self.epoch_year, self.epoch_day, self.dn, self.ddn, self.bstar, self.set_num], # line 1 time-derivative - [self.inc, self.raan, self.ecc, self.argp, self.M, self.n, self.rev_num] # line 2 orbital params + # line 1 time-derivative + self.epoch_year, self.epoch_day, self.dn, self.ddn, self.bstar, self.set_num], + # line 2 orbital params + [self.inc, self.raan, self.ecc, self.argp, self.mo, self.n, self.rev_num] ] - - + def sgp4_update(self, jd, fr): """ For a julian date (jd) and its fractional representation (fr), propagate Satrec using sgp4 """ - tsince = ((jd - self.sat_epoch) * minutes_per_day + - (fr - self.sat_epochF) * minutes_per_day) + tsince = ((jd - self.jdsatepoch) * MIN_PER_DAY + + (fr - self.jdsatepoch_f) * MIN_PER_DAY) r, v = sgp4_update(self, tsince) return self.error, r, v def error_message(self): + """ + Return error message from current self.error (when polled) + """ if self.error == self.MOTION: - return ('mean motion {0:f} is less than zero' - .format(self.n)) + return (f'mean motion {0:f} is less than zero').format(self.n) diff --git a/unit_tests/lib/sgp4_test.py b/unit_tests/lib/sgp4_test.py index 06b6a96..f180ef5 100644 --- a/unit_tests/lib/sgp4_test.py +++ b/unit_tests/lib/sgp4_test.py @@ -1,11 +1,32 @@ import unittest -import math -from tle import Satrec +from tle import Satrec, jday try: import ulab.numpy as np # For CircuitPython except ImportError: import numpy as np # For GitHub Actions / PC testing + +ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" + class PropagatorTest(unittest.TestCase): - pass + + def propagate_to_may_3(self): + sat: Satrec = Satrec.from_tle_str( + ISS_TLE + ) + + sgp4_obj = Satrec.sgp4_init(sat) + error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) + + self.assertEqual(error, 0) # if no error, + + tol = 3 # to 3 decimal places + self.assertAlmostEqual(r[0], 4698.782358, tol) + self.assertAlmostEqual(r[1], -3867.014434, tol) + self.assertAlmostEqual(r[2], 3028.549126, tol) + + tol = 6 # to 6 decimal places, these need to be more accurate + self.assertAlmostEqual(v[0], 5.281325344, tol) + self.assertAlmostEqual(v[1], 2.530170911, tol) + self.assertAlmostEqual(v[2], -4.936649541, tol) From 55ad5815eda16bf73ab40f521ca95d830ccb7763 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 4 May 2026 18:17:23 -0700 Subject: [PATCH 07/19] actually running the unit test --- artifacts/adcs_breakout_board_sim/include.json | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/artifacts/adcs_breakout_board_sim/include.json b/artifacts/adcs_breakout_board_sim/include.json index 05bdd60..0273495 100644 --- a/artifacts/adcs_breakout_board_sim/include.json +++ b/artifacts/adcs_breakout_board_sim/include.json @@ -14,7 +14,8 @@ "unit_tests": [ "lib/pin_manager_test.py:pin_manager_test.py", "lib/custom_module_mocking.py:custom_module_mocking.py", - "lib/quaternion_test.py:quaternion_test.py" + "lib/quaternion_test.py:quaternion_test.py", + "lib/sgp4_test.py:sgp4_test.py" ], "submodules":[ "Adafruit_CircuitPython_Ticks/adafruit_ticks.py:adafruit_ticks.py", From a2c21f4ece390b47644e28a8aaaadce907785ca4 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 4 May 2026 18:47:56 -0700 Subject: [PATCH 08/19] forgot the new includes --- artifacts/adcs_breakout_board_sim/include.json | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/artifacts/adcs_breakout_board_sim/include.json b/artifacts/adcs_breakout_board_sim/include.json index 0273495..cd065d3 100644 --- a/artifacts/adcs_breakout_board_sim/include.json +++ b/artifacts/adcs_breakout_board_sim/include.json @@ -9,7 +9,9 @@ "tasks/adcs/detumble.py:adcs/detumble.py", "tasks/adcs/point_to_earth.py:adcs/point_to_earth.py", "tasks/adcs/point_to_sun.py:adcs/point_to_sun.py", - "lib/datastores/adcs.py:datastores/adcs.py" + "lib/datastores/adcs.py:datastores/adcs.py", + "lib/tle.py:tle.py", + "lib/sgp4.py:sgp4.py" ], "unit_tests": [ "lib/pin_manager_test.py:pin_manager_test.py", From 31665a23cb370c5a6a39ecbe0d42a8ce490b9d25 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 4 May 2026 18:49:48 -0700 Subject: [PATCH 09/19] static analyser, self-referencing errors --- src/lib/datastores/adcs.py | 2 +- src/lib/tle.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 13d70fb..7ae8d3c 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -4,7 +4,7 @@ are set to `None` throughout this module. """ -import tle as tle +import tle class Datastore: """ diff --git a/src/lib/tle.py b/src/lib/tle.py index 045bb63..3e718dc 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -243,7 +243,7 @@ def from_tle_str(cls, string): return cls.from_tle_lines(*string.split('\n')[:3]) @classmethod - def sgp4_init(cls, tle: Satrec): + def sgp4_init(cls, tle): """ Creates a satrec object specifically modified to be used in sgp4. From 0b8f1210d6f8ab501cbbe0a26287e107bfb67b1f Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Mon, 4 May 2026 18:53:28 -0700 Subject: [PATCH 10/19] testing if the unit test even works right now --- src/lib/tle.py | 11 +++++++++++ unit_tests/lib/sgp4_test.py | 2 +- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/lib/tle.py b/src/lib/tle.py index 3e718dc..b10aa03 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -331,3 +331,14 @@ def error_message(self): """ if self.error == self.MOTION: return (f'mean motion {0:f} is less than zero').format(self.n) + +ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" + +if __name__ == "__main__": + sat: Satrec = Satrec.from_tle_str( + ISS_TLE + ) + + sgp4_obj = Satrec.sgp4_init(sat) + error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) + print(error, r, v) diff --git a/unit_tests/lib/sgp4_test.py b/unit_tests/lib/sgp4_test.py index f180ef5..20f6fbd 100644 --- a/unit_tests/lib/sgp4_test.py +++ b/unit_tests/lib/sgp4_test.py @@ -19,7 +19,7 @@ def propagate_to_may_3(self): sgp4_obj = Satrec.sgp4_init(sat) error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) - self.assertEqual(error, 0) # if no error, + self.assertEqual(error, 1) # if no error, tol = 3 # to 3 decimal places self.assertAlmostEqual(r[0], 4698.782358, tol) From c53440e01462f07e7abfc91c6e33fca560bb4b60 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Thu, 7 May 2026 11:52:54 -0700 Subject: [PATCH 11/19] entirely innaccurate compared to other sgp4 impls --- src/lib/{sgp4.py => sgp4s.py} | 15 +++++++---- src/lib/tle.py | 50 ++++++++++++++++++++++------------- unit_tests/lib/sgp4_test.py | 22 +++++++++------ 3 files changed, 56 insertions(+), 31 deletions(-) rename src/lib/{sgp4.py => sgp4s.py} (98%) diff --git a/src/lib/sgp4.py b/src/lib/sgp4s.py similarity index 98% rename from src/lib/sgp4.py rename to src/lib/sgp4s.py index 71dbd1f..0b24e14 100644 --- a/src/lib/sgp4.py +++ b/src/lib/sgp4s.py @@ -128,7 +128,7 @@ def sgp4_update(satrec, tsince): nm = satrec.xke / pow(am, 1.5) em = em - tempe - if em >= 1.0 or em < -0.001: + if em >= 1.0 or em < 0.0: satrec.error = satrec.ECCENTRICITY return False, False @@ -140,7 +140,7 @@ def sgp4_update(satrec, tsince): emsq = em * em temp = 1.0 - emsq - nodem = nodem % tau if nodem >= 0.0 else -(-nodem % tau) + nodem = nodem % tau argpm = argpm % tau xlm = xlm % tau mm = (xlm - argpm - nodem) % tau @@ -189,6 +189,10 @@ def sgp4_update(satrec, tsince): eo1 = eo1 + tem5 ktr = ktr + 1 + if ktr > 10: + satrec.error = satrec.MOTION # or a new error code + return False, False + # ------------- short period preliminary quantities ----------- ecose = axnl*coseo1 + aynl*sineo1 esine = axnl*sineo1 - aynl*coseo1 @@ -248,10 +252,11 @@ def sgp4_update(satrec, tsince): (mvt * uy + rvdot * vy) * vkmpersec, (mvt * uz + rvdot * vz) * vkmpersec) - if mrt < 1.0: - satrec.error = satrec.DECAY + if mrt < 1.0: + satrec.error = satrec.DECAY + return False, False - return r, v + return r, v def sgp4_init(satrec, satn, epoch, bstar, dn, ddn, ecc, argp, diff --git a/src/lib/tle.py b/src/lib/tle.py index b10aa03..2c673bc 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -4,7 +4,9 @@ Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, """ -from sgp4 import sgp4_init, sgp4_update +from sgp4s import sgp4_init, sgp4_update + +from sgp4.api import Satrec try: import ulab.numpy as np # For CircuitPython @@ -35,7 +37,9 @@ def _days2mdhms(year, days, round_to_microsecond=6): The floating point seconds are rounded to an even number of microseconds if ``round_to_microsecond`` is true. """ - second = days * 86400.0 + day_of_year, day_fraction = divmod(days, 1.0) + + second = day_fraction * 86400.0 if round_to_microsecond: second = round(second, round_to_microsecond) @@ -45,15 +49,15 @@ def _days2mdhms(year, days, round_to_microsecond=6): minute = int(minute) hour, minute = divmod(minute, 60) - day_of_year, hour = divmod(hour, 24) + hour = int(hour) is_leap = year % 400 == 0 or (year % 4 == 0 and year % 100 != 0) - month, day = _day_of_year_to_month_day(day_of_year, is_leap) + month, day = _day_of_year_to_month_day(int(day_of_year), is_leap) if month == 13: # behave like the original in case of overflow month = 12 day += 31 - return month, day, int(hour), int(minute), second + return month, day, hour, minute, second def _conv_year(s): """Interpret a two-digit year string.""" @@ -108,7 +112,7 @@ def jday(year, mon, day, hr, minute, sec): fr = (sec + minute * 60.0 + hr * 3600.0) / 86400.0 return jd, fr -class Satrec: +class Satrecs: """ Satellite record object Includes parameters, constants that are commonly used across the sgp4 logical flow @@ -267,18 +271,16 @@ def sgp4_init(cls, tle): self.argp = self.argp * deg2rad self.mo = self.mo * deg2rad - yr = self.epoch_year - _, fraction = divmod(self.epoch_day, 1.0) - self.jdsatepoch_f = round(fraction, 8) # exact number of digits in TLE + year = self.epoch_year + + mon, day, hr, minute, sec = _days2mdhms(year, self.epoch_day) + jd_full = _sgp4_jday(year, mon, day, hr, minute, sec) - # Build Julian Date - if yr < 57: - year = yr + 2000 - else: - year = yr + 1900 + # Split into two-part representation + self.jdsatepoch = np.floor(jd_full-0.5) + 0.5 # noon-to-noon JD boundary + self.jdsatepoch_f = jd_full - self.jdsatepoch - mon,day,hr,minute,sec = _days2mdhms(year, self.epoch_day) - self.jdsatepoch = _sgp4_jday(year,mon,day,hr,minute,sec) + print(self.jdsatepoch, self.jdsatepoch_f) sgp4_init(self, self.set_num, self.jdsatepoch-EPOCH0, self.bstar, self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, @@ -321,8 +323,12 @@ def sgp4_update(self, jd, fr): tsince = ((jd - self.jdsatepoch) * MIN_PER_DAY + (fr - self.jdsatepoch_f) * MIN_PER_DAY) + print(tsince) + r, v = sgp4_update(self, tsince) + print(tsince, _sgp4_jday(2026, 5, 3, 0, 0, 0)) + return self.error, r, v def error_message(self): @@ -333,12 +339,20 @@ def error_message(self): return (f'mean motion {0:f} is less than zero').format(self.n) ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" +s = "1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997" +t = "2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" + if __name__ == "__main__": - sat: Satrec = Satrec.from_tle_str( + sat: Satrecs = Satrecs.from_tle_str( ISS_TLE ) - sgp4_obj = Satrec.sgp4_init(sat) + satel = Satrec.twoline2rv(s, t) + + e, r_o, v_o = satel.sgp4(*jday(2026, 5, 3, 0, 0, 0)) + + sgp4_obj = Satrecs.sgp4_init(sat) error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) + print(e, r_o, v_o) print(error, r, v) diff --git a/unit_tests/lib/sgp4_test.py b/unit_tests/lib/sgp4_test.py index 20f6fbd..6566f10 100644 --- a/unit_tests/lib/sgp4_test.py +++ b/unit_tests/lib/sgp4_test.py @@ -8,18 +8,16 @@ ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" - -class PropagatorTest(unittest.TestCase): - - def propagate_to_may_3(self): - sat: Satrec = Satrec.from_tle_str( +SAT: Satrec = Satrec.from_tle_str( ISS_TLE ) - - sgp4_obj = Satrec.sgp4_init(sat) +sgp4_obj = Satrec.sgp4_init(SAT) + +class PropagatorTest(unittest.TestCase): + def propagation_accuracy(self): error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) - self.assertEqual(error, 1) # if no error, + self.assertEqual(error, 0) # if no error, tol = 3 # to 3 decimal places self.assertAlmostEqual(r[0], 4698.782358, tol) @@ -30,3 +28,11 @@ def propagate_to_may_3(self): self.assertAlmostEqual(v[0], 5.281325344, tol) self.assertAlmostEqual(v[1], 2.530170911, tol) self.assertAlmostEqual(v[2], -4.936649541, tol) + + def decay_orbit(self): + error, _, _ = sgp4_obj.sgp4_update(*jday(2060, 5, 3, 0, 0, 0)) + + self.assertEqual(error, 6) + +if __name__ == "__main__": + unittest.main() From 47d64488c49c4ccbd036113623948a82fa4224c8 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Fri, 15 May 2026 16:06:21 -0700 Subject: [PATCH 12/19] continued debugging --- src/lib/sgp4s.py | 36 +++++++++++++++++++++++++++++------- src/lib/tle.py | 8 +++++--- 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/src/lib/sgp4s.py b/src/lib/sgp4s.py index 0b24e14..caa2b46 100644 --- a/src/lib/sgp4s.py +++ b/src/lib/sgp4s.py @@ -36,6 +36,7 @@ def _initl(xke, j2, # only xke and j2 are used here so pass them in directly # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst x2o3 = 2.0 / 3.0; + tau = np.pi*2 # ------------- calculate auxillary epoch quantities ---------- eccsq = ecc * ecc; @@ -61,7 +62,18 @@ def _initl(xke, j2, posq = po * po rp = ao * (1.0 - ecc) - gsto = _gstime(epoch + 2433281.5) + ts70 = epoch - 7305.0; + ds70 = (ts70 + 1.0e-8) // 1.0; + tfrac = ts70 - ds70; + # find greenwich location at epoch + c1 = 1.72027916940703639e-2; + thgr70= 1.7321343856509374; + fk5r = 5.07551419432269442e-15; + c1p2p = c1 + tau; + gsto = (thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r) % tau + if gsto < 0.0: + gsto = gsto + tau + # gsto = _gstime(epoch + 2433281.5) return (no, ao, con41, con42, cosio, @@ -189,10 +201,6 @@ def sgp4_update(satrec, tsince): eo1 = eo1 + tem5 ktr = ktr + 1 - if ktr > 10: - satrec.error = satrec.MOTION # or a new error code - return False, False - # ------------- short period preliminary quantities ----------- ecose = axnl*coseo1 + aynl*sineo1 esine = axnl*sineo1 - aynl*coseo1 @@ -252,11 +260,11 @@ def sgp4_update(satrec, tsince): (mvt * uy + rvdot * vy) * vkmpersec, (mvt * uz + rvdot * vz) * vkmpersec) - if mrt < 1.0: + if mrt < 1.0: satrec.error = satrec.DECAY return False, False - return r, v + return r, v def sgp4_init(satrec, satn, epoch, bstar, dn, ddn, ecc, argp, @@ -342,6 +350,20 @@ def sgp4_init(satrec, satn, epoch, satrec.isimp = 1 sfour = ss qzms24 = qzms2t + perige = (rp - 1.0) * satrec.radiusearthkm; + + # - for perigees below 156 km, s and qoms2t are altered - + if perige < 156.0: + + sfour = perige - 78.0; + if perige < 98.0: + sfour = 20.0; + # sgp4fix use multiply for speed instead of pow + qzms24temp = (120.0 - sfour) / satrec.radiusearthkm; + qzms24 = qzms24temp * qzms24temp * qzms24temp * qzms24temp; + sfour = sfour / satrec.radiusearthkm + 1.0; + print(perige) + pinvsq = 1.0 / posq tsi = 1.0 / (ao - sfour) diff --git a/src/lib/tle.py b/src/lib/tle.py index 2c673bc..3dc2178 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -249,7 +249,7 @@ def from_tle_str(cls, string): @classmethod def sgp4_init(cls, tle): """ - Creates a satrec object specifically modified to be used in sgp4. + Creates a satrec object specifically modified from TLE to be used in sgp4. Changes units, activates certain new parameters, etc. """ @@ -282,9 +282,11 @@ def sgp4_init(cls, tle): print(self.jdsatepoch, self.jdsatepoch_f) - sgp4_init(self, self.set_num, self.jdsatepoch-EPOCH0, self.bstar, + sgp4_init(self, self.set_num, jd_full-EPOCH0, self.bstar, self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, self.raan) + + print(self.jdsatepoch, self.jdsatepoch_f) return self @@ -352,7 +354,7 @@ def error_message(self): e, r_o, v_o = satel.sgp4(*jday(2026, 5, 3, 0, 0, 0)) - sgp4_obj = Satrecs.sgp4_init(sat) + sgp4_obj : Satrecs= Satrecs.sgp4_init(sat) error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) print(e, r_o, v_o) print(error, r, v) From 3311e057f060a0beea15f75ac8ad97cccf1009dd Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:29:15 +0800 Subject: [PATCH 13/19] Stripped TLE without ADCS --- src/lib/sgp4s.py | 440 ----------------------------------------------- src/lib/tle.py | 188 +------------------- 2 files changed, 5 insertions(+), 623 deletions(-) delete mode 100644 src/lib/sgp4s.py diff --git a/src/lib/sgp4s.py b/src/lib/sgp4s.py deleted file mode 100644 index caa2b46..0000000 --- a/src/lib/sgp4s.py +++ /dev/null @@ -1,440 +0,0 @@ -""" -Simplified General Perturbations Model 4 Implementation for Orbit Propagation - -As outlined in https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf -""" - -# Error codes -ECCENTRICITY = 1 # eccentricity is not within 0-1 -MOTION = 2 # error in propagating mean motion -SEMIRECT = 4 # apoapsis, periapsis characteristics error -DECAY = 6 # orbit has decayed - -try: - import ulab.numpy as np # For CircuitPython -except ImportError: - import numpy as np # For GitHub Actions / PC testing - -def _gstime(jdut1): - deg2rad = np.pi / 180.0 - tau = np.pi*2 - tut1 = (jdut1 - 2451545.0) / 36525.0 - temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \ - (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841 # sec - temp = (temp * deg2rad / 240.0) % tau # 360/86400 = 1/240, to deg, to rad - - # ------------------------ check quadrants --------------------- - if temp < 0.0: - temp += tau - - return temp - -def _initl(xke, j2, - ecc, epoch, inclo, no): - # ----------------------- earth constants ---------------------- - # sgp4fix identify constants and allow alternate values - # only xke and j2 are used here so pass them in directly - # tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst - x2o3 = 2.0 / 3.0; - tau = np.pi*2 - - # ------------- calculate auxillary epoch quantities ---------- - eccsq = ecc * ecc; - omeosq = 1.0 - eccsq; - rteosq = np.sqrt(omeosq); - cosio = np.cos(inclo); - cosio2 = cosio * cosio; - - # ------------------ un-kozai the mean motion ----------------- - ak = pow(xke / no, x2o3); - d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); - del_ = d1 / (ak * ak); - adel = ak * (1.0 - del_ * del_ - del_ * - (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); - del_ = d1/(adel * adel) - no = no / (1.0 + del_) - - ao = pow(xke / no, x2o3) - sinio = np.sin(inclo) - po = ao * omeosq - con42 = 1.0 - 5.0 * cosio2 - con41 = -con42-cosio2-cosio2 - posq = po * po - rp = ao * (1.0 - ecc) - - ts70 = epoch - 7305.0; - ds70 = (ts70 + 1.0e-8) // 1.0; - tfrac = ts70 - ds70; - # find greenwich location at epoch - c1 = 1.72027916940703639e-2; - thgr70= 1.7321343856509374; - fk5r = 5.07551419432269442e-15; - c1p2p = c1 + tau; - gsto = (thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r) % tau - if gsto < 0.0: - gsto = gsto + tau - # gsto = _gstime(epoch + 2433281.5) - - return (no, - ao, con41, con42, cosio, - cosio2, omeosq, posq, - rp, rteosq,sinio , gsto - ) - -def sgp4_update(satrec, tsince): - """ - Transforms to cartesian (r, v) from a Satrec (TLE 7+) input, then performs matrix operations - to propagate the system - - Current implementation uses an analytical approach to compute the Jacobian A matrix - Possible alternative is numerical stepping (for lower accuracy) of associated partials - """ - # -- Set mathematical constants - x2o3 = 2.0 / 3.0 - tau = 2.0 * np.pi - vkmpersec = satrec.radiusearthkm * satrec.xke/60.0; - - satrec.t = tsince - - # -- Update for secular gravity and atmospheric drag - xmdf = satrec.mo + satrec.mdot * satrec.t - argpdf = satrec.argp + satrec.argpdot * satrec.t - nodedf = satrec.raan + satrec.nodedot * satrec.t - argpm = argpdf - mm = xmdf - t2 = satrec.t * satrec.t - nodem = nodedf + satrec.nodecf * t2 - tempa = 1.0 - satrec.cc1 * satrec.t - tempe = satrec.bstar * satrec.cc4 * satrec.t - templ = satrec.t2cof * t2 - - if satrec.isimp != 1: - - delomg = satrec.omgcof * satrec.t - # sgp4fix use mutliply for speed instead of pow - delmtemp = 1.0 + satrec.eta * np.cos(xmdf) - delm = satrec.xmcof * \ - (delmtemp * delmtemp * delmtemp - - satrec.delmo) - temp = delomg + delm - mm = xmdf + temp - argpm = argpdf - temp - t3 = t2 * satrec.t - t4 = t3 * satrec.t - tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \ - satrec.d4 * t4 - tempe = tempe + satrec.bstar * satrec.cc5 * (np.sin(mm) - - satrec.sinmao) - templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + - satrec.t * satrec.t5cof) - - nm = satrec.n - em = satrec.ecc - inclm = satrec.inclo - - if nm <= 0.0: - satrec.error = satrec.MOTION - return False, False - - am = pow((satrec.xke / nm),x2o3) * tempa * tempa - nm = satrec.xke / pow(am, 1.5) - em = em - tempe - - if em >= 1.0 or em < 0.0: - satrec.error = satrec.ECCENTRICITY - - return False, False - - if em < 1.0e-6: - em = 1.0e-6 - mm = mm + satrec.n * templ - xlm = mm + argpm + nodem - emsq = em * em - temp = 1.0 - emsq - - nodem = nodem % tau - argpm = argpm % tau - xlm = xlm % tau - mm = (xlm - argpm - nodem) % tau - - satrec.am = am - satrec.em = em - satrec.im = inclm - satrec.Om = nodem - satrec.om = argpm - satrec.mm = mm - satrec.nm = nm - - # ----------------- compute extra mean quantities ------------- - sinim = np.sin(inclm) - cosim = np.cos(inclm) - - # -------------------- add lunar-solar periodics -------------- - ep = em - xincp = inclm - argpp = argpm - nodep = nodem - mp = mm - sinip = sinim - cosip = cosim - - # -------------------- long period periodics ------------------ - axnl = ep * np.cos(argpp) - temp = 1.0 / (am * (1.0 - ep * ep)) - aynl = ep* np.sin(argpp) + temp * satrec.aycof - xl = mp + argpp + nodep + temp * satrec.xlcof * axnl - - # --------------------- solve kepler's equation --------------- - u = (xl - nodep) % tau - eo1 = u - tem5 = 9999.9 - ktr = 1 - - while np.fabs(tem5) >= 1.0e-12 and ktr <= 10: - - sineo1 = np.sin(eo1) - coseo1 = np.cos(eo1) - tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl - tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5 - if np.fabs(tem5) >= 0.95: - tem5 = 0.95 if tem5 > 0.0 else -0.95 - eo1 = eo1 + tem5 - ktr = ktr + 1 - - # ------------- short period preliminary quantities ----------- - ecose = axnl*coseo1 + aynl*sineo1 - esine = axnl*sineo1 - aynl*coseo1 - el2 = axnl*axnl + aynl*aynl - pl = am*(1.0-el2) - if pl < 0.0: - satrec.error = satrec.SEMIRECT - - return False, False - - else: - - rl = am * (1.0 - ecose) - rdotl = np.sqrt(am) * esine/rl - rvdotl = np.sqrt(pl) / rl - betal = np.sqrt(1.0 - el2) - temp = esine / (1.0 + betal) - sinu = am / rl * (sineo1 - aynl - axnl * temp) - cosu = am / rl * (coseo1 - axnl + aynl * temp) - su = np.atan2(sinu, cosu) - sin2u = (cosu + cosu) * sinu - cos2u = 1.0 - 2.0 * sinu * sinu - temp = 1.0 / pl - temp1 = 0.5 * satrec.j2 * temp - temp2 = temp1 * temp - - # -------------- update for short period periodics ------------ - mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \ - 0.5 * temp1 * satrec.x1mth2 * cos2u - su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u - xnode = nodep + 1.5 * temp2 * cosip * sin2u - xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u - mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / satrec.xke - rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + - 1.5 * satrec.con41) / satrec.xke - - # --------------------- orientation vectors ------------------- - sinsu = np.sin(su) - cossu = np.cos(su) - snod = np.sin(xnode) - cnod = np.cos(xnode) - sini = np.sin(xinc) - cosi = np.cos(xinc) - xmx = -snod * cosi - xmy = cnod * cosi - ux = xmx * sinsu + cnod * cossu - uy = xmy * sinsu + snod * cossu - uz = sini * sinsu - vx = xmx * cossu - cnod * sinsu - vy = xmy * cossu - snod * sinsu - vz = sini * cossu - - # --------- position and velocity (in km and km/sec) ---------- - _mr = mrt * satrec.radiusearthkm - r = (_mr * ux, _mr * uy, _mr * uz) - v = ((mvt * ux + rvdot * vx) * vkmpersec, - (mvt * uy + rvdot * vy) * vkmpersec, - (mvt * uz + rvdot * vz) * vkmpersec) - - if mrt < 1.0: - satrec.error = satrec.DECAY - return False, False - - return r, v - -def sgp4_init(satrec, satn, epoch, - bstar, dn, ddn, ecc, argp, - inclo, mo, - raan - ): - """ - (Math Heavy) initialisation of sgp4 which instantiates all important - variables in the satrec object - """ - - temp4 = 1.5e-12 - - # Near Earth Variables - satrec.isimp = 0; satrec.aycof = 0.0 - satrec.con41 = 0.0; satrec.cc1 = 0.0; satrec.cc4 = 0.0 - satrec.cc5 = 0.0; satrec.d2 = 0.0; satrec.d3 = 0.0 - satrec.d4 = 0.0; satrec.delmo = 0.0; satrec.eta = 0.0 - satrec.argpdot = 0.0; satrec.omgcof = 0.0; satrec.sinmao = 0.0 - satrec.t = 0.0; satrec.t2cof = 0.0; satrec.t3cof = 0.0 - satrec.t4cof = 0.0; satrec.t5cof = 0.0; satrec.x1mth2 = 0.0 - satrec.x7thm1 = 0.0; satrec.mdot = 0.0; satrec.nodedot = 0.0 - satrec.xlcof = 0.0; satrec.xmcof = 0.0; satrec.nodecf = 0.0 - - # Earth Constants - satrec.mu = 398600.5 # in km3 / s2 - satrec.radiusearthkm = 6378.137 # km - satrec.xke = 0.07436685317 - satrec.tumin = 13.44685108 - satrec.j2 = 0.00108262998905 - satrec.j3 = -0.00000253215306 - satrec.j4 = -0.00000161098761 - satrec.j3oj2 = -0.002338890559 - - ss = 1.012229276 - qzms2ttemp = 0.00658499496 - qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp - x2o3 = 2.0 / 3.0 - - # -- initialisation markers - satrec.t = 0.0 - - # -- - satrec.norad = satn - satrec.classification = 'U' - - # -- - satrec.bstar = bstar - satrec.dn = dn - satrec.ddn = ddn - satrec.ecc = ecc - satrec.argp = argp - satrec.inclo = inclo - satrec.mo = mo - satrec.raan = raan - - # single averaged mean elements - satrec.am = 0.0 - satrec.em = 0.0 - satrec.im = 0.0 - satrec.Om = 0.0 - satrec.mm = 0.0 - satrec.nm = 0.0 - - satrec.error = 0 - - # -- - ( - satrec.n, - ao, satrec.con41, con42, cosio, - cosio2, omeosq, posq, - rp, rteosq,sinio , satrec.gsto, - ) = _initl( - satrec.xke, satrec.j2, satrec.ecc, epoch, satrec.inclo, satrec.n - ) - satrec.a = pow( satrec.n*satrec.tumin , (-2.0/3.0) ) - satrec.alta = satrec.a*(1.0 + satrec.ecc) - 1.0 - satrec.altp = satrec.a*(1.0 - satrec.ecc) - 1.0 - - if omeosq >= 0.0 or satrec.n >= 0.0: - satrec.isimp = 0 - if rp < 220.0 / satrec.radiusearthkm + 1.0: - satrec.isimp = 1 - sfour = ss - qzms24 = qzms2t - perige = (rp - 1.0) * satrec.radiusearthkm; - - # - for perigees below 156 km, s and qoms2t are altered - - if perige < 156.0: - - sfour = perige - 78.0; - if perige < 98.0: - sfour = 20.0; - # sgp4fix use multiply for speed instead of pow - qzms24temp = (120.0 - sfour) / satrec.radiusearthkm; - qzms24 = qzms24temp * qzms24temp * qzms24temp * qzms24temp; - sfour = sfour / satrec.radiusearthkm + 1.0; - print(perige) - - pinvsq = 1.0 / posq - - tsi = 1.0 / (ao - sfour) - satrec.eta = ao * satrec.ecc * tsi - etasq = satrec.eta * satrec.eta - eeta = satrec.ecc * satrec.eta - psisq = np.fabs(1.0 - etasq) - coef = qzms24 * pow(tsi, 4.0); - coef1 = coef / pow(psisq, 3.5) - cc2 = coef1 * satrec.n * (ao * (1.0 + 1.5 * etasq + eeta * - (4.0 + etasq)) + 0.375 * satrec.j2 * tsi / psisq * satrec.con41 * - (8.0 + 3.0 * etasq * (8.0 + etasq))) - satrec.cc1 = satrec.bstar * cc2 - cc3 = 0.0 - if satrec.ecc > 1.0e-4: - cc3 = -2.0 * coef * tsi * satrec.j3oj2 * satrec.n * sinio / satrec.ecc - satrec.x1mth2 = 1.0 - cosio2 - satrec.cc4 = 2.0* satrec.n * coef1 * ao * omeosq * \ - (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecc * - (0.5 + 2.0 * etasq) - satrec.j2 * tsi / (ao * psisq) * - (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * - (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 * - (2.0 * etasq - eeta * (1.0 + etasq)) * np.cos(2.0 * satrec.argp))) - satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * - (etasq + eeta) + eeta * etasq) - cosio4 = cosio2 * cosio2 - temp1 = 1.5 * satrec.j2 * pinvsq * satrec.n - temp2 = 0.5 * temp1 * satrec.j2 * pinvsq - temp3 = -0.46875 * satrec.j4 * pinvsq * pinvsq * satrec.n - satrec.mdot = satrec.n + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 * \ - temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4) - satrec.argpdot = (-0.5 * temp1 * con42 + 0.0625 * temp2 * - (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + - temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4)) - xhdot1 = -temp1 * cosio - satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + - 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio - satrec.omgcof = satrec.bstar * cc3 * np.cos(satrec.argp) - satrec.xmcof = 0.0 - if satrec.ecc > 1.0e-4: - satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta - satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1 - satrec.t2cof = 1.5 * satrec.cc1 - - if np.fabs(cosio+1.0) > 1.5e-12: - satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio) - else: - satrec.xlcof = -0.25 * satrec.j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4 - satrec.aycof = -0.5 * satrec.j3oj2 * sinio - - delmotemp = 1.0 + satrec.eta * np.cos(satrec.mo) - satrec.delmo = delmotemp * delmotemp * delmotemp - satrec.sinmao = np.sin(satrec.mo) - satrec.x7thm1 = 7.0 * cosio2 - 1.0 - - if satrec.isimp != 1: - cc1sq = satrec.cc1 * satrec.cc1 - satrec.d2 = 4.0 * ao * tsi * cc1sq - temp = satrec.d2 * tsi * satrec.cc1 / 3.0 - satrec.d3 = (17.0 * ao + sfour) * temp - satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * \ - satrec.cc1 - satrec.t3cof = satrec.d2 + 2.0 * cc1sq - satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 * - (12.0 * satrec.d2 + 10.0 * cc1sq)) - satrec.t5cof = 0.2 * (3.0 * satrec.d4 + - 12.0 * satrec.cc1 * satrec.d3 + - 6.0 * satrec.d2 * satrec.d2 + - 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq)) - - # propagate to 0 - sgp4_update(satrec, 0) - - return True diff --git a/src/lib/tle.py b/src/lib/tle.py index 3dc2178..2c05ae3 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -1,64 +1,9 @@ """ -Functions and Variables used by ADCS to update and use TLE (two-line element data) +Functions and Variables to update and use TLE (two-line element data) -Adapted from the TLE-tools library by @FedericoStra on GitHub for ulab, +Adapted from the TLE-tools library by @FedericoStra on GitHub """ -from sgp4s import sgp4_init, sgp4_update - -from sgp4.api import Satrec - -try: - import ulab.numpy as np # For CircuitPython -except ImportError: - import numpy as np # For GitHub Actions / PC testing - -MIN_PER_DAY = 1440 -EPOCH0 = 2433281.5 # jan 0 1950 - -def _day_of_year_to_month_day(day_of_year, is_leap): - """Core logic for turning days into months, for easy testing.""" - february_bump = (2 - is_leap) * (day_of_year >= 60 + is_leap) - august = day_of_year >= 215 - month, day = divmod(2 * (day_of_year - 1 + 30 * august + february_bump), 61) - month += 1 - august - day //= 2 - day += 1 - return month, day - -def _days2mdhms(year, days, round_to_microsecond=6): - """Convert a float point number of days into the year into date and time. - - >>> days2mdhms(2000, 32.0) # February 1 - (2, 1, 0, 0, 0.0) - >>> days2mdhms(2000, 366.0) # December 31, since 2000 was a leap year - (12, 31, 0, 0, 0.0) - - The floating point seconds are rounded to an even number of - microseconds if ``round_to_microsecond`` is true. - """ - day_of_year, day_fraction = divmod(days, 1.0) - - second = day_fraction * 86400.0 - if round_to_microsecond: - second = round(second, round_to_microsecond) - - minute, second = divmod(second, 60.0) - if round_to_microsecond: - second = round(second, round_to_microsecond) - - minute = int(minute) - hour, minute = divmod(minute, 60) - hour = int(hour) - - is_leap = year % 400 == 0 or (year % 4 == 0 and year % 100 != 0) - month, day = _day_of_year_to_month_day(int(day_of_year), is_leap) - if month == 13: # behave like the original in case of overflow - month = 12 - day += 31 - - return month, day, hour, minute, second - def _conv_year(s): """Interpret a two-digit year string.""" if isinstance(s, int): @@ -86,36 +31,9 @@ def _parse_float(s): """ return float(s[0] + '.' + s[1:6] + 'e' + s[6:8]) -def _sgp4_jday(year, mon, day, hr, minute, sec): - """ - Converts jdsatepoch into a compatible number for sgp4 - """ - - return (367.0 * year - - 7.0 * (year + ((mon + 9.0) // 12.0)) * 0.25 // 1.0 + - 275.0 * mon // 9.0 + - day + 1721013.5 + - ((sec / 60.0 + minute) / 60.0 + hr) / 24.0 - ) - -def jday(year, mon, day, hr, minute, sec): - """ - From a date, return a Julian date in its date + fractional form - - Used to build the jd, fr parameters in the function call of sgp4_update - """ - jd = (367.0 * year - - 7 * (year + ((mon + 9) // 12.0)) * 0.25 // 1.0 - + 275 * mon / 9.0 // 1.0 - + day - + 1721013.5) - fr = (sec + minute * 60.0 + hr * 3600.0) / 86400.0 - return jd, fr - -class Satrecs: +class Satrec: """ Satellite record object - Includes parameters, constants that are commonly used across the sgp4 logical flow In this implementation, built from TLE data @@ -161,12 +79,6 @@ class Satrecs: Revolution number. """ - # Error codes - ECCENTRICITY = 1 # eccentricity is not within 0-1 - MOTION = 2 # error in propagating mean motion - SEMIRECT = 4 # apoapsis, periapsis characteristics error - DECAY = 6 # orbit has decayed - def __init__(self, name:str, # ID parameters, Line 1 norad:str, classification:str, int_desig:str, @@ -202,10 +114,6 @@ def __init__(self, name:str, self.tle_str = tle_str - # for sgp - self.jdsatepoch = 0 - self.jdsatepoch_f = 0 - @classmethod def from_tle_lines(cls, name, line1, line2): """Parse a TLE from its constituent lines. @@ -238,7 +146,7 @@ def from_tle_lines(cls, name, line1, line2): def from_tle_file(cls, filename): """Load TLE from a file.""" if isinstance(filename, str): - with open(filename) as fp: + with open(filename, encoding="utf-8") as fp: return cls.from_tle_lines(*fp.readlines[:2]) @classmethod @@ -246,54 +154,10 @@ def from_tle_str(cls, string): """Load TLE from a string.""" return cls.from_tle_lines(*string.split('\n')[:3]) - @classmethod - def sgp4_init(cls, tle): - """ - Creates a satrec object specifically modified from TLE to be used in sgp4. - - Changes units, activates certain new parameters, etc. - """ - - self = tle - - # constants for unit change - deg2rad = np.pi / 180.0 # 0.0174532925199433 - xpdotp = 1440.0 / (2.0 *np.pi) # 229.1831180523293 - - # ---- convert to sgp4 units ---- - self.n = self.n / xpdotp - self.dn = self.dn / (xpdotp*1440.0) - self.ddn= self.ddn / (xpdotp*1440.0*1440) - - # ---- find standard orbital elements ---- - self.inc = self.inc * deg2rad - self.raan = self.raan * deg2rad - self.argp = self.argp * deg2rad - self.mo = self.mo * deg2rad - - year = self.epoch_year - - mon, day, hr, minute, sec = _days2mdhms(year, self.epoch_day) - jd_full = _sgp4_jday(year, mon, day, hr, minute, sec) - - # Split into two-part representation - self.jdsatepoch = np.floor(jd_full-0.5) + 0.5 # noon-to-noon JD boundary - self.jdsatepoch_f = jd_full - self.jdsatepoch - - print(self.jdsatepoch, self.jdsatepoch_f) - - sgp4_init(self, self.set_num, jd_full-EPOCH0, self.bstar, - self.dn, self.ddn, self.ecc, self.argp, self.inc, self.n, - self.raan) - - print(self.jdsatepoch, self.jdsatepoch_f) - - return self - def to_array(self): """ Return 2D array of TLE values - + Indexed as [line, col] @@ -316,45 +180,3 @@ def to_array(self): # line 2 orbital params [self.inc, self.raan, self.ecc, self.argp, self.mo, self.n, self.rev_num] ] - - def sgp4_update(self, jd, fr): - """ - For a julian date (jd) and its fractional representation (fr), - propagate Satrec using sgp4 - """ - - tsince = ((jd - self.jdsatepoch) * MIN_PER_DAY + - (fr - self.jdsatepoch_f) * MIN_PER_DAY) - print(tsince) - - r, v = sgp4_update(self, tsince) - - print(tsince, _sgp4_jday(2026, 5, 3, 0, 0, 0)) - - return self.error, r, v - - def error_message(self): - """ - Return error message from current self.error (when polled) - """ - if self.error == self.MOTION: - return (f'mean motion {0:f} is less than zero').format(self.n) - -ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" -s = "1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997" -t = "2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" - - -if __name__ == "__main__": - sat: Satrecs = Satrecs.from_tle_str( - ISS_TLE - ) - - satel = Satrec.twoline2rv(s, t) - - e, r_o, v_o = satel.sgp4(*jday(2026, 5, 3, 0, 0, 0)) - - sgp4_obj : Satrecs= Satrecs.sgp4_init(sat) - error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) - print(e, r_o, v_o) - print(error, r, v) From 30d2b5c1614a0fc3deb23a9c31cebd2b845453ec Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:32:12 +0800 Subject: [PATCH 14/19] removed unit testing --- .../adcs_breakout_board_sim/include.json | 6 +-- unit_tests/lib/sgp4_test.py | 38 ------------------- 2 files changed, 2 insertions(+), 42 deletions(-) delete mode 100644 unit_tests/lib/sgp4_test.py diff --git a/artifacts/adcs_breakout_board_sim/include.json b/artifacts/adcs_breakout_board_sim/include.json index cd065d3..6d84119 100644 --- a/artifacts/adcs_breakout_board_sim/include.json +++ b/artifacts/adcs_breakout_board_sim/include.json @@ -10,14 +10,12 @@ "tasks/adcs/point_to_earth.py:adcs/point_to_earth.py", "tasks/adcs/point_to_sun.py:adcs/point_to_sun.py", "lib/datastores/adcs.py:datastores/adcs.py", - "lib/tle.py:tle.py", - "lib/sgp4.py:sgp4.py" + "lib/tle.py:tle.py" ], "unit_tests": [ "lib/pin_manager_test.py:pin_manager_test.py", "lib/custom_module_mocking.py:custom_module_mocking.py", - "lib/quaternion_test.py:quaternion_test.py", - "lib/sgp4_test.py:sgp4_test.py" + "lib/quaternion_test.py:quaternion_test.py" ], "submodules":[ "Adafruit_CircuitPython_Ticks/adafruit_ticks.py:adafruit_ticks.py", diff --git a/unit_tests/lib/sgp4_test.py b/unit_tests/lib/sgp4_test.py deleted file mode 100644 index 6566f10..0000000 --- a/unit_tests/lib/sgp4_test.py +++ /dev/null @@ -1,38 +0,0 @@ -import unittest -from tle import Satrec, jday - -try: - import ulab.numpy as np # For CircuitPython -except ImportError: - import numpy as np # For GitHub Actions / PC testing - - -ISS_TLE = "ISS (ZARYA)\n1 25544U 98067A 26121.81277072 .00006771 00000-0 13067-3 0 9997\n2 25544 51.6311 169.6452 0007227 12.7206 347.3964 15.49051775564564" -SAT: Satrec = Satrec.from_tle_str( - ISS_TLE - ) -sgp4_obj = Satrec.sgp4_init(SAT) - -class PropagatorTest(unittest.TestCase): - def propagation_accuracy(self): - error, r, v = sgp4_obj.sgp4_update(*jday(2026, 5, 3, 0, 0, 0)) - - self.assertEqual(error, 0) # if no error, - - tol = 3 # to 3 decimal places - self.assertAlmostEqual(r[0], 4698.782358, tol) - self.assertAlmostEqual(r[1], -3867.014434, tol) - self.assertAlmostEqual(r[2], 3028.549126, tol) - - tol = 6 # to 6 decimal places, these need to be more accurate - self.assertAlmostEqual(v[0], 5.281325344, tol) - self.assertAlmostEqual(v[1], 2.530170911, tol) - self.assertAlmostEqual(v[2], -4.936649541, tol) - - def decay_orbit(self): - error, _, _ = sgp4_obj.sgp4_update(*jday(2060, 5, 3, 0, 0, 0)) - - self.assertEqual(error, 6) - -if __name__ == "__main__": - unittest.main() From fb4aa65d544695b57c3efc53cd0ea6b64aea5e23 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:36:14 +0800 Subject: [PATCH 15/19] datastore errors --- src/lib/datastores/adcs.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 7ae8d3c..4c45ab6 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -11,6 +11,11 @@ class Datastore: Datastore class for adcs processes. Holds time, sensor, and attitude data to be used system-wide """ + # TLE String + TLE = """ISS (ZARYA)\n +1 25544U 98067A 26166.51237796 .00007685 00000-0 14626-3 0 9999\n +2 25544 51.6337 308.3821 0004850 189.0196 171.0706 15.49243792571497""" + # Action types DETUMBLE = 0 POINT_TO_SUN = 1 @@ -35,7 +40,7 @@ def __init__(self): None # Quaternion representing attitude from body frame to inertial frame ) self.mode = self.DETUMBLE - self.tle: tle.TLE = tle.TLE() + # self.satrecs: tle.Satrec = tle.Satrec.from_tle_str(TLE) class AdcsTime: """ From 4039b065bb191f01f2afebc983ed5bfe0ce0dd97 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:37:27 +0800 Subject: [PATCH 16/19] unused import --- src/lib/datastores/adcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib/datastores/adcs.py b/src/lib/datastores/adcs.py index 4c45ab6..f85a11c 100644 --- a/src/lib/datastores/adcs.py +++ b/src/lib/datastores/adcs.py @@ -4,7 +4,7 @@ are set to `None` throughout this module. """ -import tle +# import tle class Datastore: """ From c1125bf640718fd6a77b05a30072114af6aba85c Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:44:38 +0800 Subject: [PATCH 17/19] Made default instance blank, from_tle_lines as new instantiation --- src/lib/tle.py | 98 +++++++++++++++++++++++--------------------------- 1 file changed, 44 insertions(+), 54 deletions(-) diff --git a/src/lib/tle.py b/src/lib/tle.py index 2c05ae3..9b408b2 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -79,40 +79,31 @@ class Satrec: Revolution number. """ - def __init__(self, name:str, - # ID parameters, Line 1 - norad:str, classification:str, int_desig:str, - # time (derivative) parameters, line 1 - epoch_year:int, epoch_day:float, dn:float, ddn:float, bstar:float, - ephtype: str, set_num:int, - # keplerian parameters, line 2 - inc:float, raan:float, ecc:float, argp:float, mo:float, n:float, rev_num:int, - # for the purposes of keeping the tle around as future-proofing - tle_str:str ): - - self.name = str.strip(name) - - self.norad = str.strip(norad) - self.classification = classification - self.int_desig = str.strip(int_desig) - - self.epoch_year = _conv_year(epoch_year) - self.epoch_day = epoch_day - self.dn = dn - self.ddn = ddn - self.bstar = bstar - self.ephtype = ephtype - self.set_num = int(set_num) - - self.inc = inc - self.raan = raan - self.ecc = ecc - self.argp = argp - self.mo = mo - self.n = n # mean motion - self.rev_num = int(rev_num) - - self.tle_str = tle_str + def __init__(self): + + self.name = "" + + self.norad = "" + self.classification = "" + self.int_desig = "" + + self.epoch_year = 0 + self.epoch_day = 0.0 + self.dn = 0.0 + self.ddn = 0.0 + self.bstar = 0.0 + self.ephtype = "" + self.set_num = 0 + + self.inc = 0.0 + self.raan = 0.0 + self.ecc = 0.0 + self.argp = 0.0 + self.mo = 0.0 + self.n = 0.0 # mean motion + self.rev_num = 0 + + self.tle_str = "" @classmethod def from_tle_lines(cls, name, line1, line2): @@ -121,26 +112,25 @@ def from_tle_lines(cls, name, line1, line2): All the attributes parsed from the TLE are expressed in the same units that are used in the TLE format. """ - return cls( - name=name, - norad=line1[2:7], - classification=line1[7] or 'U', - int_desig=line1[9:17], - epoch_year=line1[18:20], - epoch_day=float(line1[20:32]), - dn=float(line1[33:43]), - ddn=_parse_float(line1[44:52]), - bstar=_parse_float(line1[53:61]), - ephtype = line1[62], - set_num=line1[64:68], - inc=float(line2[8:16]), - raan=float(line2[17:25]), - ecc=_parse_decimal(line2[26:33]), - argp=float(line2[34:42]), - mo=float(line2[43:51]), - n=float(line2[52:63]), - rev_num=line2[63:68], - tle_str=name+line1+line2) + cls.name=name + cls.norad=line1[2:7] + cls.classification=line1[7] or 'U' + cls.int_desig=line1[9:17] + cls.epoch_year=_conv_year(line1[18:20]) + cls.epoch_day=float(line1[20:32]) + cls.dn=float(line1[33:43]) + cls.bstar=_parse_float(line1[53:61]) + cls.ddn=_parse_float(line1[44:52]) + cls.ephtype = line1[62] + cls.set_num=line1[64:68] + cls.inc=float(line2[8:16]) + cls.raan=float(line2[17:25]) + cls.ecc=_parse_decimal(line2[26:33]) + cls.argp=float(line2[34:42]) + cls.mo=float(line2[43:51]) + cls.n=float(line2[52:63]) + cls.rev_num=line2[63:68] + cls.tle_str=name+line1+line2 @classmethod def from_tle_file(cls, filename): From e0116bf7651c2adfacc59b84b8f7f157b14bb625 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Tue, 16 Jun 2026 09:46:15 +0800 Subject: [PATCH 18/19] static analyser --- src/lib/tle.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lib/tle.py b/src/lib/tle.py index 9b408b2..1f9b30b 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -139,6 +139,8 @@ def from_tle_file(cls, filename): with open(filename, encoding="utf-8") as fp: return cls.from_tle_lines(*fp.readlines[:2]) + return None + @classmethod def from_tle_str(cls, string): """Load TLE from a string.""" From 99cece04810cb6d30bba7868dc4bf901304de289 Mon Sep 17 00:00:00 2001 From: TinFoiKa Date: Sat, 27 Jun 2026 13:04:11 +0800 Subject: [PATCH 19/19] initialiser was implicit - made explicit --- src/lib/tle.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lib/tle.py b/src/lib/tle.py index 1f9b30b..b0bb76f 100644 --- a/src/lib/tle.py +++ b/src/lib/tle.py @@ -112,6 +112,8 @@ def from_tle_lines(cls, name, line1, line2): All the attributes parsed from the TLE are expressed in the same units that are used in the TLE format. """ + + cls() cls.name=name cls.norad=line1[2:7] cls.classification=line1[7] or 'U'