Fix: refract() on array should return same angles

If the refract() routine was given an array of altitudes, and some of
them converged to a result more quickly than the rest, then it would
over-refine those early results as it kept looping.  Now, it uses a mask
that lets each result stop converging as soon as that result is ready.

This tightens agreement with NOVAS to 0.00001 arcseconds even for
positions that include the effects of refraction.

Author: brandon-rhodes

builders/build_novas_tests.py

@@ -40,7 +40,7 @@ def main():
         from skyfield.api import Topos, load
         from skyfield.constants import AU_KM, AU_M
         from skyfield.data import hipparcos
-        from skyfield.functions import BytesIO, length_of
+        from skyfield.functions import A, BytesIO, length_of
         from .fixes import low_precision_ERA
 
         OLD_AU_KM = 149597870.691
@@ -237,15 +237,35 @@ def test_refraction{i}():
                 compare(r, {r!r}, 1e-9 * arcsecond)
             """)
 
+    location = novas.make_on_surface(0.0, 0.0, 0, 10.0, 1010.0)
+    angle_in = [-5, -1, 15, 89.95]
+    angle_out = [novas.refract(location, 90 - a, 2) for a in angle_in]
+
+    output(locals(), """\
+        def test_refraction_with_angle_array():
+            r = earthlib.refraction(A{angle_in}, 10.0, 1010.0)
+            compare(r, {angle_out!r}, 1e-9 * arcsecond)
+        """)
+
     northpole = novas.make_on_surface(90.0, 0.0, 0.0, 10.0, 1010.0)
 
-    for i, angle in enumerate([-90, -2, -1, 0, 1, 3, 9, 90]):
-        alt, az = altaz_maneuver(T0, northpole, 0.0, angle, ref=2)
+    alt_in = []
+    alt_out = []
+    for i, alt in enumerate([-90, -2, -1, 0, 1, 3, 9, 90]):
+        alt2, az = altaz_maneuver(T0, northpole, 0.0, alt, ref=2)
         output(locals(), """\
             def test_refract{i}():
-                alt = earthlib.refract({angle!r}, 10.0, 1010.0)
-                compare(alt, {alt!r}, 1e-9 * arcsecond)
+                alt = earthlib.refract({alt!r}, 10.0, 1010.0)
+                compare(alt, {alt2!r}, 1e-9 * arcsecond)
             """)
+        alt_in.append(alt)
+        alt_out.append(alt2)
+
+    output(locals(), """\
+        def test_refract_with_angle_array():
+            r = earthlib.refract(A{alt_in}, 10.0, 1010.0)
+            compare(r, {alt_out!r}, 1e-9 * arcsecond)
+        """)
 
     usno = novas.make_on_surface(38.9215, -77.0669, 92.0, 10.0, 1010.0)
 
@@ -408,12 +428,12 @@ def test_{slug}_topocentric_date{i}(de405):
             compare(az.degrees, {az!r}, 0.00001 * arcsecond)
 
             alt, az, distance = apparent.altaz('standard')
-            compare(alt.degrees, {alt2!r}, 0.0005 * arcsecond)
-            compare(az.degrees, {az2!r}, 0.0005 * arcsecond)
+            compare(alt.degrees, {alt2!r}, 0.00001 * arcsecond)
+            compare(az.degrees, {az2!r}, 0.00001 * arcsecond)
 
             alt, az, distance = apparent.altaz(10.0, 1010.0)
-            compare(alt.degrees, {alt3!r}, 0.0005 * arcsecond)
-            compare(az.degrees, {az3!r}, 0.0005 * arcsecond)
+            compare(alt.degrees, {alt3!r}, 0.00001 * arcsecond)
+            compare(az.degrees, {az3!r}, 0.00001 * arcsecond)
 
         """)
 

documentation/index.rst

@@ -16,7 +16,7 @@ and satellites in orbit around the Earth.
 Its results should agree
 with the positions generated by the United States Naval Observatory
 and their *Astronomical Almanac*
-to within 0.0005 arcseconds (half a “mas” or milliarcsecond).
+to within 0.00001 arcseconds (one-hundredth of a “mas” or “milliarcsecond”).
 
 * Written in pure Python.
 * Installs without any compilation.

skyfield/earthlib.py

@@ -1,7 +1,9 @@
 """Formulae for specific earth behaviors and effects."""
 
-from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos, minimum,
-                   nan_to_num, pi, sin, sqrt, tan, where, zeros_like)
+from numpy import (
+    abs, arcsin, arccos, arctan2, array, bool, clip, cos, minimum,
+    nan_to_num, ones_like, pi, sin, sqrt, tan, where, zeros_like,
+)
 
 from .constants import (AU_M, ANGVEL, DAY_S, DEG2RAD, ERAD,
                         IERS_2010_INVERSE_EARTH_FLATTENING, RAD2DEG, T0, tau)
@@ -148,10 +150,10 @@ def refraction(alt_degrees, temperature_C, pressure_mbar):
 def refract(alt_degrees, temperature_C, pressure_mbar):
     """Given an unrefracted `alt` determine where it will appear in the sky."""
     alt = alt_degrees
-    while True:
-        alt1 = alt
-        alt = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
-        converged = nan_to_num(abs(alt - alt1)).max() <= 3.0e-5
-        if converged:
-            break
+    mask = ones_like(alt, dtype=bool)
+    while mask.sum():
+        alt_previous = alt
+        alt_new = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
+        alt = mask * alt_new + ~mask * alt_previous
+        mask = nan_to_num(abs(alt - alt_previous)) > 3.0e-5
     return alt