Central meridian longitude of planets and Jovian satellites?

[davetyp]

This is tangentially related to issue # 272. How can one extract the central meridian longitude (CML) (the body-centric longitude pointing at Earth) for a planet or a Jovian satellite? The information should be in de430 for the planets and jup310 for the Jovian satellites. What call do we make to skyfield to obtain this information? Note that for the gas giants, Jupiter in specific, one normally needs CML of System III, which is the planetary magnetosphere, not the CML of system I or II (low latitude and high latitude clouds). These three CML Systems rotate at different speeds.

This is related to issue 272 because eight Jovian satellites are tidally locked. That means if you know the CML of, e.g., Io, you know its position in orbit around Jupiter after you make a tiny adjustment for parallax. This of course doesn't help for any of the Jovian satellites that are not tidally locked.

So, how do we get CML out of Skyfield?

[brandon-rhodes]

The information should be in de430 for the planets and jup310 for the Jovian satellites.

At least as delivered in .bsp files, those ephemerides list no orientation data, alas:

$ python -m jplephem spk de430t.bsp
File type DAF/SPK and format LTL-IEEE with 15 segments:
2287184.50..2688976.50  Type 2  Venus Barycenter (2) -> Venus (299)
2287184.50..2688976.50  Type 2  Mercury Barycenter (1) -> Mercury (199)
2287184.50..2688976.50  Type 2  Earth Barycenter (3) -> Earth (399)
2287184.50..2688976.50  Type 2  Earth Barycenter (3) -> Moon (301)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Sun (10)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
2287184.50..2688976.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
2287184.50..2688976.50  Type 2  Tdb (1000000000) -> Tt (1000000001)

If you can find what kind of files the CMS information is wrapped up in by the NASA folks, however, we can look at loading that data into Skyfield. My guess would be that they distribute it as a Planetary Reference Frame?

https://rhodesmill.org/skyfield/planetary.html

Here is where the Moon kernels live:

https://naif.jpl.nasa.gov/pub/naif/LADEE/kernels/fk/

You might find Jupiter and its moons somewhere else in that file hierarchy. I will be interested to hear what you find!

[davetyp]

The information should be in de430 for the planets and jup310 for the Jovian satellites.

At least as delivered in .bsp files, those ephemerides list no orientation data, alas:

You are correct! Alas indeed.

You might find Jupiter and its moons somewhere else in that file hierarchy. I will be interested to hear what you find!

I got in touch with the author of Horizons. Turns out there is nothing too fancy for these orientations. They can be found by using the equations listed in Archinal, et al., "Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009." A PDF is available here: https://astropedia.astrogeology.usgs.gov/alfresco/d/d/workspace/SpacesStore/28fd9e81-1964-44d6-a58b-fbbf61e64e15/WGCCRE2009reprint.pdf This is what Horizons uses as a basis for calculation of orientation parameters.

Apparently that is the only solution for adding the body orientations to skyfield. At least for the outer planets and their satellites. I know the Moon kind of has its own thing going (as you mentioned). Not sure about Venus or Mercury.

For what it's worth, we are using skyfield to find transit times and fixed-elevation times (e.g., x° above the horizon and 180° AZ). The find_discrete() function is an awesomely powerful utility and works great. We ended up using astroquery's Horizon's module to pull the orientation data from JPL. Not as nice as being able to pull it from a local resource (via skyfield), but it works.

[brandon-rhodes]

Not as nice as being able to pull it from a local resource (via skyfield), but it works.

I have good news! Happily, it looks like the orientation data is available in the same pck00008.tpc same that Skyfield loads when computing Moon orientation. Actually pck00010.tpc, available a few directories up and over from the NAIF link I provided, looks to be an even more up-to-date version:

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc

It describes Jupiter’s rotation like this:

        BODY599_POLE_RA        = (   268.056595     -0.006499       0. )
        BODY599_POLE_DEC       = (    64.495303      0.002413       0. )
        BODY599_PM             = (   284.95        870.5360000      0. )
...

It looks like all Skyfield needs to learn to do is plug those values (and all of the precession numbers that follow over almost 3 dozen more lines!) into the right formulae, and it could predict where Jupiter’s pole and prime meridian is pointing at any given moment.

Let’s keep this issue open until I have time to code something up!

[davetyp]

Let’s keep this issue open until I have time to code something up!

That's excellent, glad you found that file. Sure appreciate the consideration of adding this functionality. We'd certainly use it. For our application, we would need not only Jupiter, but the eight Jovian satellites listed in that file as well. Sincerest thanks!

[psbaltar]

This is actually one of the problems that I was hoping to solve when started looking for something like Skyfield. I do some planetary imaging and one of the things I like to know is the orientation of the planet I'm looking at. In amateur planetary imaging, the de facto standard tool for getting this type of info is WinJUPOS. However, it's not easy to export the data.

The script below is what I ended up with for calculating the latitude and longitude we're looking at. (CLat and CM in WinJUPOS terminology, ObsSub-LAT and ObsSub-LON in Horizons.) It's effectively the planetographic coordinates of the center of the observed disk. I did some spot checks to make sure that WinJUPOS and Horizons agreed, then started using Horizons to verify my work.

My overall approach was to find the observer's position relative to the target, then convert to latlon coordinates.

1. get position of target relative to observer
2. reverse the position vector
3. convert reversed vector to body-fixed frame of target
4. convert to latlon

from skyfield.api import load, PlanetaryConstants
#from skyfield.api import PlanetaryOrientation #Needs PR #798

import numpy as np

'''
Returns rotation matrix to convert from ICRF to body-fixed frame

References:
Burmeister, S. (2017). Determining rotational elements of planetary bodies
https://doi.org/10.14279/depositonce-5796

Burmeister, Steffi & Willner, Konrad & Schmidt, Valentina & Oberst, Jürgen. (2018). Determination of Phobos' Rotational Parameters by an Inertial Frame Bundle Block Adjustment. Journal of Geodesy. 92.
http://dx.doi.org/10.1007/s00190-018-1112-8
'''
def Rbff(orientation):
  a0, d0, W = orientation
  Wda = np.deg2rad(np.array([W,90-d0,90+a0]))

  c1, c2, c3 = np.cos(Wda)
  s1, s2, s3 = np.sin(Wda)

  # simplified equation 2.3 (Burmeister, 2017)
  R = np.array([ [c1*c3-c2*s1*s3, c2*c3*s1+c1*s3, s1*s2],
                 [-c3*s1-c1*c2*s3, c1*c2*c3-s1*s3, c1*s2],
                 [s2*s3, -c3*s2, c2] ])

  return R

# output angles in degrees
def cart2latlon(xyz):
  x, y, z = xyz
  r = np.linalg.norm(xyz)
  theta = np.rad2deg(np.arccos(z/r))
  phi = np.rad2deg(np.arctan2(y,x))

  lat = 90-theta
  lon = phi%360
  return [lat, lon, r]

'''
Convert planetocentric to planetographic
 - input is degrees, output is degrees
 - throws exception for triaxial ellipsoid shape models

References:
Archinal, B.A., Acton, C.H., A’Hearn, M.F. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015. Celest Mech Dyn Astr 130, 22 (2018).
https://doi.org/10.1007/s10569-017-9805-5

Synder, J. (1987) Map projections: A working manual
https://doi.org/10.3133/pp1395

Clynch, J. (2002) Geodetic Coordinate Conversions. Naval Postgraduate School.
https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf
'''
def centric2graphic(centric, shape, retrograde=False):
  if shape[0] != shape[1]:
    e = 'Equatorial radii are not equal, this function only works for spheroids'
    raise ValueError(e)

  polar_radius = shape[0]
  eq_radius = shape[2]

  #f = 1-a/b
  #e2 = 2*f-f**2
  e2 = 1-((polar_radius**2)/(eq_radius**2))

  #equation 3-28 (Snyder, 1987)
  graphic_lat = np.arctan((1-e2)*np.tan(np.deg2rad(centric[0])))
  graphic_lat = np.rad2deg(graphic_lat)

  # Bodies with prograde rotation measure planetographic longitude westward,
  # and bodies with retrograde motion measure planetographic lon eastward.
  # Planetocentric longitude is always measured eastward.
  # (Archinal et al, 2018) Section 7 Definition of cartographic coordinate...
  if not retrograde:
    graphic_lon = 360-centric[1] # flip planetocentric lon to be +W
  else:
    graphic_lon = centric[1] # planetocentric lon is already +E

  return [graphic_lat, graphic_lon]

# need to load both main and planet-specific ephemerides
# https://github.com/skyfielders/python-skyfield/issues/145
ephemeris = load('de430.bsp')
jup_ephemeris = load('jup365.bsp')
jup_ephemeris.segments.extend(ephemeris.segments)

e = ephemeris['EARTH']
jup = jup_ephemeris['JUPITER']

ts = load.timescale()
t = ts.utc(2000,1,1)
print('TT:'+str(t.tt)+' TDB:'+str(t.tdb))

obs = e.at(t)
target = obs.observe(jup)

'''
#Needs PR #798
pc = PlanetaryConstants()
pc.read_text(load('pck00010.tpc'))
po = PlanetaryOrientation(pc, '599')
po.at(t-target.light_time)
orientation = po.orientation
radii = po.radii
'''

radii = [71492, 71492, 66854] #Jupiter
#Jupiter orientation at TDB=2451544.4740979616 (2451544.500742869-light_time)
orientation = [268.05720438797204, 64.49580994559105, -172.8666569354981]

v = np.array(target.position.au)*-1 # vector for target->observer
#v = (obs-target).position.au # why isn't this the same?

R = Rbff(orientation) # rotation matrix for ICRF to body-fixed frame
vbff = np.matmul(v, np.transpose(R)) # transpose so we can multiply row vector

lat, lon, d = cart2latlon(vbff) # get planetocentric (spherical) latlon
pg_lat, pg_lon = centric2graphic([lat, lon], radii) # convert to planetographic

print('ObsSub-LAT: ', round(pg_lat,6))
print('ObsSub-LON: ', round(pg_lon, 6))
print('distance: ', round(d,6))

Output of the above script:

TT:2451544.5007428704 TDB:2451544.500742869
ObsSub-LAT:  3.310617
ObsSub-LON:  340.18225
distance:  4.613423

Output from Horizons:

*******************************************************************************
Ephemeris / WWW_USER Mon Sep 26 17:58:23 2022 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Jupiter (599)                   {source: jup365_merged}
Center body name: Earth (399)                     {source: DE441}
Center-site name: GEOCENTRIC
*******************************************************************************
Start time      : A.D. 2000-Jan-01 00:00:00.0000 UT      
Stop  time      : A.D. 2000-Jan-01 00:00:00.5000 UT      
Step-size       : 0 steps
*******************************************************************************
[SNIP]                    
*******************************************************************************
 Date__(UT)__HR:MN:SC.fff Date_________JDUT     ObsSub-LON ObsSub-LAT
*********************************************************************
$$SOE
 2000-Jan-01 00:00:00.000 2451544.500000000     340.186180   3.310427
$$EOE
*******************************************************************************

There's a couple issues that I haven't figured out yet:

1. (obs-target).position.au This isn't the same as reversing the obs->target vector and I don't understand why. Negating the obs->target vector produces the correct result, but obs-target doesn't.
2. I haven't verified it's accuracy on moons yet. However, some of the IAU shape models for moons are ellipsoids. Archinal, et al mentions inconsistencies in whether those bodies should be treated as triaxial ellipsoids or spheroids. So, the planetographic conversion function may need to be modified to handle them. The one I have above only handles spheroids and throws an exception for ellipsoids. I also haven't looked into how they define the X and Y axes for ellipsoidal bodies.

[brandon-rhodes]

Thank you for working out how this calculation looks in Python! In October I should have time to add a routine that does this to Skyfield, based on your example here. Here are a few quick notes, though, while I have a few minutes:

• Subtracting obs-target does not, I don't think, produce a useful vector? The obs vector is Solar System barycenter → earth, target is earth → juptier, so subtracting them produces—the position that Jupiter would have in the Solar System if when viewed from Earth it were 180° away from its current position?

• Happily, you don't need to wrap np.array() around target.position.au, as it's already a NumPy array.

• Instead of matmul(), you can do matrix × vector and have all the array broadcasting done correctly with:

    from skyfield.functions import mxv
    vbff = mxv(R, v)

• There might be a few other routines in skyfield.functions that you might be able to use; for example, one of them might do most of the work of cart2latlon()?

[psbaltar]

I made some changes to use functions in skyfield.functions as well as vectorize things. The main thing was updating Rbff() for it.

from skyfield.api import load, PlanetaryConstants
#from skyfield.api import PlanetaryOrientation #Needs PR #798

from skyfield.functions import mxm, mxv, rot_x, rot_z, to_spherical

import numpy as np

'''
Returns rotation matrix to convert from ICRF to body-fixed frame

References:
Burmeister, S. (2017). Determining rotational elements of planetary bodies
https://doi.org/10.14279/depositonce-5796

Burmeister, Steffi & Willner, Konrad & Schmidt, Valentina & Oberst, Jürgen. (2018). Determination of Phobos' Rotational Parameters by an Inertial Frame Bundle Block Adjustment. Journal of Geodesy. 92.
http://dx.doi.org/10.1007/s00190-018-1112-8
'''
def Rbff(orientation):
  a0, d0, W = orientation
  adW = np.deg2rad(np.array([90+a0, 90-d0, W]))

  # equation 2.3 (Burmeister, 2017)
  # arguments negated to work with Skyfield rotation matrix functions
  R = mxm(rot_z(-adW[2]), mxm(rot_x(-adW[1]), rot_z(-adW[0])))
  return R

'''
Convert planetocentric to planetographic
 - input is degrees, output is degrees
 - throws exception for triaxial ellipsoid shape models

References:
Archinal, B.A., Acton, C.H., A’Hearn, M.F. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015. Celest Mech Dyn Astr 130, 22 (2018).
https://doi.org/10.1007/s10569-017-9805-5

Synder, J. (1987) Map projections: A working manual
https://doi.org/10.3133/pp1395

Clynch, J. (2002) Geodetic Coordinate Conversions. Naval Postgraduate School.
https://www.oc.nps.edu/oc2902w/coord/coordcvt.pdf
'''
def centric2graphic(lat, lon, shape, retrograde=False):
  if shape[0] != shape[1]:
    e = 'Equatorial radii are not equal, this function only works for spheroids'
    raise ValueError(e)

  polar_radius = shape[0]
  eq_radius = shape[2]

  #f = 1-a/b
  #e2 = 2*f-f**2
  e2 = 1-((polar_radius**2)/(eq_radius**2))

  #equation 3-28 (Snyder, 1987)
  graphic_lat = np.arctan((1-e2)*np.tan(np.deg2rad(lat)))
  graphic_lat = np.rad2deg(graphic_lat)

  # Bodies with prograde rotation measure planetographic longitude westward,
  # and bodies with retrograde rotation measure planetographic lon eastward.
  # Planetocentric longitude is always measured eastward.
  # (Archinal et al, 2018) Section 7 Definition of cartographic coordinate...
  if not retrograde:
    graphic_lon = 360-lon # flip planetocentric lon to be +W
  else:
    graphic_lon = lon # planetocentric lon is already +E

  return [graphic_lat, graphic_lon]

# need to load both main and planet-specific ephemerides
# https://github.com/skyfielders/python-skyfield/issues/145
ephemeris = load('de430.bsp')
jup_ephemeris = load('jup365.bsp')
jup_ephemeris.segments.extend(ephemeris.segments)

e = ephemeris['EARTH']
jup = jup_ephemeris['JUPITER']

ts = load.timescale()
t = ts.ut1(2000)
#t = ts.ut1(range(2000,2020))

obs = e.at(t)
target = obs.observe(jup)

'''
#Needs PR #798
pc = PlanetaryConstants()
pc.read_text(load('pck00010.tpc'))
po = PlanetaryOrientation(pc, '599')
po.at(ts.tt_jd(t.tt-target.light_time))
orientation = po.orientation
radii = po.radii
'''

radii = [71492, 71492, 66854] #Jupiter
#Jupiter orientation at TDB=2451544.4740979616 (2451544.500742869-light_time)
orientation = [268.05720438797204, 64.49580994559105, -172.8666569354981]

v = -1*target.position.au # vector for target->observer

R = Rbff(orientation)
vbff = mxv(R, v)

lat, lon = np.rad2deg(to_spherical(vbff)[1:]) # planetocentric coordinates
pg_lat, pg_lon = centric2graphic(lat, lon, radii) # planetographic

with np.printoptions(suppress=True):
  data = np.array([ t.ut1, np.round(pg_lon,6), np.round(pg_lat,6) ]).T
  print(data)

[brandon-rhodes]

Thank you for the update! In a couple of weeks, I will hopefully be able to reply with a plan for getting this logic into Skyfield.