diff --git a/skyfield/api.py b/skyfield/api.py index 1b59b0399..c0bcceeb9 100644 --- a/skyfield/api.py +++ b/skyfield/api.py @@ -10,7 +10,7 @@ from .constants import B1950, T0, pi, tau from .constellationlib import load_constellation_map, load_constellation_names from .iokit import Loader, load_file -from .planetarylib import PlanetaryConstants +from .planetarylib import PlanetaryConstants, PlanetaryOrientation from .positionlib import position_from_radec, position_of_radec from .starlib import Star from .sgp4lib import EarthSatellite @@ -27,7 +27,8 @@ __all__ = [ 'Angle', 'B1950', 'Distance', 'E', 'EarthSatellite', 'GREGORIAN_START', 'GREGORIAN_START_ENGLAND', - 'Loader', 'PlanetaryConstants', 'N', 'S', 'Star', 'W', + 'Loader', 'PlanetaryConstants', 'PlanetaryOrientation', + 'N', 'S', 'Star', 'W', 'T0', 'Time', 'Timescale', 'Topos', 'Velocity', 'datetime', 'iers2010', 'load', 'load_constellation_map', 'load_constellation_names', 'load_file', diff --git a/skyfield/constants.py b/skyfield/constants.py index 64a8cf951..5f4ccbb2b 100644 --- a/skyfield/constants.py +++ b/skyfield/constants.py @@ -28,5 +28,6 @@ # Time. T0 = 2451545.0 B1950 = 2433282.4235 +D_JC = 36525 # days in a Julian century C_AUDAY = C * DAY_S / AU_M diff --git a/skyfield/planetarylib.py b/skyfield/planetarylib.py index 5f0ef48a7..5655c97be 100644 --- a/skyfield/planetarylib.py +++ b/skyfield/planetarylib.py @@ -1,9 +1,9 @@ # -*- coding: utf-8 -*- """Open a BPC file, read its angles, and produce rotation matrices.""" -from numpy import array, cos, nan, sin +from numpy import array, cos, nan, sin, deg2rad, ndarray, pad from jplephem.pck import DAF, PCK -from .constants import ASEC2RAD, AU_KM, DAY_S, tau +from .constants import ASEC2RAD, AU_KM, DAY_S, tau, T0, D_JC from .data import text_pck from .functions import _T, mxv, mxm, mxmxm, rot_x, rot_y, rot_z from .units import Angle, Distance @@ -251,3 +251,150 @@ def rotation_at(self, t): # azimuth reads north-east rather than the other direction. R[1] *= -1 return R + +class PlanetaryOrientation(): + """Planet/moon orientation and shape using models from .tpc file + + Based on the models presented in pck00010.tpc, which is based on + + Archinal, B.A., A’Hearn, M.F., Bowell, E. et al. Report of the IAU + Working Group on Cartographic Coordinates and Rotational Elements: 2009. + Celest Mech Dyn Astr 109, 101–135 (2011). + https://doi.org/10.1007/s10569-010-9320-4 + + >>> from skyfield.api import load, PlanetaryConstants, PlanetaryOrientation + >>> pc = PlanetaryConstants() + >>> pc.read_text(load('pck00010.tpc')) + >>> jup_orientation = PlanetaryOrientation(pc, '599') + >>> t = load.timescale().utc(2020,1,1) + >>> jup_orientation.at(t).orientation # [RA, DEC, prime_meridian] + [268.05773343067165, 64.49711711016498, 6359115.859072568] + + """ + def __init__(self, pc, naifid): + self.a0 = None + self.d0 = None + self.W = None + self.orientation = [self.a0, self.d0, self.W] + + pck = pc.variables + + # make sure that naifid is a string number + naifid = str(naifid) + if naifid.upper().startswith('BODY'): + naifid = naifid[4:] + + self.naifid = naifid + + # not all variables are available for every body + # set to None if they aren't + self.pole_ra = self._getPCKvar(pck, naifid, 'POLE_RA') + self.pole_dec = self._getPCKvar(pck, naifid, 'POLE_DEC') + self.pm = self._getPCKvar(pck, naifid, 'PM') + self.long_axis = self._getPCKvar(pck, naifid, 'LONG_AXIS') + self.nut_prec_ra = self._getPCKvar(pck, naifid, 'NUT_PREC_RA') + self.nut_prec_dec = self._getPCKvar(pck, naifid, 'NUT_PREC_DEC') + self.nut_prec_pm = self._getPCKvar(pck, naifid, 'NUT_PREC_PM') + # NUT_PREC_ANGLES are listed under main planet barycenter number + self.nut_prec_angles = self._getPCKvar(pck, naifid[0], 'NUT_PREC_ANGLES') + self.radii = self._getPCKvar(pck, naifid, 'RADII') + + # parse NUT_PREC_ANGLES into Nx2 array for easier use later + if self.nut_prec_angles is not None: + self.nut_prec_angles_array = array([ self.nut_prec_angles[::2], + self.nut_prec_angles[1::2] ]) + + # check if POLE and PM values (orientation) are available + # many moons don't have them + if (self.pole_ra is not None and + self.pole_dec is not None and + self.pm is not None): + self.orientation_model = True + else: + self.orientation_model = False + + # check if _NUT_PREC_* values (nutation) are available + # not all bodies have it, most notably Earth + if (self.nut_prec_ra is not None and + self.nut_prec_dec is not None and + self.nut_prec_pm is not None and + self.nut_prec_angles is not None): + self.nutation_model = True + + # some bodies need to pad _NUT_PREC_ arrays with zeroes so they are the + # same length as self.nut_prec_angles_array + num_npa = self.nut_prec_angles_array.shape[1] + self.nut_prec_ra = pad(self.nut_prec_ra, (0,num_npa-len(self.nut_prec_ra))) + self.nut_prec_dec = pad(self.nut_prec_dec, (0,num_npa-len(self.nut_prec_dec))) + self.nut_prec_pm = pad(self.nut_prec_pm, (0,num_npa-len(self.nut_prec_pm))) + else: + self.nutation_model = False + + # check if RADII values (shape) are available + if (self.radii is not None): + self.shape_model = True + else: + self.shape_model = False + + # check that at least one model is available + if (not self.orientation_model and + not self.nutation_model and + not self.shape_model): + e = 'No models found for BODY' + self.naifid + \ + '. Verify that NAIF ID is correct' + raise KeyError(e) + + def _getPCKvar(self, pckvars, naifid, var): + """Do .variables[key] but return None instead of throwing an exception.""" + try: + key = 'BODY'+naifid+'_'+var # e.g. BODY399_POLE_RA + v = pckvars[key] + except KeyError: + v = None + return v + + def at(self, t): + self.t = t + self._calcOrientation() + return self + + def _calcOrientation(self, t=None): + if not self.orientation_model: + e = 'BODY' + self.naifid + ' does not have an orientation model' + raise AttributeError(e) + + if t is None: + t = self.t + else: + self.t = t + + d = t.tdb - T0 # Interval in days from the standard epoch + T = d/D_JC # Interval in Julian centuries from the standard epoch + + nut_prec_ra_sum = 0 + nut_prec_dec_sum = 0 + nut_prec_pm_sum = 0 + + # only do the nutation calculations if the necessary values are available + if self.nutation_model: + a1, a2 = self.nut_prec_angles_array + + # if time is an array, need to flip row into column + if isinstance(T, ndarray): + npa = deg2rad(a1 + a2*T.reshape(-1,1)) + else: + npa = deg2rad(a1 + a2*T) + + sum_axis = npa.ndim-1 # 0 for single time, 1 for time array + nut_prec_ra_sum = (self.nut_prec_ra*sin(npa)).sum(axis=sum_axis) + nut_prec_dec_sum = (self.nut_prec_dec*cos(npa)).sum(axis=sum_axis) + nut_prec_pm_sum = (self.nut_prec_pm*sin(npa)).sum(axis=sum_axis) + + a0 = self.pole_ra[0] + self.pole_ra[1]*T + nut_prec_ra_sum + d0 = self.pole_dec[0] + self.pole_dec[1]*T + nut_prec_dec_sum + W = self.pm[0] + self.pm[1]*d + nut_prec_pm_sum + + self.orientation = array([a0, d0, W]) + self.a0, self.d0, self.W = self.orientation + + return self.orientation