Coverage for python / lsst / pipe / tasks / ssp / ephem.py: 0%

1from astropy.time import Time 

2import astropy.units as u 

3 

4import jax 

5from jorbit import Particle 

6 

7import numpy as np 

8 

9jax.config.update("jax_enable_x64", True) 

10 

11 

12def _aux_compute_ephemerides(provID, ephTimes, mpcorb): 

13 

14 kms = u.km / u.s 

15 

16 # Ephemerides 

17 (packed, a, e, i, node, argperi, M, epoch, H, G) = mpcorb.query( 

18 "unpacked_primary_provisional_designation == @provID", engine="python" 

19 )["packed_primary_provisional_designation a e i node argperi mean_anomaly epoch_mjd h g".split()].iloc[0] 

20 

21 # FIXME: hack until jorbit is fixed (https://github.com/ben-cassese/jorbit/issues/26) 

22 p = Particle.from_horizons(name=packed, time=Time(epoch, format="mjd", scale="tdb")) 

23 

24 eph, xx, vv, obs = p.ephemeris(times=ephTimes, observer="rubin") 

25 

26 xx, vv, obs = np.array(xx).T, np.array(vv).T, np.array(obs).T 

27 vv = (vv * u.au / u.day).to(kms).value # convert from AU/day to km/s 

28 

29 # rate of motion 

30 dt = 1.0 / (3600.0 + 24.0) * u.s 

31 ephTimes2 = ephTimes + dt 

32 eph2, _, _, _ = p.ephemeris(times=ephTimes2, observer="rubin") 

33 dlon, dlat = eph.spherical_offsets_to(eph2) # small offsets on tangent plane 

34 mu_lon = (dlon / dt).to(u.deg / u.day) # ≈ d(RA*cosDec)/dt 

35 mu_lat = (dlat / dt).to(u.deg / u.day) # ≈ d(Dec)/dt 

36 mu = (eph.separation(eph2) / dt).to(u.deg / u.day) # total rate of motion 

37 

38 return eph, (H, G), xx, vv, obs, mu_lon, mu_lat, mu