skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Discrepancy between TEME and ITRF velocity vector magnitudes #139

Closed leblancfg closed 7 years ago

leblancfg commented 7 years ago

Not sure if it's the TEME_to_ITRF coordinate conversion method that's causing the issue, but I get different magnitudes (~1%) for the vTEME and vITRF velocity vectors -- but not for the altitude!

Will need to investigate that method further, but wanted to bring it up first. Seems related to the fact that the method takes units per day, not per second.

from math import sqrt, degrees
import requests
import skyfield.api as sf

from skyfield.constants import ERAD, AU_KM, DAY_S
from skyfield.sgp4lib import TEME_to_ITRF

# Download TLE and get `satellite` instance for RADARSAT-2
tle_url = 'http://www.celestrak.com/cgi-bin/TLE.pl?CATNR=32382'
html = requests.get(tle_url).text
tle = [line for line in html.split('\n')[13:16]]
satellite = sf.EarthSatellite(tle[1], tle[2], tle[0])

now = ts.now()

# First, get the TEME coordinates
rTEME, vTEME, error = satellite._position_and_velocity_TEME_km(now)
altitudeTEME = sqrt(sum(n**2 for n in rTEME))
speedITEME = sqrt(sum(n**2 for n in vTEME))

# Now convert them to ITRF coordinates
vTEME *= DAY_S  # The velocity should be provided in units per day, not per second
rITRF, vITRF = TEME_to_ITRF(now.ut1, rTEME, vTEME)
vITRF /= DAY_S  # Re-transform velocity back to km/s
altitudeITRF = sqrt(sum(n**2 for n in rITRF))
speedITRF = sqrt(sum(n**2 for n in vITRF))

print('TEME Altitude: ', altitudeTEME, 'km from Center of Earth')
print('ITRF Altitude: ', altitudeITRF, 'km from Center of Earth')
print('---')
print('TEME Speed:    ', speedTEME, 'km/s')
print('ITRF Speed:    ', speedITRF, 'km/s')

Returns

TEME Altitude:  7177.3443534706275 km from Center of Earth
ITRF Altitude:  7177.3443534706275 km from Center of Earth
---
TEME Speed:     7.4535983257758955 km/s
ITRF Speed:     7.531937780480272 km/s
brandon-rhodes commented 7 years ago

Unless I am mistaken, the TEME and ITRF reference frames do not rotate at exactly the same rate, so a velocity expressed with respect to them will indeed usually have a different magnitude in each.

Could you double-check me on this? Open the official PDF on the SGP4 algorithm:

http://ww.celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf

And compute the magnitude of the vectors vTEME and vITRF there in the official example from the publication. I get different magnitudes:

In [7]: np.linalg.norm([-4.746131487, 0.785818041, 5.531931288])
Out[7]: 7.3311348276107466
In [8]: np.linalg.norm([-3.225636520, -2.872451450, +5.531924446])
Out[8]: 7.0183969941872917

Do you get the same result?

leblancfg commented 7 years ago

That makes sense, and I confirm I get the same numbers.

Thank you very much for looking into it; I never would have figured that out by myself. I had dismissed the Earth's rotation (0.45 km/s at the equator vs. the ~0.08 km/s discrepancy) as the source of "error", and still don't quite understand how that number shows up. But it does go to show how little I know about coordinate systems and their conversions :expressionless: