Make all magnitude routines support vectors (#210)

CHANGELOG.rst

@@ -5,6 +5,13 @@ Changelog
 .. TODO After finding how to test TIRS reference frame, add it to changelog.
         And double-check the constellation boundaries array.
 
+v1.39 — ?
+---------
+
+* The prototype :func:`~skyfield.magnitudelib.planetary_magnitude()`
+  function now works not only when given a single position, but when
+  given a vector of several positions.
+
 v1.38 — 2021 April 3
 --------------------
 

builders/build_magnitude_tests.py

@@ -39,9 +39,6 @@ def main(argv):
             answer = fields[-1]
             tests[planet].append((args, answer))
 
-    # from pprint import pprint
-    # pprint(tests)
-
     for planet, test_list in sorted(tests.items()):
         if not test_list:
             continue
@@ -51,9 +48,6 @@ def main(argv):
             print(f'    mag = m._{planet}_magnitude({joined})')
             print(f'    assert abs({answer} - mag) < 0.0005')
 
-        if planet != 'venus':
-            continue
-
         print()
         print('    args = [')
         for arg_vector in zip(*[args for args, answer in test_list]):

skyfield/magnitudelib.py

@@ -99,23 +99,20 @@ def _earth_magnitude(r, delta, ph_ang):
 def _jupiter_magnitude(r, delta, ph_ang):
     distance_mag_factor = 5 * log10(r * delta)
     geocentric_phase_angle_limit = 12.0
-
-    if ph_ang <= geocentric_phase_angle_limit:  # TODO: time arrays
-        ph_ang_factor = -3.7E-04 * ph_ang + 6.16E-04 * ph_ang**2
-    else:
-        ph_ang_factor = -2.5 * log10(
-            1.0 - 1.507 * (ph_ang / 180.)
-            - 0.363 * (ph_ang / 180.)**2
-            - 0.062 * (ph_ang / 180.)**3
-            + 2.809 * (ph_ang / 180.)**4
-            - 1.876 * (ph_ang / 180.)**5
-        )
-
-    if ph_ang <= geocentric_phase_angle_limit:
-        ap_mag = -9.395 + distance_mag_factor + ph_ang_factor
-    else:
-        ap_mag = -9.428 + distance_mag_factor + ph_ang_factor
-
+    ph_ang_pi = ph_ang / 180.0
+    ph_ang_factor = where(
+        ph_ang <= geocentric_phase_angle_limit,
+        (6.16E-04 * ph_ang - 3.7E-04) * ph_ang,
+        -2.5 * log10(
+            ((((- 1.876 * ph_ang_pi + 2.809) * ph_ang_pi - 0.062) * ph_ang_pi
+              - 0.363) * ph_ang_pi - 1.507) * ph_ang_pi + 1.0
+        ),
+    )
+    ap_mag = where(
+        ph_ang <= geocentric_phase_angle_limit,
+        -9.395 + distance_mag_factor + ph_ang_factor,
+        -9.428 + distance_mag_factor + ph_ang_factor,
+    )
     return ap_mag
 
 def _uranus_magnitude(r, delta, ph_ang,
@@ -126,8 +123,11 @@ def _uranus_magnitude(r, delta, ph_ang,
     sub_lat_factor = -0.00084 * sub_lat_planetog
     geocentric_phase_angle_limit = 3.1
     ap_mag = -7.110 + distance_mag_factor + sub_lat_factor
-    if ph_ang > geocentric_phase_angle_limit:  # TODO: time arrays
-        ap_mag += 6.587e-3 * ph_ang + 1.045e-4 * ph_ang**2
+    ap_mag += where(
+        ph_ang > geocentric_phase_angle_limit,
+        (1.045e-4 * ph_ang + 6.587e-3) * ph_ang,
+        0.0,
+    )
     return ap_mag
 
 _FUNCTIONS = {

skyfield/tests/test_magnitudes_raw.py

@@ -10,6 +10,15 @@ def test_earth_magnitude_function():
     mag = m._earth_magnitude(0.983356467727, 0.62933287342927, 175.6869)
     assert abs(1.122 - mag) < 0.0005
 
+    args = [
+        A[0.983331936476, 0.983356079811, 0.983356467727],
+        A[1.41317594650699, 0.26526856764726, 0.62933287342927],
+        A[8.7897, 4.1369, 175.6869],
+    ]
+    magnitudes = m._earth_magnitude(*args)
+    expected = [-3.269, -6.909, 1.122]
+    assert all(magnitudes - expected < 0.0005)
+
 def test_jupiter_magnitude_function():
     mag = m._jupiter_magnitude(5.446231815414, 6.44985867459088, 0.2446)
     assert abs(-1.667 - mag) < 0.0005
@@ -18,6 +27,15 @@ def test_jupiter_magnitude_function():
     mag = m._jupiter_magnitude(5.227587855371, 5.23501920009381, 147.0989)
     assert abs(0.790 - mag) < 0.0005
 
+    args = [
+        A[5.446231815414, 4.957681473205, 5.227587855371],
+        A[6.44985867459088, 3.95393078136013, 5.23501920009381],
+        A[0.2446, 0.3431, 147.0989],
+    ]
+    magnitudes = m._jupiter_magnitude(*args)
+    expected = [-1.667, -2.934, 0.790]
+    assert all(magnitudes - expected < 0.0005)
+
 def test_mercury_magnitude_function():
     mag = m._mercury_magnitude(0.310295423552, 1.32182643625754, 1.1677)
     assert abs(-2.477 - mag) < 0.0005
@@ -26,6 +44,15 @@ def test_mercury_magnitude_function():
     mag = m._mercury_magnitude(0.448947624811, 0.56004973217883, 178.7284)
     assert abs(7.167 - mag) < 0.0005
 
+    args = [
+        A[0.310295423552, 0.413629222334, 0.448947624811],
+        A[1.32182643625754, 0.92644808718613, 0.56004973217883],
+        A[1.1677, 90.1662, 178.7284],
+    ]
+    magnitudes = m._mercury_magnitude(*args)
+    expected = [-2.477, 0.181, 7.167]
+    assert all(magnitudes - expected < 0.0005)
+
 def test_uranus_magnitude_function():
     mag = m._uranus_magnitude(18.321003215845, 17.3229728525108, 0.0410, -20.29, -20.28)
     assert abs(5.381 - mag) < 0.0005
@@ -34,6 +61,17 @@ def test_uranus_magnitude_function():
     mag = m._uranus_magnitude(19.38003071775, 11.1884243801383, 161.7728, -71.16, 55.11)
     assert abs(8.318 - mag) < 0.0005
 
+    args = [
+        A[18.321003215845, 20.096361095266, 19.38003071775],
+        A[17.3229728525108, 21.0888470145276, 11.1884243801383],
+        A[0.0410, 0.0568, 161.7728],
+        A[-20.29, 1.02, -71.16],
+        A[-20.28, 0.97, 55.11],
+    ]
+    magnitudes = m._uranus_magnitude(*args)
+    expected = [5.381, 6.025, 8.318]
+    assert all(magnitudes - expected < 0.0005)
+
 def test_venus_magnitude_function():
     mag = m._venus_magnitude(0.722722540169, 1.71607489554051, 1.3232)
     assert abs(-3.917 - mag) < 0.0005