skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.41k stars 211 forks source link

improve observe() speed by making light time correction optional #993

Open smroid opened 1 month ago

smroid commented 1 month ago

Hi,

I'm writing code to compute RA/Dec for a bunch of asteroids from MPCORB.DAT and I'm finding that the earth.observe(asteroid) call is my bottleneck.

Investigating, I find that if I nerf _correct_for_light_travel_time() in vectorlib.py, by skipping its for loop and just doing the one target._at(t) call, I triple my performance with only a small change in RA/Dec values.

Proposal: allow an optional argument to observe() to skip the light travel time correction?

EndlessDex commented 4 weeks ago

This is already a feature!

You can find information on it on the earth satellite doc page. Here for how to do it: https://rhodesmill.org/skyfield/earth-satellites.html#satellite-altitude-azimuth-and-distance And here for why its good for satellites and bad for others: https://rhodesmill.org/skyfield/earth-satellites.html#avoid-calling-the-observe-method

Basically any ephemeris (represented as a VectorSum) can be differenced and "at()"-ed to get a rough vector. See example below.

Note that this excludes both the light-travel correction and the light-bending correction that you get from observe and apparent, respectively.

import skyfield.api as sf
from skyfield.vectorlib import VectorSum
from skyfield.positionlib import Apparent, Astrometric, Barycentric, Geocentric

ts: sf.Timescale = sf.load.timescale()
t: sf.Time = ts.utc(2024, 8, 15)
planets = sf.load("de421.bsp")
earth: VectorSum = planets["earth"]
mars: VectorSum = planets["mars"]

barycentric: Barycentric = earth.at(t)
astrometric: Astrometric = barycentric.observe(mars)
apparent: Apparent = astrometric.apparent()
ra, dec, dist = apparent.radec()

print(f"Mars is {dist.au:.9f} au from Earth at ({ra.hours:.9f}, {dec.degrees:.9f})")
# Result = "Mars is 1.526189175 au from Earth at (5.029890987, 22.378234247)"

vector: VectorSum = mars - earth
geocentric: Geocentric = vector.at(t)
ra, dec, dist = geocentric.radec()
print(f"Mars is {dist.au:.9f} au from Earth at ({ra.hours:.9f}, {dec.degrees:.9f})")
# result = "Mars is 1.526279313 au from Earth at (5.030301116, 22.378944314)"
smroid commented 4 weeks ago

Works, and more than triples the speed. For my purposes, the accuracy is fine:

For the asteroid (1) Ceres, I get:

using observe(): (<Angle 18h 34m 26.10s>, <Angle -30deg 54' 56.2">, <Distance 2.13274 au>)

using the vector difference method: (<Angle 18h 34m 26.98s>, <Angle -30deg 54' 57.7">, <Distance 2.13272 au>)

However, without your assistance, I would never have been able to discover this. Can the documentation for observe() have a note pointing us towards the more quick-n-dirty approach?

EndlessDex commented 4 weeks ago

Awesome! Glad it works for you :D

Just make sure you are aware that the further (and/or faster) an object is, the worse the error will be. Anything extra-solar will likely be very wrong.