Support other observers

[jdtsmith]

Brandon Hensley suggested I take a look at zodipy. We are performing a large archival reanalysis of Spitzer Space Telescope spectroscopy, with zodiacal foreground factoring heavily into our process. Spitzer is/was on an Earth-trailing orbit and so had a unique vantage of the Earth-dust-cloud over time (different from Earth). There are ".bsp" ephemeris data for the Spitzer spacecraft available here but it isn't clear how the observer obs can be updated in ZodiPy to anything other than one of the strings in supported_observers. I can verify that jplephem can read in spk_030825_200134_220101.bsp but couldn't get much further.

Is there something we are missing? And if not, would it be straightforward to allow the use of arbitrary ephemeris positions such as a spacecraft? I imagine this will be relevant for JWST as well.

[MetinSa]

If you have the position of Spitzer available you can use the obs_pos keyword instead of obs which takes in an explicit heliocentric position. I realize that there are no examples of this keyword being used in the documentation as off now.

Regarding the ephemeris, ZodiPy currently uses the astropy.coordinates.solar_system_ephemeris which does accept urls to .bsp files, but from a quick look it doesn't look like the solar_system_ephemeris.bodies (ZodiPy's supported_observer) are updated.

I agree that it would great if arbitrary spacecrafts could be selected with the obs keyword. I do have an issue open at Astropy regarding the addition of L2 to support JWST, but I don't think this is a priority for the Astropy team. I am going to look more into this and see if there is a way to support arbitrary ephemeris positions in a nice way.

Thanks!

[jdtsmith]

Thanks. I'm certainly not an ephemeris expert. For astropy.coordinates, I believe you can use a 2-tuple of ints instead of a name, ala (10,-79):

2457598.44..2457601.50  Sun (10) -> Sirtf (-79) data_type=1

But I didn't get much further. Using jplephem and vector data might be easier. I found out the bsp file for Spitzer is type 1, which jplephem doesn't support out of the box, but the add-on spktype01 does. So something like:

from astropy.time import Time
from spktype01 import SPKType01
kernel = SPKType01.open('spk_030825_200134_220101.bsp')
pos, vel = kernel.compute_type01(10, -79, Time('2004-01-31',format='iso').to_value('jd'))
print(pos)
[-8.79526323e+07  1.10242904e+08  4.84143238e+07]

works well for vector data. It agrees perfectly (better than a mm) with the JPL Horizons web interface, but only if you select the following non-default reference frame (see this related issue):

What is the appropriate reference frame/plane to use for inputting via obs_pos? By default Horizons uses the ecliptic.

[MetinSa]

Thank you for looking more into this. I arrive at the same conclusion. My guess is that Astropy's solar system ephemeris calculations are meant mainly for major bodies, and as such only supports jplephem-friendly ephemerides.

This does look like a pretty straight-forward way to compute the position:

from astropy.time import Time
from spktype01 import SPKType01
kernel = SPKType01.open('spk_030825_200134_220101.bsp')
pos, vel = kernel.compute_type01(10, -79, Time('2004-01-31',format='iso').to_value('jd'))
print(pos)
[-8.79526323e+07  1.10242904e+08  4.84143238e+07]

but perhaps it would be difficult to generalise to arbitrary satellites? I am also not an ephemeris expert but it looks like we would need a specific implementation for the various types.

The appropriate reference frame for obs_pos is ecliptic heliocentric coordinates in AU, yes.

[jdtsmith]

The appropriate reference frame for obs_pos is ecliptic heliocentric coordinates in AU, yes.

In that case I need to rotate into the ecliptic:

from astropy.time import Time
from astropy.coordinates import SkyCoord, HeliocentricMeanEcliptic
import astropy.units as u
from spktype01 import SPKType01

kernel = SPKType01.open('/Users/jdsmith/projects/spitzer/orbit/spk_030825_200134_220101.bsp')
date = Time('2004-01-31', format='iso').to_value('jd')
pos, vel = kernel.compute_type01(10, -79, date)  # sun-centered ICRS frame, in km

# Rotate into Ecliptic frame
c = SkyCoord(*pos, unit='km', representation_type='cartesian')
c2 = c.transform_to(HeliocentricMeanEcliptic)
c2.representation_type = 'cartesian'

vec = c2.data.xyz.to(u.AU).value  # convert to AU
print(vec)
# [-0.5807906   0.80764672  0.00358533]

These agree a bit less well with Horizons, but still <30km difference.

[MetinSa]

Closing this for now since I don't have any plans in the near future to add support for Solar system bodies or satellites not supported by astropy.coordinates.solar_system_ephemeris.