1# -*- coding: utf-8 -*- 2"""An interface between Skyfield and the Python ``sgp4`` library.""" 3 4from numpy import ( 5 array, concatenate, identity, multiply, ones_like, repeat, 6) 7from sgp4.api import SGP4_ERRORS, Satrec 8 9from .constants import AU_KM, DAY_S, T0, tau 10from .functions import _T, mxm, mxv, rot_x, rot_y, rot_z 11from .searchlib import _find_discrete, find_maxima 12from .timelib import compute_calendar_date 13from .vectorlib import VectorFunction 14 15_identity = identity(3) 16 17class EarthSatellite(VectorFunction): 18 """An Earth satellite loaded from a TLE file and propagated with SGP4. 19 20 An earth satellite object is a Skyfield vector function, so you can 21 either call its ``at()`` method to generate its position in the sky 22 or else use addition and subtraction to combine it with other 23 vectors. 24 25 Satellite parameters are generally only accurate for a week or two 26 around the *epoch* of the parameters, the date for which they were 27 generated, which is available as an attribute: 28 29 ``epoch`` 30 A Skyfield :class:`~skyfield.timelib.Time` giving the exact 31 epoch moment for these satellite orbit parameters. 32 ``name`` 33 Satellite name 34 35 When building a satellite, use the arguments ``line1`` and ``line2`` 36 to provide the two data lines from a TLE file as separate strings. 37 Optional ``name`` lets you give a name to the satellite, accessible 38 later through the ``name`` attribute. ``ts`` is a 39 :class:`~skyfield.timelib.Timescale` object, used to generate the 40 ``epoch`` value; if it is not provided, the satellite will use a 41 built in ``Timescale`` object. 42 43 If you are interested in the catalog entry details, the SGP4 model 44 parameters for a particular satellite can be accessed through its 45 ``model`` attribute: 46 47 ``model.satnum`` 48 The unique satellite NORAD catalog number given in the TLE file. 49 ``model.classification`` 50 Satellite classification, or else ``'U'`` for “Unknown” 51 ``model.intldesg`` 52 International designator 53 ``model.epochyr`` 54 Full four-digit year of this element set's epoch moment. 55 ``model.epochdays`` 56 Fractional days into the year of the epoch moment. 57 ``model.jdsatepoch`` 58 Julian date of the epoch (computed from ``epochyr`` and ``epochdays``). 59 ``model.ndot`` 60 First time derivative of the mean motion (ignored by SGP4). 61 ``model.nddot`` 62 Second time derivative of the mean motion (ignored by SGP4). 63 ``model.bstar`` 64 Ballistic drag coefficient B* in inverse earth radii. 65 ``model.ephtype`` 66 Ephemeris type (ignored by SGP4 as determination now automatic) 67 ``model.elnum`` 68 Element number 69 ``model.inclo`` 70 Inclination in radians. 71 ``model.nodeo`` 72 Right ascension of ascending node in radians. 73 ``model.ecco`` 74 Eccentricity. 75 ``model.argpo`` 76 Argument of perigee in radians. 77 ``model.mo`` 78 Mean anomaly in radians. 79 ``model.no_kozai`` 80 Mean motion in radians per minute. 81 ``model.revnum`` 82 Revolution number at epoch [Revs] 83 84 """ 85 center = 399 86 ts = None # see __init__() 87 88 def __init__(self, line1, line2, name=None, ts=None): 89 if ts is None: 90 ts = self.ts 91 if ts is None: 92 from .api import load # avoid import loop 93 ts = EarthSatellite.ts = load.timescale() 94 95 self.name = None if name is None else name.strip() 96 satrec = Satrec.twoline2rv(line1, line2) 97 self.model = satrec 98 99 two_digit_year = satrec.epochyr 100 if two_digit_year < 57: 101 year = two_digit_year + 2000 102 else: 103 year = two_digit_year + 1900 104 105 self.epoch = ts.utc(year, 1, satrec.epochdays) 106 107 self._setup(satrec) 108 109 def _setup(self, satrec): 110 # If only I had not made __init__() specific to TLE lines, but 111 # had put them in an alternate construtor instead, this would 112 # simply have lived in __init__(). Alas! I was so young then. 113 114 self.target = -100000 - satrec.satnum 115 116 @classmethod 117 def from_satrec(cls, satrec, ts): 118 """Build an EarthSatellite from a raw sgp4 Satrec object. 119 120 This lets you provide raw numeric orbital elements instead of 121 the text of a TLE set. See :ref:`from-satrec` for detais. 122 123 """ 124 self = cls.__new__(cls) 125 self.model = satrec 126 self.name = None 127 128 # TODO: once sgp4 starts filling in epochyr and epochdays in 129 # sgp4init(), the separate epoch code here and in __init__() can 130 # be unified to always use epochyr and epochdays. 131 whole, fraction = divmod(satrec.jdsatepoch, 1.0) 132 year, month, day = compute_calendar_date(whole) 133 day += 0.5 # convert integer Julian day into Julian date float 134 self.epoch = ts.utc(year, month, day + fraction + satrec.jdsatepochF) 135 136 self._setup(satrec) 137 return self 138 139 def __str__(self): 140 return self.target_name 141 142 @property 143 def target_name(self): 144 return '{0}{1}catalog #{2} epoch {3}'.format( 145 self.name or '', 146 ' ' if self.name else '', 147 self.model.satnum, 148 self.epoch.utc_strftime(), 149 ) 150 151 def _position_and_velocity_TEME_km(self, t): 152 """Return the raw true equator mean equinox (TEME) vectors from SGP4. 153 154 Returns a tuple of NumPy arrays ``([x y z], [xdot ydot zdot])`` 155 expressed in kilometers and kilometers per second. Note that we 156 assume the TLE epoch to be a UTC date, per AIAA 2006-6753. 157 158 """ 159 sat = self.model 160 jd = t.whole 161 fraction = t.tai_fraction - t._leap_seconds() / DAY_S 162 163 if getattr(jd, 'shape', None): 164 e, r, v = sat.sgp4_array(jd, fraction) 165 messages = [SGP4_ERRORS[error] if error else None for error in e] 166 return r.T, v.T, messages 167 else: 168 error, position, velocity = sat.sgp4(jd, fraction) 169 message = SGP4_ERRORS[error] if error else None 170 return array(position), array(velocity), message 171 172 def ITRF_position_velocity_error(self, t): 173 """Deprecated: use the TEME and ITRS frame objects instead.""" 174 # TODO: can we teach frame objects to figure out that the 175 # transform TEME -> ITRS can not only skip the t.M rotation, but 176 # can also subtract the angles of their two competing z-axis 177 # rotations and call rot_z() only once instead of twice? 178 rTEME, vTEME, error = self._position_and_velocity_TEME_km(t) 179 rTEME /= AU_KM 180 vTEME /= AU_KM 181 vTEME *= DAY_S 182 rITRF, vITRF = TEME_to_ITRF(t.whole, rTEME, vTEME, 0.0, 0.0, 183 t.ut1_fraction) 184 return rITRF, vITRF, error 185 186 def _at(self, t): 187 """Compute this satellite's GCRS position and velocity at time `t`.""" 188 r, v, error = self._position_and_velocity_TEME_km(t) 189 r /= AU_KM 190 v /= AU_KM 191 v *= DAY_S 192 R = _T(TEME.rotation_at(t)) 193 r = mxv(R, r) 194 v = mxv(R, v) 195 return r, v, None, error 196 197 def find_events(self, topos, t0, t1, altitude_degrees=0.0): 198 """Return the times at which the satellite rises, culminates, and sets. 199 200 Searches between ``t0`` and ``t1``, which should each be a 201 Skyfield :class:`~skyfield.timelib.Time` object, for passes of 202 this satellite above the location ``topos`` that reach at least 203 ``altitude_degrees`` above the horizon. 204 205 Returns a tuple ``(t, events)`` whose first element is a 206 :class:`~skyfield.timelib.Time` array and whose second element 207 is an array of events: 208 209 * 0 — Satellite rose above ``altitude_degrees``. 210 * 1 — Satellite culminated and started to descend again. 211 * 2 — Satellite fell below ``altitude_degrees``. 212 213 Note that multiple culminations in a row are possible when, 214 without setting, the satellite reaches a second peak altitude 215 after descending partway down the sky from the first one. 216 217 """ 218 # First, we find the moments of maximum altitude over the time 219 # period. Some of these maxima will be negative, meaning the 220 # satellite failed to crest the horizon. 221 222 ts = t0.ts 223 at = (self - topos).at 224 half_second = 0.5 / DAY_S 225 orbits_per_minute = self.model.no_kozai / tau 226 orbits_per_day = 24 * 60 * orbits_per_minute 227 228 # Note the protection against zero orbits_per_day. 229 # TODO: why isn't 3 samples per orbit enough? 230 step_days = 0.05 / max(orbits_per_day, 1.0) 231 232 # Long-period satellites might rise each day not because of 233 # their own motion, but because the Earth rotates under them, so 234 # check position at least each quarter-day. We might need to 235 # tighten this even further if experience someday shows it 236 # missing a pass of a particular satellite. 237 if step_days > 0.25: 238 step_days = 0.25 239 240 def cheat(t): 241 """Avoid computing expensive values that cancel out anyway.""" 242 t.gast = t.tt * 0.0 243 t.M = t.MT = _identity 244 245 def altitude_at(t): 246 cheat(t) 247 return at(t).altaz()[0].degrees 248 249 altitude_at.step_days = step_days 250 tmax, altitude = find_maxima(t0, t1, altitude_at, half_second, 12) 251 if not tmax: 252 return tmax, ones_like(tmax) 253 254 # Next, filter out the maxima that are not high enough. 255 256 keepers = altitude >= altitude_degrees 257 jdmax = tmax.tt[keepers] 258 ones = ones_like(jdmax, 'uint8') 259 260 # Finally, find the rising and setting that bracket each maximum 261 # altitude. We guess that the satellite will be back below the 262 # horizon in between each pair of adjancent maxima. 263 264 def below_horizon_at(t): 265 cheat(t) 266 return at(t).altaz()[0].degrees < altitude_degrees 267 268 # The `jdo` array are the times of maxima, with their averages 269 # in between them. The start and end times are thrown in too, 270 # in case a rising or setting is lingering out between a maxima 271 # and the ends of our range. Could this perhaps still miss a 272 # stubborn rising or setting near the ends? 273 doublets = repeat(concatenate(((t0.tt,), jdmax, (t1.tt,))), 2) 274 jdo = (doublets[:-1] + doublets[1:]) / 2.0 275 276 trs, rs = _find_discrete(t0.ts, jdo, below_horizon_at, half_second, 8) 277 278 jd = concatenate((jdmax, trs.tt)) 279 v = concatenate((ones, rs * 2)) 280 281 i = jd.argsort() 282 return ts.tt_jd(jd[i]), v[i] 283 284class TEME(object): 285 """The SGP4-specific True Equator Mean Equinox frame of reference. 286 287 Described in AIAA 2006-6753 Appendix C. See :ref:`reference_frames` 288 for a guide to using Skyfield reference frames like this one. 289 290 """ 291 @staticmethod 292 def rotation_at(t): 293 theta, theta_dot = theta_GMST1982(t.whole, t.ut1_fraction) 294 angle = theta - t.gast / 24.0 * tau 295 return mxm(rot_z(angle), t.M) 296 297 # TODO: Are there any applications that will need us to include the 298 # tiny affect on velocity of the rate of change of the difference 299 # between GMST1982 and GAST? 300 301def theta_GMST1982(jd_ut1, fraction_ut1=0.0): 302 """Return the angle of Greenwich Mean Standard Time 1982 given the JD. 303 304 This angle defines the difference between the idiosyncratic True 305 Equator Mean Equinox (TEME) frame of reference used by SGP4 and the 306 more standard Pseudo Earth Fixed (PEF) frame of reference. The UT1 307 time should be provided as a Julian date. Theta is returned in 308 radians, and its velocity in radians per day of UT1 time. 309 310 From AIAA 2006-6753 Appendix C. 311 312 """ 313 t = (jd_ut1 - T0 + fraction_ut1) / 36525.0 314 g = 67310.54841 + (8640184.812866 + (0.093104 + (-6.2e-6) * t) * t) * t 315 dg = 8640184.812866 + (0.093104 * 2.0 + (-6.2e-6 * 3.0) * t) * t 316 theta = (jd_ut1 % 1.0 + fraction_ut1 + g / DAY_S % 1.0) % 1.0 * tau 317 theta_dot = (1.0 + dg / (DAY_S * 36525.0)) * tau 318 return theta, theta_dot 319 320_zero_zero_minus_one = array((0.0, 0.0, -1.0)) 321_cross120 = array((1,2,0)) 322_cross201 = array((2,0,1)) 323 324def _cross(a, b): 325 # Nearly 4x speedup over numpy cross(). TODO: Maybe move to .functions? 326 return a[_cross120] * b[_cross201] - a[_cross201] * b[_cross120] 327 328def TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp=0.0, yp=0.0, fraction_ut1=0.0): 329 """Deprecated: use the TEME and ITRS frame objects instead.""" 330 theta, theta_dot = theta_GMST1982(jd_ut1, fraction_ut1) 331 angular_velocity = multiply.outer(_zero_zero_minus_one, theta_dot) 332 333 R = rot_z(-theta) 334 335 if len(rTEME.shape) == 1: 336 rPEF = (R).dot(rTEME) 337 vPEF = (R).dot(vTEME) + _cross(angular_velocity, rPEF) 338 else: 339 rPEF = mxv(R, rTEME) 340 vPEF = mxv(R, vTEME) + _cross(angular_velocity, rPEF) 341 342 if xp == 0.0 and yp == 0.0: 343 rITRF = rPEF 344 vITRF = vPEF 345 else: 346 W = (rot_x(yp)).dot(rot_y(xp)) 347 rITRF = (W).dot(rPEF) 348 vITRF = (W).dot(vPEF) 349 return rITRF, vITRF 350