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)
 
         """)
 

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