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