skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

Vector sum for satellites does not work #928

Closed lpmrfentazis closed 5 months ago

lpmrfentazis commented 7 months ago

Hi, I was looking for an alternative to pyorbital and came across skyfield. But as soon as I decided to calculate the satellite position I ran into a problem:

Following https://rhodesmill.org/skyfield/earth-satellites.html difference = satellite - bluffton image

I get ValueError: you can only subtract two vectors if they both start at the same center, by experimentation I found that I can avoid the problem by doing .at for both satellite and observer. But in this case the calculations are not correct, and suspiciously strongly resemble the coordinates of the sun. distance is 1 au

image

What am I doing wrong?

image

full code as text: https://pastebin.com/CxqdQ85T

Python 3.11.4 skyfield 1.46

MichaelBonnet commented 6 months ago

You stumbled upon the root cause:

I get ValueError: you can only subtract two vectors if they both start at the same center, by experimentation I found that I can avoid the problem by doing .at for both satellite and observer.

The issue is that your observer's coordinate system is centered on the solar system barycenter. So, your position vector elements are seemingly colossal:

print(observer.at(t).position.km)

[3.82582172e+07 1.30718801e+08 5.66792475e+07]

While your satellite position vector elements are orders of magnitude smaller, as they are centered on Earth and refer to a satellite orbiting it (I used the ISS here for illustration):

print(sats["ISS (ZARYA)"].at(t).position.km)

[ 378.7947757 -1244.90480596 1311.47538942]

Without transformations between the two reference frames, you're going to get wildly inaccurate results.


My edited version of your code:

from skyfield.api import load, wgs84, N, E

stations_url = "http://celestrak.org/NORAD/elements/stations.txt"
satellites = load.tle_file(stations_url)
sats = {sat.name: sat for sat in satellites}
ISS = sats['ISS (ZARYA)']

objects = load('de421.bsp')
ts = load.timescale()
t = ts.utc(2030, 12, 7, 10, 0, 0)

earth, sun = objects['earth'], objects['sun']

observer = earth + wgs84.latlon(55.6713361 * N, 37.6253844 * E, 180)

astrometric = observer.at(t).observe(sun)
el, az, d = astrometric.apparent().altaz()

print(f"Observer-Sun Elevation : {el.degrees}")
print(f"Observer-Sun Azimuth   : {az.degrees}")
print(f"Observer-Sun Distance  : {d.km}")
print(f"Sat Position (km)      : {ISS.at(t).position.km}")
print(f"Observer Position (km) : {observer.at(t).position.km}")

print("=======================================================================")
difference = ISS.at(t) - observer.at(t)
el, az, d = difference.altaz()

print(el.degrees)
print(az.degrees)
print(d.km)
brandon-rhodes commented 5 months ago

Thank you for pointing out that problem, @lpmrfentazis, and for pointing out the correct workaround, @MichaelBonnet! I have just landed a commit (see the link just above this comment) that should protect users in the future from subtracting the wrong vectors. The same check was already in place for vector functions; I'm not sure why it was never added to the position class.