(ANSWERED) Observe a planet from a satellite

[craigim]

I have a satellite with a camera on it. I'd like to take a picture of Jupiter. I cannot figure out what RA/Dec direction I should look in.

I start by picking a time:

from dateutil.parser import parse
import skyfield.api as sf
import numpy as np

iso_time_string = "2021-03-14T15:09:26.5359Z"
ts = sf.load.timescale(builtin=True)
epoch = parse(iso_time_string, fuzzy=True)
t = ts.from_datetime(epoch)

I then set an inertial position of the satellite (I don't have a TLE to initialize an EarthSatellite object, just an ECI position and time)

mean_rad = 6371008.8
altitude =  540000.0
eci_m = np.random.normal(size=(3))
eci_m /= np.sqrt(np.sum(eci_m*eci_m))
eci_m *= (mean_rad + altitude)
eci_x_m, eci_y_m, eci_z_m = eci_m
p = sf.Distance(m=[eci_x_m, eci_y_m, eci_z_m])
sat_position = Geocentric(p.au, t=t)
observer = sf.wgs84.subpoint(sat_position).at(t)

I get the ephemeris of Jupiter

planets = sf.load("de440.bsp")
target = planets["JUPITER BARYCENTER"]

Here is what I'd like to do:

direction = sat_position.observe(target).apparent()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-68-f3693bf3c06e> in <module>
----> 1 sat_position.observe(target).apparent()

AttributeError: 'Geocentric' object has no attribute 'observe'

direction = observer.observe(target).apparent()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-54-e4f28b2d8a31> in <module>
----> 1 observer.observe(target).apparent()

AttributeError: 'Geocentric' object has no attribute 'observe'

The best I can do is

direction = (target - observer).radec()

but that gives me the instantaneous position. Given that it takes 43 minutes for light to get from Jupiter to my camera, my camera has moved almost half an orbit and Jupiter will have moved with respect to the background stars.

Of course, this doesn't work either.

direction = (target - observer).apparent()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-70-08e4737a0eb3> in <module>
----> 1 direction = (target - observer).apparent()

AttributeError: 'Geometric' object has no attribute 'apparent'

Do I have to convert my inertial position to Earth fixed? If so, what is the best way to do that? Ultimately, isn't RA/Dec an inertial coordinate? I've tried the same thing using the ICRF class and get the same results.

[JoshPaterson]

Hi! I think I was able to figure out what is going wrong. Here is how I edited your example:

from dateutil.parser import parse
import skyfield.api as sf
import numpy as np
from skyfield.positionlib import Geocentric, Barycentric

planets = sf.load("de440.bsp")
jupiter = planets["JUPITER BARYCENTER"]
earth = planets["EARTH"]

iso_time_string = "2021-03-14T15:09:26.5359Z"
ts = sf.load.timescale(builtin=True)
epoch = parse(iso_time_string, fuzzy=True)
t = ts.from_datetime(epoch)

mean_rad = 6371008.8
altitude =  540000.0
eci_m = np.random.normal(size=(3))
eci_m /= np.sqrt(np.sum(eci_m*eci_m))
eci_m *= (mean_rad + altitude)
eci_x_m, eci_y_m, eci_z_m = eci_m
p = sf.Distance(m=[eci_x_m, eci_y_m, eci_z_m])
geocentric_sat_position = Geocentric(p.au, t=t)
observer = sf.wgs84.subpoint(geocentric_sat_position).at(t)

barycentric_sat_position = Barycentric(earth.at(t).position.au + p.au, t=t)
barycentric_sat_position._ephemeris = planets
direction = barycentric_sat_position.observe(jupiter).apparent()

The observe method is not available on on Geocentric positions, so we need a Barycentric position instead. After doing that I started getting an error about a missing _ephemeris attribute, so I also added the line barycentric_sat_position._ephemeris = planets. I think this was missing because we made the Barycentric position manually.

Let me know if that produces the numbers you are expecting!

[craigim]

When I run your code, I get the following result:

direction.radec()

.../skyfield/functions.py:85: RuntimeWarning: invalid value encountered in double_scalars
  phi = arctan2(y, x) % tau

(<Angle nan>, <Angle nan>, <Distance nan au>)

[JoshPaterson]

I should have made sure actual numbers were being produced instead of just checking that direction had the right object type! I'm not entirely certain why this happens, but I did find out that getting rid of the .apparent() method call makes direction have actual numbers:

direction = barycentric_sat_position.observe(jupiter)
print(direction.radec(epoch='date'))

produces this for me:

(<Angle 21h 29m 01.91s>, <Angle -15deg 30' 04.5">, <Distance 5.85848 au>)

This removes the correction for gravitational deflection, but that should be quite small, so I would guess this is enough to get Jupiter in the frame. I'm trying to think of why the apparent method wouldn't work in this case and nothing is coming to mind. I think Skyfield's maintainer (@brandon-rhodes) is out of his office for a bit, but when he comes back he might know what needs to change in the script to allow deflection to be taken into account. Also if I think of anything else I'll post it here.

[brandon-rhodes]

Thanks for working to answer this question in my absence, @JoshPaterson, you made helpful progress!

@craigim, the only remaining snag with the script is that it provides no velocity for the satellite, and with a nan unknown velocity, apparent() can't apply the affect of aberration to the coordinates. While it's probably not what you want, setting a zero velocity at least results in coordinates being returned:

import numpy as np
import skyfield.api as sf
from dateutil.parser import parse
from skyfield.positionlib import Barycentric

planets = sf.load("de440.bsp")
jupiter = planets["JUPITER BARYCENTER"]
earth = planets["EARTH"]

iso_time_string = "2021-03-14T15:09:26.5359Z"
ts = sf.load.timescale(builtin=True)
epoch = parse(iso_time_string, fuzzy=True)
t = ts.from_datetime(epoch)

mean_rad = 6371008.8
altitude =  540000.0
eci_m = np.random.normal(size=(3))
eci_m /= np.sqrt(np.sum(eci_m*eci_m))
eci_m *= (mean_rad + altitude)
eci_x_m, eci_y_m, eci_z_m = eci_m
p_au = sf.Distance(m=[eci_x_m, eci_y_m, eci_z_m]).au
v_au_per_d = p_au * 0.0
e = earth.at(t)

barycentric_sat_position = Barycentric(
    e.position.au + p_au,
    e.velocity.au_per_d + v_au_per_d,
    t=t,
)
barycentric_sat_position._ephemeris = planets
direction = barycentric_sat_position.observe(jupiter).apparent()
print(direction.radec())

The result:

(<Angle 21h 27m 51.84s>, <Angle -15deg 35' 36.4">, <Distance 5.85842 au>)

You will presumably want to supply the actual satellite velocity, if you care about an effect of the size of aberration (about 20 arcseconds).

Let us know whether this gets your script working!

[craigim]

Sorry for the late reply. Adding in the velocities is exactly what was needed. Thanks!