Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
smroid opened this issue Aug 12, 2024 · 3 comments
Open

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

smroid opened this issue Aug 12, 2024 · 3 comments

Comments

@smroid
Copy link

smroid commented Aug 12, 2024

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
Copy link

EndlessDex commented Aug 16, 2024

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
Copy link
Author

smroid commented Aug 16, 2024

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
Copy link

EndlessDex commented Aug 16, 2024 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants