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