1# -*- coding: utf-8 -*-
2"""Open a BPC file, read its angles, and produce rotation matrices."""
3
4from numpy import array, cos, nan, sin
5from jplephem.pck import DAF, PCK
6from .constants import ASEC2RAD, AU_KM, DAY_S, tau
7from .data import text_pck
8from .functions import _T, mxv, mxm, mxmxm, rot_x, rot_y, rot_z
9from .units import Angle, Distance
10from .vectorlib import VectorFunction
11
12_TEXT_MAGIC_NUMBERS = b'KPL/FK', b'KPL/PCK'
13_NAN3 = array((nan, nan, nan))
14_halftau = tau / 2.0
15_quartertau = tau / 4.0
16
17class PlanetaryConstants(object):
18    """Planetary constants manager.
19
20    You can use this class to build working models of Solar System
21    bodies by loading text planetary constant files and binary
22    orientation kernels.  For a full description of how to use this, see
23    :doc:`planetary`.
24
25    """
26    def __init__(self):
27        self.variables = {}
28        self._binary_files = []
29        self._segment_map = {}
30
31    @property
32    def assignments(self):  # deprecated original name for the variables dict
33        return self.variables
34
35    def read_text(self, file):
36        """Read frame variables from a KPL/FK file.
37
38        Appropriate files will typically have the extension ``.tf`` or
39        ``.tpc`` and will define a series of names and values that will
40        be loaded into this object's ``.variables`` dictionary.
41
42        >>> from skyfield.api import load
43        >>> pc = PlanetaryConstants()
44        >>> pc.read_text(load('moon_080317.tf'))
45        >>> pc.variables['FRAME_31006_NAME']
46        'MOON_PA_DE421'
47
48        """
49        file.seek(0)
50        try:
51            if not file.read(7).startswith(_TEXT_MAGIC_NUMBERS):
52                raise ValueError('file must start with one of the patterns:'
53                                 ' {0}'.format(_TEXT_MAGIC_NUMBERS))
54            text_pck.load(file, self.variables)
55        finally:
56            file.close()
57
58    def read_binary(self, file):
59        """Read binary segments descriptions from a DAF/PCK file.
60
61        Binary segments live in ``.bpc`` files and predict how a body
62        like a planet or moon will be oriented on a given date.
63
64        """
65        file.seek(0)
66        if file.read(7) != b'DAF/PCK':
67            raise ValueError('file must start with the bytes "DAF/PCK"')
68        pck = PCK(DAF(file))
69        self._binary_files.append(pck)
70        for segment in pck.segments:
71            self._segment_map[segment.body] = segment
72
73    def _get_assignment(self, key):
74        """Do .variables[key] but with a pretty exception on failure."""
75        try:
76            return self.variables[key]
77        except KeyError:
78            e = ValueError(_missing_name_message.format(key))
79            e.__cause__ = None
80            raise e
81
82    def build_frame_named(self, name):
83        """Given a frame name, return a :class:`Frame` object."""
84        integer = self._get_assignment('FRAME_{0}'.format(name))
85        return self.build_frame(integer)
86
87    def build_frame(self, integer, _segment=None):
88        """Given a frame integer code, return a :class:`Frame` object."""
89        center = self._get_assignment('FRAME_{0}_CENTER'.format(integer))
90        spec = self.variables.get('TKFRAME_{0}_SPEC'.format(integer))
91        if spec is None:
92            matrix = None
93        else:
94            if spec == 'ANGLES':
95                angles = self.variables['TKFRAME_{0}_ANGLES'.format(integer)]
96                axes = self.variables['TKFRAME_{0}_AXES'.format(integer)]
97                units = self.variables['TKFRAME_{0}_UNITS'.format(integer)]
98                scale = _unit_scales[units]
99                matrix = 1,0,0, 0,1,0, 0,0,1
100                matrix = array(matrix)
101                matrix.shape = 3, 3
102                for angle, axis in list(zip(angles, axes)):
103                    rot = _rotations[axis]
104                    matrix = mxm(rot(angle * scale), matrix)
105            elif spec == 'MATRIX':
106                matrix = self.variables['TKFRAME_{0}_MATRIX'.format(integer)]
107                matrix = array(matrix)
108                matrix.shape = 3, 3
109            else:
110                raise NotImplementedError('spec %r not yet implemented' % spec)
111            relative = self.variables['TKFRAME_{0}_RELATIVE'.format(integer)]
112            integer = self.variables['FRAME_{0}'.format(relative)]
113
114        if _segment is None:
115            segment = self._segment_map.get(integer)
116        else:
117            segment = _segment
118
119        if segment is None:
120            raise LookupError('you have not yet loaded a binary PCK file that'
121                              ' has a segment for frame {0}'.format(integer))
122        assert segment.frame == 1  # base frame should be ITRF/J2000
123        return Frame(center, segment, matrix)
124
125    def build_latlon_degrees(self, frame, latitude_degrees, longitude_degrees,
126                             elevation_m=0.0):
127        """Build an object representing a location on a body's surface."""
128        lat = Angle.from_degrees(latitude_degrees)
129        lon = Angle.from_degrees(longitude_degrees)
130        radii = self._get_assignment('BODY{0}_RADII'.format(frame.center))
131        if not radii[0] == radii[1] == radii[2]:
132            raise ValueError('only spherical bodies are supported,'
133                             ' but the radii of this body are: %s' % radii)
134        au = (radii[0] + elevation_m * 1e-3) / AU_KM
135        distance = Distance(au)
136        return PlanetTopos.from_latlon_distance(frame, lat, lon, distance)
137
138_rotations = None, rot_x, rot_y, rot_z
139_unit_scales = {'ARCSECONDS': ASEC2RAD}
140_missing_name_message = """unknown planetary constant {0!r}
141
142You should either use this object's `.read_text()` method to load an
143additional "*.tf" PCK text file that defines the missing name, or
144manually provide a value by adding the name and value to the this
145object's `.variables` dictionary."""
146
147class Frame(object):
148    """Planetary constants frame, for building rotation matrices."""
149
150    def __init__(self, center, segment, matrix):
151        self.center = center
152        self._segment = segment
153        self._matrix = matrix
154
155    def rotation_at(self, t):
156        """Return the rotation matrix for this frame at time ``t``."""
157        ra, dec, w = self._segment.compute(t.tdb, 0.0, False)
158        R = mxm(rot_z(-w), mxm(rot_x(-dec), rot_z(-ra)))
159        if self._matrix is not None:
160            R = mxm(self._matrix, R)
161        return R
162
163    def rotation_and_rate_at(self, t):
164        """Return rotation and rate matrices for this frame at time ``t``.
165
166        The rate matrix returned is in units of angular motion per day.
167
168        """
169        components, rates = self._segment.compute(t.whole, t.tdb_fraction, True)
170        ra, dec, w = components
171        radot, decdot, wdot = rates
172
173        R = mxm(rot_z(-w), mxm(rot_x(-dec), rot_z(-ra)))
174
175        zero = w * 0.0
176        ca = cos(w)
177        sa = sin(w)
178        u = cos(dec)
179        v = -sin(dec)
180
181        domega0 = wdot + u * radot
182        domega1 = ca * decdot - sa * v * radot
183        domega2 = sa * decdot + ca * v * radot
184
185        drdtrt = array((
186            (zero, domega0, domega2),
187            (-domega0, zero, domega1),
188            (-domega2, -domega1, zero),
189        ))
190
191        dRdt = mxm(drdtrt, R)
192
193        if self._matrix is not None:
194            R = mxm(self._matrix, R)
195            dRdt = mxm(self._matrix, dRdt)
196
197        return R, dRdt * DAY_S
198
199class PlanetTopos(VectorFunction):
200    """Location that rotates with the surface of another Solar System body.
201
202    The location can either be on the surface of the body, or in some
203    other fixed position that rotates with the body's surface.
204
205    """
206    def __init__(self, frame, position_au):
207        self.center = frame.center
208        self._frame = frame
209        self._position_au = position_au
210
211    @classmethod
212    def from_latlon_distance(cls, frame, latitude, longitude, distance):
213        r = array((distance.au, 0.0, 0.0))
214        r = mxv(rot_z(longitude.radians), mxv(rot_y(-latitude.radians), r))
215
216        self = cls(frame, r)
217        self.latitude = latitude
218        self.longitude = longitude
219        return self
220
221    @property
222    def target(self):
223        # When used as a vector function, this planetary geographic
224        # location computes positions from the planet's center to
225        # itself.  (This is a property, rather than an attribute, to
226        # avoid a circular reference that delays garbage collection.)
227        return self
228
229    def _at(self, t):
230        # Since `_position_au` has zero velocity in this reference
231        # frame, velocity includes a `dRdt` term but not an `R` term.
232        R, dRdt = self._frame.rotation_and_rate_at(t)
233        r = mxv(_T(R), self._position_au)
234        v = mxv(_T(dRdt), self._position_au)
235        return r, v, None, None
236
237    def rotation_at(self, t):
238        """Compute the altazimuth rotation matrix for this location’s sky."""
239        R = mxmxm(
240            # TODO: Figure out how to produce this rotation directly
241            # from _position_au, to support situations where we were not
242            # given a latitude and longitude.  If that is not feasible,
243            # then at least cache the product of these first two matrices.
244            rot_y(_quartertau - self.latitude.radians),
245            rot_z(_halftau - self.longitude.radians),
246            self._frame.rotation_at(t),
247        )
248        # TODO:
249        # Can clockwise be turned into counterclockwise through any
250        # possible rotation?  For now, flip the sign of y so that
251        # azimuth reads north-east rather than the other direction.
252        R[1] *= -1
253        return R
254