Implement basic Mars magnitude function for #210

builders/build_magnitude_tests.py

@@ -32,6 +32,8 @@ def main(argv):
                 args = fields[4], fields[6], fields[8]
             elif planet == 'venus':
                 args = fields[4], fields[6], fields[10]
+            elif planet == 'mars':
+                args = fields[10], fields[12], fields[16]
             elif planet == 'uranus':
                 args = fields[8], fields[10], fields[12], fields[7], fields[5]
             else:
@@ -40,13 +42,20 @@ def main(argv):
             tests[planet].append((args, answer))
 
     for planet, test_list in sorted(tests.items()):
+        tolerance = (
+            # Mars rotation effects are not yet written up.
+            '0.1' if planet == 'mars'
+            # Other planets should match to high precision.
+            else '0.0005'
+        )
+
         if not test_list:
             continue
         print('\ndef test_{}_magnitude_function():'.format(planet))
         for args, answer in test_list:
             joined = ', '.join(str(arg) for arg in args)
             print(f'    mag = m._{planet}_magnitude({joined})')
-            print(f'    assert abs({answer} - mag) < 0.0005')
+            print(f'    assert abs({answer} - mag) < {tolerance}')
 
         print()
         print('    args = [')
@@ -57,7 +66,7 @@ def main(argv):
         print(f'    magnitudes = m._{planet}_magnitude(*args)')
         joined = ', '.join(answer for args, answer in test_list)
         print(f'    expected = [{joined}]')
-        print('    assert all(magnitudes - expected < 0.0005)')
+        print(f'    assert all(magnitudes - expected < {tolerance})')
 
 if __name__ == '__main__':
     main(sys.argv[1:])

skyfield/magnitudelib.py

@@ -92,10 +92,50 @@ def _venus_magnitude(r, delta, ph_ang):
     return -4.384 + distance_mag_factor + ph_ang_factor
 
 def _earth_magnitude(r, delta, ph_ang):
-    distance_mag_factor = 5 * log10 (r * delta)
+    distance_mag_factor = 5 * log10(r * delta)
     ph_ang_factor = -1.060e-03 * ph_ang + 2.054e-04 * ph_ang**2
     return -3.99 + distance_mag_factor + ph_ang_factor
 
+def _mars_magnitude(r, delta, ph_ang):
+    r_mag_factor = 2.5 * log10(r * r)
+    delta_mag_factor = 2.5 * log10(delta * delta)
+    distance_mag_factor = r_mag_factor + delta_mag_factor
+
+    geocentric_phase_angle_limit = 50.0
+
+    condition = ph_ang <= geocentric_phase_angle_limit
+    a = where(condition, 2.267E-02, - 0.02573)
+    b = where(condition, - 1.302E-04, 0.0003445)
+    ph_ang_factor = a * ph_ang + b * ph_ang**2
+
+    # Compute the effective central meridian longitude
+    # eff_CM = ( sub_earth_long + sub_sun_long ) / 2.
+    # if ( abs ( sub_earth_long - sub_sun_long ) > 180. ):
+    #     Eff_CM = Eff_CM + 180.
+    # if ( Eff_CM > 360. ):
+    #     Eff_CM = Eff_CM - 360.
+
+    # ! Use Stirling interpolation to determine the magnitude correction
+    # call Mars_Stirling ( 'R', eff_CM, mag_corr_rot )
+
+    # Convert the ecliptic longitude to Ls
+    # Ls = h_ecl_long + Ls_offset
+    # if ( Ls > 360. ) Ls = Ls - 360.
+    # if ( Ls <   0. ) Ls = Ls + 360.
+
+    # Use Stirling interpolation to determine the magnitude correction
+    # call Mars_Stirling ( 'O', Ls, mag_corr_orb )
+
+    # Until effects from Mars rotation are written up:
+    mag_corr_rot = 0.0
+    mag_corr_orb = 0.0
+
+    # Add factors to determine the apparent magnitude
+    ap_mag = where(ph_ang <= geocentric_phase_angle_limit, -1.601, -0.367)
+    ap_mag += distance_mag_factor + ph_ang_factor + mag_corr_rot + mag_corr_orb
+
+    return ap_mag
+
 def _jupiter_magnitude(r, delta, ph_ang):
     distance_mag_factor = 5 * log10(r * delta)
     geocentric_phase_angle_limit = 12.0

skyfield/tests/test_magnitudes_raw.py

@@ -36,6 +36,23 @@ def test_jupiter_magnitude_function():
     expected = [-1.667, -2.934, 0.790]
     assert all(magnitudes - expected < 0.0005)
 
+def test_mars_magnitude_function():
+    mag = m._mars_magnitude(1.381191244505, 0.37274381097911, 4.8948)
+    assert abs(-2.862 - mag) < 0.1
+    mag = m._mars_magnitude(1.664150453905, 2.58995164518460, 11.5877)
+    assert abs(1.788 - mag) < 0.1
+    mag = m._mars_magnitude(1.591952180003, 3.85882552272013, 167.9000)
+    assert abs(8.977 - mag) < 0.1
+
+    args = [
+        A[1.381191244505, 1.664150453905, 1.591952180003],
+        A[0.37274381097911, 2.58995164518460, 3.85882552272013],
+        A[4.8948, 11.5877, 167.9000],
+    ]
+    magnitudes = m._mars_magnitude(*args)
+    expected = [-2.862, 1.788, 8.977]
+    assert all(magnitudes - expected < 0.1)
+
 def test_mercury_magnitude_function():
     mag = m._mercury_magnitude(0.310295423552, 1.32182643625754, 1.1677)
     assert abs(-2.477 - mag) < 0.0005