add_deflection() needs to be smarter for non-Earth observers

[psbaltar]

I've come across an oddity and I'm trying to figure out if it's an issue/bug in Skyfield or if it's a real effect that I can't figure out.

I'm calculating positions of Saturn observed from Jupiter. Astrometric positions look fine and match JPL Horizons. However, there's occasional errant points in the apparent positions. They were very obvious when I plotted the apparent positions, so I checked with Horizons and they don't match. I'm inclined to trust Horizons, but then I don't know what's going on with Skyfield.

Here's my test program that shows the issue:

from skyfield.api import load
import numpy as np
import matplotlib.pyplot as plt

# 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)
sat_ephemeris = load('sat368-named_only.bsp')
sat_ephemeris.segments.extend(ephemeris.segments)

sat = sat_ephemeris['SATURN']
jup = jup_ephemeris['JUPITER']

ts = load.timescale()

# times checked in Horizons
# marked times have errant results in Skyfield
times = [
2451874.5,
2451884.5, #
2451894.5,
2454544.5,
2454554.5, #
2454564.5,
2454574.5,
2454584.5, #
2454594.5,
2454794.5,
2454804.5, #
2454814.5 ]

# this array includes the above times
# comment it out to use the above list only
times = np.arange(2451544.5, 2454824.5, 10)

positions_astrometric = []
positions_apparent = []
for time in times:
  t = ts.tt_jd(time)

  j = jup.at(t)
  s = j.observe(sat)
  s_apparent = s.apparent()

  p_astrometric = [t.tt, *s.position.au]
  positions_astrometric.append(p_astrometric)

  p_apparent = [t.tt, *s_apparent.position.au]
  positions_apparent.append(p_apparent)
  print(p_apparent)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

T, X, Y, Z = list(zip(*positions_astrometric))
ax.plot(X, Y, Z, color='b')

T, X, Y, Z = list(zip(*positions_apparent))
ax.plot(X, Y, Z, color='r')

# results from Horizons
horizons_apparent = [
[2451874.5,2.816119885341581E+00,2.954881942270835E+00,9.948520616306117E-01],
[2451884.5,2.836025414613544E+00,2.950965039194668E+00,9.932114145858306E-01],
[2451894.5,2.856196951962868E+00,2.947748000814186E+00,9.918753685856518E-01],
[2454544.5,-9.199022886059632E+00,8.118064024004788E+00,3.809435955566133E+00],
[2454554.5,-9.298386979694197E+00,8.056312677992311E+00,3.786610180397977E+00],
[2454564.5,-9.397263345924927E+00,7.993454756839015E+00,3.763291783355337E+00],
[2454574.5,-9.495643256191899E+00,7.929483920056447E+00,3.739483553767437E+00],
[2454584.5,-9.593515311231004E+00,7.864416314379952E+00,3.715186799430429E+00],
[2454594.5,-9.690847561984429E+00,7.798245756057534E+00,3.690399660673702E+00],
[2454794.5,-1.150102167816221E+01,6.250982693174371E+00,3.094485848225297E+00],
[2454804.5,-1.158345042561566E+01,6.162888510398925E+00,3.059843118585931E+00],
[2454814.5,-1.166498799983735E+01,6.073817017832502E+00,3.024757926531031E+00] ]

T, X, Y, Z = list(zip(*horizons_apparent))
ax.plot(X, Y, Z, color='g')

plt.show()

The plot of the positions shows little spikes at a few dates. The code also includes a list of times that I ran through Horizons and I marked the errant ones. Here's the apparent XYZ coordinates the program gives me:

[2451874.5, 2.8160334634027753, 2.954952212968695, 0.9948879654471273]
[2451884.5, 2.795828156792425, 3.0058680033596077, 0.9448669500585696] # errant position
[2451894.5, 2.8563517140668604, 2.9476222331573827, 0.9918034776864725]
[2454544.5, -9.19903992139708, 8.118048704866258, 3.809427545846457]
[2454554.5, -9.327177298373325, 7.970098094246509, 3.8993412066926294] # errant position
[2454564.5, -9.397407407156875, 7.993315771684726, 3.7632273451043403]
[2454574.5, -9.495604307476361, 7.92952264211106, 3.7395004295608585]
[2454584.5, -9.583386601090348, 7.856212794574348, 3.758707128248735] # errant position
[2454594.5, -9.690852134393351, 7.798241628351907, 3.690396459709577]
[2454794.5, -11.501002341703467, 6.2510122655146105, 3.0944980717849733]
[2454804.5, -11.587438894385166, 6.143281940915633, 3.0842342923239476] # errant position
[2454814.5, -11.66539082262049, 6.0731908240432055, 3.024461953602981]

Here's the results for the same times from Horizons:

*******************************************************************************
Ephemeris / WWW_USER Sat Sep 17 10:37:05 2022 Pasadena, USA      / Horizons
*******************************************************************************
Target body name: Saturn (699)                    {source: sat441l}
Center body name: Jupiter (599)                   {source: jup365_merged}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2000-Nov-26 00:00:00.0000 TT 
Stop  time      : A.D. 2008-Dec-14 00:00:00.0000 TT 
Step-size       : DISCRETE TIME-LIST
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {W-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {W-lon(deg),Dxy(km),Dz(km)}
Center radii    : 71492.0 x 71492.0 x 66854.0 km  {Equator, meridian, pole}    
Output units    : AU-D
Output type     : LT+S CORRECTED cartesian states
Output format   : 1 (position only)
Reference frame : ICRF
*******************************************************************************
JDTT 
   X     Y     Z
*******************************************************************************
$$SOE
2451874.500000000 = A.D. 2000-Nov-26 00:00:00.0000 TT  
 X = 2.816119885341581E+00 Y = 2.954881942270835E+00 Z = 9.948520616306117E-01
2451884.500000000 = A.D. 2000-Dec-06 00:00:00.0000 TT  
 X = 2.836025414613544E+00 Y = 2.950965039194668E+00 Z = 9.932114145858306E-01
2451894.500000000 = A.D. 2000-Dec-16 00:00:00.0000 TT  
 X = 2.856196951962868E+00 Y = 2.947748000814186E+00 Z = 9.918753685856518E-01
2454544.500000000 = A.D. 2008-Mar-19 00:00:00.0000 TT  
 X =-9.199022886059632E+00 Y = 8.118064024004788E+00 Z = 3.809435955566133E+00
2454554.500000000 = A.D. 2008-Mar-29 00:00:00.0000 TT  
 X =-9.298386979694197E+00 Y = 8.056312677992311E+00 Z = 3.786610180397977E+00
2454564.500000000 = A.D. 2008-Apr-08 00:00:00.0000 TT  
 X =-9.397263345924927E+00 Y = 7.993454756839015E+00 Z = 3.763291783355337E+00
2454574.500000000 = A.D. 2008-Apr-18 00:00:00.0000 TT  
 X =-9.495643256191899E+00 Y = 7.929483920056447E+00 Z = 3.739483553767437E+00
2454584.500000000 = A.D. 2008-Apr-28 00:00:00.0000 TT  
 X =-9.593515311231004E+00 Y = 7.864416314379952E+00 Z = 3.715186799430429E+00
2454594.500000000 = A.D. 2008-May-08 00:00:00.0000 TT  
 X =-9.690847561984429E+00 Y = 7.798245756057534E+00 Z = 3.690399660673702E+00
2454794.500000000 = A.D. 2008-Nov-24 00:00:00.0000 TT  
 X =-1.150102167816221E+01 Y = 6.250982693174371E+00 Z = 3.094485848225297E+00
2454804.500000000 = A.D. 2008-Dec-04 00:00:00.0000 TT  
 X =-1.158345042561566E+01 Y = 6.162888510398925E+00 Z = 3.059843118585931E+00
2454814.500000000 = A.D. 2008-Dec-14 00:00:00.0000 TT  
 X =-1.166498799983735E+01 Y = 6.073817017832502E+00 Z = 3.024757926531031E+00
$$EOE
*******************************************************************************

[brandon-rhodes]

Well. Drat.

The astronomy library that Skyfield's core logic is based on, the NOVAS library from the USNO, is mostly designed to produce observations for an Earth observer. When writing Skyfield, I generalized the logic so that users could choose other vantage points in the Solar System.

Problem: I didn't fully generalize the deflection logic! So what I think is happening here is that it's trying to properly account for deflection, but doesn't realize that neither Saturn's gravitational deflection of light nor Jupiter's matters when we are considering Jupiter looking at Saturn. Specifically, the add_deflection() routine is going to need more information about the bodies involved in a calculation, so that it doesn't try to have bodies deflecting their own images.

Anyway, that's definitely the cause of the spikes you're seeing; if I comment out the call to add_deflection(), the spikes disappear.

Thanks for the pleasantly complete bug report, that—once I'd found all the .bsp files—let me see the same output you are seeing. I'll see if in the next week I can get some kind of fix worked out.