entrouvert
/
debian-skyfield
17
0
Fork 0
Code Issues Pull Requests Releases Activity

Make TEME→ITRF test 20,000 times more sensitive 

It turns out I was somehow being hamfisted in how I built the Julian
date to pass to `TEME_to_ITRF()` in the test.  By building the Julian
date manually, without going through a `Time` object, I can produce far
more exactly the quantities specified in AIAA-2006-6753 Appendix C.
This commit is contained in:
Brandon Rhodes 2020-07-14 15:44:55 -04:00
parent af622b7f69
commit da2e1fd735
2 changed files with 17 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
skyfield/sgp4lib.py
View File

@ -263,8 +263,6 @@ class EarthSatellite(VectorFunction):
i = jd.argsort()
return ts.tt_jd(jd[i]), v[i]
_second = 1.0 / (24.0 * 60.0 * 60.0)
def theta_GMST1982(jd_ut1):
"""Return the angle of Greenwich Mean Standard Time 1982 given the JD.
@ -278,8 +276,8 @@ def theta_GMST1982(jd_ut1):
t = (jd_ut1 - T0) / 36525.0
g = 67310.54841 + (8640184.812866 + (0.093104 + (-6.2e-6) * t) * t) * t
dg = 8640184.812866 + (0.093104 * 2.0 + (-6.2e-6 * 3.0) * t) * t
theta = (jd_ut1 % 1.0 + g * _second % 1.0) * tau
theta_dot = (1.0 + dg * _second / 36525.0) * tau
theta = (jd_ut1 % 1.0 + g / DAY_S % 1.0) % 1.0 * tau
theta_dot = (1.0 + dg / (DAY_S * 36525.0)) * tau
return theta, theta_dot
_zero_zero_minus_one = array((0.0, 0.0, -1.0))

26
skyfield/tests/test_earth_satellites.py
View File

@ -5,6 +5,7 @@ from skyfield import api
from skyfield.api import EarthSatellite, load
from skyfield.constants import AU_KM, AU_M
from skyfield.sgp4lib import TEME_to_ITRF
from skyfield.timelib import julian_date
iss_tle0 = """\
1 25544U 98067A 18184.80969102 .00001614 00000-0 31745-4 0 9993
@ -63,7 +64,7 @@ from ..constants import DEG2RAD
arcminute = DEG2RAD / 60.0
arcsecond = arcminute / 60.0
second = 1.0 / (24.0 * 60.0 * 60.0)
seconds_per_day = 86400.0
# Note that the following test is based specifically on Revision 2 of
# "Revisiting Spacetrack Report #3" AIAA 2006-6753 (earlier versions of
@ -76,25 +77,28 @@ def test_appendix_c_conversion_from_TEME_to_ITRF():
vTEME = array([-4.746131487, 0.785818041, 5.531931288])
vTEME = vTEME * 24.0 * 60.0 * 60.0 # km/s to km/day
ts = api.load.timescale()
jd_ut1 = ts.tt(2004, 4, 6, 7, 51, 28.386 - 0.439961).tt
jd_utc = julian_date(2004, 4, 6, 7, 51, 28.386)
d_ut1 = -0.439961
jd_ut1 = jd_utc + d_ut1 / 86400.0
xp = -0.140682 * arcsecond
yp = 0.333309 * arcsecond
rITRF, vITRF = TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp, yp)
meter = 1e-3
epsilon = 5e-8 # Why not 1e-8, which would match all of their digits?
assert abs(-1033.47938300 - rITRF[0]) < 0.1 * meter
assert abs(+7901.29527540 - rITRF[1]) < 0.1 * meter
assert abs(+6380.35659580 - rITRF[2]) < 0.1 * meter
assert abs(-1033.47938300 - rITRF[0]) < epsilon
assert abs(+7901.29527540 - rITRF[1]) < epsilon
assert abs(+6380.35659580 - rITRF[2]) < epsilon
vITRF_per_second = vITRF * second
vITRF_per_second = vITRF / seconds_per_day
assert abs(-3.225636520 - vITRF_per_second[0]) < 1e-4 * meter
assert abs(-2.872451450 - vITRF_per_second[1]) < 1e-4 * meter
assert abs(+5.531924446 - vITRF_per_second[2]) < 1e-4 * meter
epsilon = 7e-8 # Why not 1e-9, which would match all of their digits?
assert abs(-3.225636520 - vITRF_per_second[0]) < epsilon
assert abs(-2.872451450 - vITRF_per_second[1]) < epsilon
assert abs(+5.531924446 - vITRF_per_second[2]) < epsilon
def test_appendix_c_satellite():
ts = api.load.timescale(builtin=True)