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

Make the Earth rotation angle vastly more precise

This commit is contained in:
Brandon Rhodes 2020-07-20 12:01:25 -04:00
parent 0099c353e7
commit d2e0549601
7 changed files with 55 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
builders/build_novas_tests.py
View File

@ -41,6 +41,7 @@ def main():
from skyfield.constants import AU_KM, AU_M
from skyfield.data import hipparcos
from skyfield.functions import length_of
from .fixes import low_precision_ERA
OLD_AU_KM = 149597870.691 # TODO: load from de405
OLD_AU = AU_KM / OLD_AU_KM
@ -225,7 +226,8 @@ def output_subroutine_tests(dates):
def test_sidereal_time_with_nonzero_delta_t_on_date{i}():
jd = load.timescale(delta_t=99.9).tt_jd({jd!r} + 99.9 * one_second)
compare(earthlib.sidereal_time(jd), {result2!r}, 1e-13)
with low_precision_ERA():
compare(earthlib.sidereal_time(jd), {result2!r}, 1e-13)
""")
p, v = novas.starvectors(novas.make_cat_entry(
@ -291,7 +293,8 @@ def output_subroutine_tests(dates):
output(locals(), """\
def test_ITRF_to_GCRS_conversion_on_date{i}():
jd = load.timescale(delta_t={delta_t!r}).tt_jd({tt!r})
position = positionlib.ITRF_to_GCRS(jd, {vector!r})
with low_precision_ERA():
position = positionlib.ITRF_to_GCRS(jd, {vector!r})
compare(position, {result!r}, 1e-13)
""")

6
skyfield/documentation/earth-satellites.rst
View File

@ -292,7 +292,7 @@ than those of the old J2000 system.)
.. testoutput::
[-3918.87651974 -1887.64832658 5209.08802573]
[-3918.8765203 -1887.6483254 5209.08802574]
.. would love to be able to do this someday - see the SPICE source file
nearpt.f
@ -365,7 +365,7 @@ just as you did for the position measured from the Earth’s center:
.. testoutput::
[ 331.6188371 392.18468546 1049.76011302]
[ 331.61883722 392.18468536 1049.76011302]
But the most popular approach is to ask the topocentric position
for its altitude and azimuth coordinates,
@ -665,7 +665,7 @@ that are limiting this TLE set’s predictions:
mean eccentricity is outside the range 0.0 to 1.0
After:
[-5021.82658136 742.71506482 3831.57403957]
[-5021.82658191 742.71506112 3831.57403957]
mrt is less than 1.0 which indicates the satellite has decayed
If you use a ``Time`` array to ask about an entire range of dates,

9
skyfield/earthlib.py
View File

@ -133,7 +133,7 @@ def sidereal_time(t):
# Compute the Earth Rotation Angle. Time argument is UT1.
theta = earth_rotation_angle(t.ut1)
theta = earth_rotation_angle(t.whole, t.ut1_fraction)
# The equinox method. See Circular 179, Section 2.6.2.
# Precession-in-RA terms in mean sidereal time taken from third
@ -152,16 +152,15 @@ def sidereal_time(t):
return (st / 54000.0 + theta * 24.0) % 24.0
def earth_rotation_angle(jd_ut1):
def earth_rotation_angle(jd_ut1, fraction_ut1=0.0):
"""Return the value of the Earth Rotation Angle (theta) for a UT1 date.
Uses the expression from the note to IAU Resolution B1.8 of 2000.
Returns a fraction between 0.0 and 1.0 whole rotations.
"""
thet1 = 0.7790572732640 + 0.00273781191135448 * (jd_ut1 - T0)
thet3 = jd_ut1 % 1.0
return (thet1 + thet3) % 1.0
th = 0.7790572732640 + 0.00273781191135448 * (jd_ut1 - T0 + fraction_ut1)
return (th % 1.0 + jd_ut1 % 1.0 + fraction_ut1) % 1.0
def refraction(alt_degrees, temperature_C, pressure_mbar):

2
skyfield/positionlib.py
View File

@ -743,7 +743,7 @@ def ITRF_to_GCRS(t, rITRF):
# Todo: wobble
spin = rot_z(t.gast * tau / 24.0)
spin = rot_z(t.gast / 24.0 * tau)
position = mxv(spin, array(rITRF))
return mxv(t.MT, position)

22
skyfield/tests/fixes.py
View File

@ -1,6 +1,7 @@
"""Helpers for making Skyfield tests stable."""
import datetime as dt
from skyfield import earthlib
import skyfield.api
import skyfield.timelib
@ -27,3 +28,24 @@ def teardown():
skyfield.api.load = skyfield.iokit.Loader('.')
skyfield.timelib.Timescale._utcnow = dt.datetime.utcnow
dt.datetime = _real_datetime_class
class low_precision_ERA(object):
"""Compute the Earth rotation angle with only a single float for UT1.
Skyfield now uses two floats ``t.whole`` and ``t.ut1_fraction`` to
represent the UT1 Julian date, supplying an additional 16 digits of
precision. For the Earth rotation angle, which moves very quickly
per unit time compared to most other astronomical quantities, this
knocks Skyfield out of agreement with other libraries like NOVAS and
SOFA that round UT1 to a single float. Various tests use this
context manager to make Skyfield match the lower-precision output.
"""
def __enter__(self):
self.saved = earthlib.earth_rotation_angle
earthlib.earth_rotation_angle = self.era
def __exit__(self, *args):
earthlib.earth_rotation_angle = self.saved
def era(self, whole, fraction):
rounded_single_float = whole + fraction
return self.saved(rounded_single_float)

25
skyfield/tests/test_against_novas.py
View File

@ -7,6 +7,7 @@ from skyfield.api import Topos, load
from skyfield.constants import AU_KM, AU_M
from skyfield.data import hipparcos
from skyfield.functions import length_of
from .fixes import low_precision_ERA
OLD_AU_KM = 149597870.691 # TODO: load from de405
OLD_AU = AU_KM / OLD_AU_KM
@ -246,7 +247,8 @@ def test_sidereal_time_on_date0():
def test_sidereal_time_with_nonzero_delta_t_on_date0():
jd = load.timescale(delta_t=99.9).tt_jd(2440423.345833333 + 99.9 * one_second)
compare(earthlib.sidereal_time(jd), 16.195436229760602, 1e-13)
with low_precision_ERA():
compare(earthlib.sidereal_time(jd), 16.195436229760602, 1e-13)
def test_sidereal_time_on_date1():
jd = load.timescale(delta_t=0.0).tt_jd(2448031.5)
@ -254,7 +256,8 @@ def test_sidereal_time_on_date1():
def test_sidereal_time_with_nonzero_delta_t_on_date1():
jd = load.timescale(delta_t=99.9).tt_jd(2448031.5 + 99.9 * one_second)
compare(earthlib.sidereal_time(jd), 15.825907462991848, 1e-13)
with low_precision_ERA():
compare(earthlib.sidereal_time(jd), 15.825907462991848, 1e-13)
def test_sidereal_time_on_date2():
jd = load.timescale(delta_t=0.0).tt_jd(2451545.0)
@ -262,7 +265,8 @@ def test_sidereal_time_on_date2():
def test_sidereal_time_with_nonzero_delta_t_on_date2():
jd = load.timescale(delta_t=99.9).tt_jd(2451545.0 + 99.9 * one_second)
compare(earthlib.sidereal_time(jd), 18.69737482966941, 1e-13)
with low_precision_ERA():
compare(earthlib.sidereal_time(jd), 18.69737482966941, 1e-13)
def test_sidereal_time_on_date3():
jd = load.timescale(delta_t=0.0).tt_jd(2456164.5)
@ -270,7 +274,8 @@ def test_sidereal_time_on_date3():
def test_sidereal_time_with_nonzero_delta_t_on_date3():
jd = load.timescale(delta_t=99.9).tt_jd(2456164.5 + 99.9 * one_second)
compare(earthlib.sidereal_time(jd), 22.2439084998698, 1e-13)
with low_precision_ERA():
compare(earthlib.sidereal_time(jd), 22.2439084998698, 1e-13)
def test_star_vector():
star = starlib.Star(ra_hours=2.530301028, dec_degrees=89.264109444,
@ -455,22 +460,26 @@ def test_from_altaz_7(earth):
def test_ITRF_to_GCRS_conversion_on_date0():
jd = load.timescale(delta_t=39.707).tt_jd(2440423.345833333)
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
with low_precision_ERA():
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
compare(position, (0.5701172053658128, -1.5232987806096392, 1.3017400651201707), 1e-13)
def test_ITRF_to_GCRS_conversion_on_date1():
jd = load.timescale(delta_t=57.1136).tt_jd(2448031.5)
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
with low_precision_ERA():
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
compare(position, (0.41362649279562963, -1.5741081933652488, 1.3004216700893525), 1e-13)
def test_ITRF_to_GCRS_conversion_on_date2():
jd = load.timescale(delta_t=63.8285).tt_jd(2451545.0)
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
with low_precision_ERA():
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
compare(position, (1.3757008573963405, -0.8702954291925735, 1.3000126987400913), 1e-13)
def test_ITRF_to_GCRS_conversion_on_date3():
jd = load.timescale(delta_t=66.7846).tt_jd(2456164.5)
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
with low_precision_ERA():
position = positionlib.ITRF_to_GCRS(jd, [1.1, 1.2, 1.3])
compare(position, (1.5243574049688486, 0.5755748855663746, 1.2980940077752074), 1e-13)
def test_tdb_minus_tt_on_date0():

4
skyfield/tests/test_positions.py
View File

@ -5,6 +5,7 @@ from skyfield.earthlib import earth_rotation_angle
from skyfield.functions import length_of
from skyfield.positionlib import ICRF, ITRF_to_GCRS2, _GIGAPARSEC_AU
from skyfield.starlib import Star
from .fixes import low_precision_ERA
def test_subtraction():
p0 = ICRF((10,20,30), (40,50,60))
@ -201,7 +202,8 @@ def test_cirs_sofa():
for ((ra_icrs, dec_icrs, tdb), (ra_sofa, dec_sofa)) in zip(test_data, sofa_results):
ss = Star(ra_hours=(ra_icrs / 15.0), dec_degrees=dec_icrs)
st = ts.tdb(jd=tdb)
ra_cirs, dec_cirs, _ = earth.at(st).observe(ss).apparent().cirs_radec(st)
with low_precision_ERA():
ra_cirs, dec_cirs, _ = earth.at(st).observe(ss).apparent().cirs_radec(st)
assert np.allclose(ra_cirs._degrees, ra_sofa, rtol=0.0, atol=tol)
assert np.allclose(dec_cirs._degrees, dec_sofa, rtol=0.0, atol=tol)