From fac654f47480f3f1900c9c9db1643d44f20577a7 Mon Sep 17 00:00:00 2001 From: Philip Baltar Date: Fri, 23 Sep 2022 17:17:28 -0700 Subject: [PATCH 1/5] first pass at planetary orientation --- skyfield/api.py | 5 +- skyfield/constants.py | 1 + skyfield/planetarylib.py | 101 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 103 insertions(+), 4 deletions(-) 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..61f073024 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 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,100 @@ def rotation_at(self, t): # azimuth reads north-east rather than the other direction. R[1] *= -1 return R + +class PlanetaryOrientation(): + """Planet/moon orientation using models from .tpc file + + >>> 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 + + self.naifid = naifid + + self.pole_ra = pck['BODY'+naifid+'_POLE_RA'] + self.pole_dec = pck['BODY'+naifid+'_POLE_DEC'] + self.pm = pck['BODY'+naifid+'_PM'] + self.long_axis = self._getPCKvar(pck, 'BODY'+naifid+'_LONG_AXIS') + self.nut_prec_ra = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_RA') + self.nut_prec_dec = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_DEC') + self.nut_prec_pm = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_PM') + # NUT_PREC_ANGLES are listed under main planet barycenter number + self.nut_prec_angles = self._getPCKvar(pck, 'BODY'+naifid[0]+'_NUT_PREC_ANGLES') + self.radii = pck['BODY'+naifid+'_RADII'] + + # check if _NUT_PREC_* values are available + 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 = True + else: + self.nutation = False + + + # parse nut_prec_angles into list of tuples + if self.nutation: + self.nut_prec_angles_pairs = list(zip( self.nut_prec_angles[::2], + self.nut_prec_angles[1::2] )) + + def _getPCKvar(self, pck, var_name): + try: + v = pck[var_name] + except: + v = None + return v + + def at(self, t): + self.t = t + self._calcOrientation() + return self + + def _calcOrientation(self, t=None): + 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 + + if self.nutation: + for i in range(len(self.nut_prec_ra)): + a1, a2 = self.nut_prec_angles_pairs[i] + npa = a1 + a2*T + npa_rad = deg2rad(npa) + + nut_prec_ra_sum += self.nut_prec_ra[i]*sin(npa_rad) + nut_prec_dec_sum += self.nut_prec_dec[i]*cos(npa_rad) + nut_prec_dec_sum += self.nut_prec_pm[i]*sin(npa_rad) + + + 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 = [a0, d0, W] + self.a0, self.d0, self.W = self.orientation + + return self.orientation + + + + From a5e62f8b2f09f90c056a57d1287be6a191186a43 Mon Sep 17 00:00:00 2001 From: Philip Baltar Date: Sat, 24 Sep 2022 11:43:34 -0700 Subject: [PATCH 2/5] more comments and error handling --- skyfield/planetarylib.py | 93 +++++++++++++++++++++++++++++----------- 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/skyfield/planetarylib.py b/skyfield/planetarylib.py index 61f073024..b485d5b40 100644 --- a/skyfield/planetarylib.py +++ b/skyfield/planetarylib.py @@ -253,7 +253,14 @@ def rotation_at(self, t): return R class PlanetaryOrientation(): - """Planet/moon orientation using models from .tpc file + """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() @@ -272,38 +279,70 @@ def __init__(self, pc, naifid): 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 - self.pole_ra = pck['BODY'+naifid+'_POLE_RA'] - self.pole_dec = pck['BODY'+naifid+'_POLE_DEC'] - self.pm = pck['BODY'+naifid+'_PM'] - self.long_axis = self._getPCKvar(pck, 'BODY'+naifid+'_LONG_AXIS') - self.nut_prec_ra = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_RA') - self.nut_prec_dec = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_DEC') - self.nut_prec_pm = self._getPCKvar(pck, 'BODY'+naifid+'_NUT_PREC_PM') + # 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, 'BODY'+naifid[0]+'_NUT_PREC_ANGLES') - self.radii = pck['BODY'+naifid+'_RADII'] + self.nut_prec_angles = self._getPCKvar(pck, naifid[0], 'NUT_PREC_ANGLES') + self.radii = self._getPCKvar(pck, naifid, 'RADII') + + # 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 are available + # 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 = True + self.nutation_model = True else: - self.nutation = False + self.nutation_model = False - - # parse nut_prec_angles into list of tuples - if self.nutation: + # 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) + + # parse NUT_PREC_ANGLES into list of tuples for easier use later + if self.nutation_model: self.nut_prec_angles_pairs = list(zip( self.nut_prec_angles[::2], self.nut_prec_angles[1::2] )) - def _getPCKvar(self, pck, var_name): + def _getPCKvar(self, pckvars, naifid, var): + """Do .variables[key] but return None instead of throwing an exception.""" try: - v = pck[var_name] - except: + key = 'BODY'+naifid+'_'+var # e.g. BODY399_POLE_RA + v = pckvars[key] + except KeyError: v = None return v @@ -313,6 +352,10 @@ def at(self, t): 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: @@ -325,7 +368,12 @@ def _calcOrientation(self, t=None): nut_prec_dec_sum = 0 nut_prec_pm_sum = 0 - if self.nutation: + # only do the nutation calculations if the necessary values are available + if self.nutation_model: + # Is keying to NUT_PREC_RA really the correct thing to do? + # This keeps everything in a single loop and there doesn't seem to be + # any cases NUT_PREC variables for the same body having different numbers + # of values. for i in range(len(self.nut_prec_ra)): a1, a2 = self.nut_prec_angles_pairs[i] npa = a1 + a2*T @@ -335,7 +383,6 @@ def _calcOrientation(self, t=None): nut_prec_dec_sum += self.nut_prec_dec[i]*cos(npa_rad) nut_prec_dec_sum += self.nut_prec_pm[i]*sin(npa_rad) - 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 @@ -344,7 +391,3 @@ def _calcOrientation(self, t=None): self.a0, self.d0, self.W = self.orientation return self.orientation - - - - From a270d71c47b28c0bd9578e560734c897e9586ea9 Mon Sep 17 00:00:00 2001 From: Philip Baltar Date: Mon, 26 Sep 2022 18:23:25 -0700 Subject: [PATCH 3/5] corrected PM nutation sum --- skyfield/planetarylib.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skyfield/planetarylib.py b/skyfield/planetarylib.py index b485d5b40..00ecefe9e 100644 --- a/skyfield/planetarylib.py +++ b/skyfield/planetarylib.py @@ -381,7 +381,7 @@ def _calcOrientation(self, t=None): nut_prec_ra_sum += self.nut_prec_ra[i]*sin(npa_rad) nut_prec_dec_sum += self.nut_prec_dec[i]*cos(npa_rad) - nut_prec_dec_sum += self.nut_prec_pm[i]*sin(npa_rad) + nut_prec_pm_sum += self.nut_prec_pm[i]*sin(npa_rad) 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 From ce05d630802a591062bce1ee0ac6e89ce9065e23 Mon Sep 17 00:00:00 2001 From: Philip Baltar Date: Thu, 29 Sep 2022 17:22:10 -0700 Subject: [PATCH 4/5] vectorized time-related calculations --- skyfield/planetarylib.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/skyfield/planetarylib.py b/skyfield/planetarylib.py index 00ecefe9e..cffdb3e3e 100644 --- a/skyfield/planetarylib.py +++ b/skyfield/planetarylib.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """Open a BPC file, read its angles, and produce rotation matrices.""" -from numpy import array, cos, nan, sin, deg2rad +from numpy import array, cos, nan, sin, deg2rad, ndarray from jplephem.pck import DAF, PCK from .constants import ASEC2RAD, AU_KM, DAY_S, tau, T0, D_JC from .data import text_pck @@ -332,10 +332,10 @@ def __init__(self, pc, naifid): '. Verify that NAIF ID is correct' raise KeyError(e) - # parse NUT_PREC_ANGLES into list of tuples for easier use later + # parse NUT_PREC_ANGLES into Nx2 array for easier use later if self.nutation_model: - self.nut_prec_angles_pairs = list(zip( self.nut_prec_angles[::2], - self.nut_prec_angles[1::2] )) + self.nut_prec_angles_array = array([ self.nut_prec_angles[::2], + self.nut_prec_angles[1::2] ]) def _getPCKvar(self, pckvars, naifid, var): """Do .variables[key] but return None instead of throwing an exception.""" @@ -370,24 +370,24 @@ def _calcOrientation(self, t=None): # only do the nutation calculations if the necessary values are available if self.nutation_model: - # Is keying to NUT_PREC_RA really the correct thing to do? - # This keeps everything in a single loop and there doesn't seem to be - # any cases NUT_PREC variables for the same body having different numbers - # of values. - for i in range(len(self.nut_prec_ra)): - a1, a2 = self.nut_prec_angles_pairs[i] - npa = a1 + a2*T - npa_rad = deg2rad(npa) - - nut_prec_ra_sum += self.nut_prec_ra[i]*sin(npa_rad) - nut_prec_dec_sum += self.nut_prec_dec[i]*cos(npa_rad) - nut_prec_pm_sum += self.nut_prec_pm[i]*sin(npa_rad) + 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 = [a0, d0, W] + self.orientation = array([a0, d0, W]) self.a0, self.d0, self.W = self.orientation return self.orientation From 6c7ba4ddb647f0a1e9ba8e608db0eecccc90cb09 Mon Sep 17 00:00:00 2001 From: Philip Baltar Date: Mon, 3 Oct 2022 21:47:35 -0700 Subject: [PATCH 5/5] account for truncated *_NUT_PREC_* arrays --- skyfield/planetarylib.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/skyfield/planetarylib.py b/skyfield/planetarylib.py index cffdb3e3e..5655c97be 100644 --- a/skyfield/planetarylib.py +++ b/skyfield/planetarylib.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """Open a BPC file, read its angles, and produce rotation matrices.""" -from numpy import array, cos, nan, sin, deg2rad, ndarray +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, T0, D_JC from .data import text_pck @@ -299,6 +299,11 @@ def __init__(self, pc, naifid): 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 @@ -315,6 +320,13 @@ def __init__(self, pc, naifid): 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 @@ -332,11 +344,6 @@ def __init__(self, pc, naifid): '. Verify that NAIF ID is correct' raise KeyError(e) - # parse NUT_PREC_ANGLES into Nx2 array for easier use later - if self.nutation_model: - self.nut_prec_angles_array = array([ self.nut_prec_angles[::2], - self.nut_prec_angles[1::2] ]) - def _getPCKvar(self, pckvars, naifid, var): """Do .variables[key] but return None instead of throwing an exception.""" try: @@ -377,7 +384,7 @@ def _calcOrientation(self, t=None): 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)