Observing celestials from Gaia (and other satellites), with light deflection?

[thunkii]

Hello. For a project I am doing, I need to calculate apparent positions of various asteroids from the Gaia spacecraft.

Right now, to get a working .observer(), I have basically downloaded the list of barycentric positions at each time I need from Horizons, and put them in a Barycentric object like so:

gaiacsv=pd.read_csv('gaia_ephem.txt',index_col=False)
gaiatimes=ts.tdb_jd(gaiacsv['time_jd'].to_numpy())
gaia_position=Barycentric(gaiacsv[['x','y','z']].to_numpy().T,gaiacsv[['vx','vy','vz']].to_numpy().T,gaiatimes,0)

Now I need to observe asteroids from Gaia, which are provided as a list of Barycentric positions at the times we need (borrowing ThingInSpace from a previous discussion answer):

class ThingInSpace(VectorFunction):
    def __init__(self, barypos):
        # Value 0 corresponds to the solar system barycenter according to this: https://rhodesmill.org/skyfield/api-vectors.html#skyfield.vectorlib.VectorFunction.center
        self.center = 0
        self.barypos = barypos

    @property
    def target(self):
        # This is a property to avoid circular references as described here: https://github.com/skyfielders/python-skyfield/blob/master/skyfield/planetarylib.py#L222
        return self

    def _at(self, t):
        cantimes=self.barypos.t
        #print(t,cantimes)
        #we're getting asked for times not quite at the canonical times, so add a velocity fudge factor (first order integration is fine, it's like 30 minutes)
        # Velocity vector
        v = self.barypos.velocity.au_per_d
        # Position vector
        p = self.barypos.xyz.au+(t-cantimes)*v
        return p, v, self.center, "Thing in space"

(... many lines ...)

        smallbody_gaia=build_position(asteroid_pos[i,6198:10318,:].T, asteroid_vel[i,6198:10318,:].T, gaiatimes,0)
        gaia_apparentpos=gaia_position.observe(ThingInSpace(smallbody_gaia)).apparent(deflectors=())

This does work to get the correct apparent positions with aberration, but not gravitational deflection. That is because there is no Ephemeris object associated with the list of Gaia positions. While this is not important for my application, it will be if I need high-precision positions at a later date.

How do I enable the list of deflectors for this constructed Gaia position, and in general what is the best way to do so? Is there anything I should keep in mind (I'm not a skyfield expert)?

[brandon-rhodes]

Skyfield's usual technique for supporting deflection is to secretly stash the ephemeris that was used to generate a position. It's saved as the ._ephemeris attribute of the Barycentric object. From there, the Barycentric.observe() method saves the ephemeris to the ._ephemeris attribute of the Astrometric position that it returns. Finally, when Astrometric.apparent() is called, it can use the secretly stashed ephemeris object to find out where the Sun and Jupiter and Saturn are.

In your case, ._ephemeris is still its default value of None, so, as you have discovered, deflection doesn't work.

Try loading an ephemeris and telling your barycentric position about it:

eph = load('de421.bsp')
gaia_position._ephemeris = eph

Hopefully that will get deflection working for you!

I wonder if Skyfield should somehow improve how this works someday, for users in your situation. Should .ephemeris become a public attribute, for example? I'm sure I probably should have structured the classes a bit differently here, in retrospect; but I'm not going to be able to just rework everything at this point, so the solution will doubtless remain a bit clunky.

Let me know if this solutions works for you!

[thunkii]

Thanks, that works!