[TEST] Add new comprehensive tests for magnetic benchmarks including edge cases and symmetry validation

- Introduced new test file `test_magnetic_benchmark_II.py` covering multiple scenarios:
  - Symmetry checks for TMI profiles.
  - Edge cases for zero and negative susceptibilities.
  - Observation distance decay validation.
  - Kernel property tests for symmetry and resolution.
- Enhanced and streamlined existing tests in `test_magnetic_benchmark.py`:
  - Adjusted formatting for clearer readability.
  - Updated comments and formulas for analytical benchmarks to improve documentation consistency.

Author: Leguark

tests/test_common/test_modules/test_magnetic_benchmark.py

@@ -15,19 +15,19 @@
 # =============================================================================
 
 @pytest.mark.parametrize("inclination,declination,expected_lx,expected_ly,expected_lz", [
-    (0.0, 0.0, 1.0, 0.0, 0.0),      # Horizontal north
-    (90.0, 0.0, 0.0, 0.0, 1.0),     # Vertical down (north pole)
-    (-90.0, 0.0, 0.0, 0.0, -1.0),   # Vertical up (south pole)
-    (0.0, 90.0, 0.0, 1.0, 0.0),     # Horizontal east
-    (45.0, 0.0, 0.707107, 0.0, 0.707107),  # 45° down to north
-    (45.0, 45.0, 0.5, 0.5, 0.707107),      # 45° down, 45° east
+        (0.0, 0.0, 1.0, 0.0, 0.0),  # Horizontal north
+        (90.0, 0.0, 0.0, 0.0, 1.0),  # Vertical down (north pole)
+        (-90.0, 0.0, 0.0, 0.0, -1.0),  # Vertical up (south pole)
+        (0.0, 90.0, 0.0, 1.0, 0.0),  # Horizontal east
+        (45.0, 0.0, 0.707107, 0.0, 0.707107),  # 45° down to north
+        (45.0, 45.0, 0.5, 0.5, 0.707107),  # 45° down, 45° east
 ])
 def test_direction_cosines(inclination, declination, expected_lx, expected_ly, expected_lz):
     """Test direction cosines computation for various field geometries."""
     l = _direction_cosines(inclination, declination)
 
     # Check components match expectations
-    np.testing.assert_allclose(l[0], expected_lx, atol=1e-5, 
+    np.testing.assert_allclose(l[0], expected_lx, atol=1e-5,
                                err_msg=f"l_x mismatch for I={inclination}, D={declination}")
     np.testing.assert_allclose(l[1], expected_ly, atol=1e-5,
                                err_msg=f"l_y mismatch for I={inclination}, D={declination}")
@@ -55,11 +55,11 @@ def test_direction_cosines_wraparound():
 # =============================================================================
 
 @pytest.mark.parametrize("inclination,declination,intensity", [
-    (90.0, 0.0, 50000.0),    # Vertical field (pole)
-    (0.0, 0.0, 45000.0),     # Horizontal field (equator)
-    (60.0, 10.0, 48000.0),   # Mid-latitude field
-    (45.0, 45.0, 50000.0),   # Oblique field
-    (-30.0, -15.0, 35000.0), # Southern hemisphere
+        (90.0, 0.0, 50000.0),  # Vertical field (pole)
+        (0.0, 0.0, 45000.0),  # Horizontal field (equator)
+        (60.0, 10.0, 48000.0),  # Mid-latitude field
+        (45.0, 45.0, 50000.0),  # Oblique field
+        (-30.0, -15.0, 35000.0),  # Southern hemisphere
 ])
 def test_path_equivalence(inclination, declination, intensity):
     """Test that pre-projected and raw V computation paths give identical results."""
@@ -82,8 +82,8 @@ def test_path_equivalence(inclination, declination, intensity):
 
     # Should match to floating point precision
     np.testing.assert_allclose(
-        tmi_path1, 
-        tmi_path2[0], 
+        tmi_path1,
+        tmi_path2[0],
         rtol=1e-10,
         err_msg=f"Paths differ for I={inclination}, D={declination}, F={intensity}"
     )
@@ -103,19 +103,19 @@ def test_path_equivalence(inclination, declination, intensity):
 def test_magnetics_sphere_analytical_benchmark_induced_only(inclination, declination, intensity_nT, rtol):
     """
     Benchmark comparing induced-only TMI against analytical solution for a uniformly
-    susceptible sphere in a vertical inducing field.
+    susceptible sphere with observation points along the vertical axis.
 
-    Geometry mirrors gravity sphere benchmark: observe along vertical line above center.
+    For a sphere with induced magnetization M = χ·B₀/μ₀, the magnetic field is that
+    of a magnetic dipole with moment m = (4/3)πR³·M.
 
-    Analytical (outside sphere, along axis):
-        m = V * M = (4/3)π R^3 * (χ * B0 / μ0)
-        Bz = μ0/(4π) * (2 m) / r^3  (dipole field on axis)
-        ΔT = B · l = Bz (since l = z-hat for vertical field)
-        => ΔT = (2/3) * R^3 * χ * B0 / r^3  [Tesla]
-        Convert to nT by × 1e9 if B0 in Tesla. Here intensity_nT is in nT, so use directly:
-        ΔT[nT] = (2/3) * R^3 * χ * intensity_nT / r^3
-
-    We accept ~15–20% error due to voxelization discretization.
+    The TMI anomaly for observation on the z-axis above the sphere center is:
+        ΔT = (R³·χ·F)/(3r³) · [3·lz² - 1 + 2·(lx² + ly²)]
+    
+    Where l = [lx, ly, lz] are the direction cosines of the IGRF field.
+    
+    Special cases:
+    - Vertical field (I=90°): l=[0,0,1] → ΔT = (2/3)·R³·χ·F/r³  (your current formula)
+    - Horizontal field (I=0°): l=[±1,0,0] or [0,±1,0] → ΔT = -(1/3)·R³·χ·F/r³
     """
 
     # Sphere parameters
@@ -140,7 +140,7 @@ def test_magnetics_sphere_analytical_benchmark_induced_only(inclination, declina
         radius=np.array([200.0, 200.0, 800.0]),
     )
 
-    # Magnetic kernel (pre-projected TMI kernel per voxel recommended by plan)
+    # Magnetic kernel
     igrf_params = {
         "inclination": inclination,
         "declination": declination,
@@ -152,12 +152,10 @@ def test_magnetics_sphere_analytical_benchmark_induced_only(inclination, declina
             geophysics_grid, igrf_params, compute_tmi=True, units_nT=True
         )
     except TypeError:
-        # Some implementations may return the kernel directly rather than a dict; handle both
         mag_kern_out = calculate_magnetic_gradient_tensor(
             geophysics_grid, igrf_params
         )
 
-    # Support both dict output with key 'tmi_kernel' or direct array output
     if isinstance(mag_kern_out, dict):
         tmi_kernel = mag_kern_out.get("tmi_kernel", None)
         if tmi_kernel is None:
@@ -175,39 +173,70 @@ def test_magnetics_sphere_analytical_benchmark_induced_only(inclination, declina
         vc = voxel_centers[sl]
         inside = (np.linalg.norm(vc - center, axis=1) <= R).astype(float)
         chi_vox = inside * chi
-        # Forward model: sum(chi * tmi_kernel)
         numerical_tmi.append(np.sum(chi_vox * tmi_kernel))
 
     numerical_tmi = -np.array(numerical_tmi)
 
-    # Analytical TMI along axis (in nT)
+    # Get field direction cosines for analytical solution
+    l = _direction_cosines(inclination, declination)
+    lx, ly, lz = l[0], l[1], l[2]
+
+    # Analytical TMI for observation on z-axis above sphere center
+    # For a magnetic dipole with moment m = (4/3)πR³·(χ·F):
+    # The field components at (0, 0, r) are:
+    #   Bx = (μ₀/4π) · (m/r³) · 3·mx·mz / |m|
+    #   By = (μ₀/4π) · (m/r³) · 3·my·mz / |m|  
+    #   Bz = (μ₀/4π) · (m/r³) · (3·mz² - |m|²) / |m|
+    # Where m = [mx, my, mz] ∝ [lx, ly, lz]
+    #
+    # For induced magnetization: m ∝ [lx, ly, lz]
+    # TMI = Bx·lx + By·ly + Bz·lz
+    #
+    # After simplification:
+    # ΔT = (R³·χ·F) / (3·r³) · [3·lz² - 1 + 2·(lx² + ly²)]
+    # Since lx² + ly² + lz² = 1:
+    # ΔT = (R³·χ·F) / (3·r³) · [3·lz² - 1 + 2·(1 - lz²)]
+    # ΔT = (R³·χ·F) / (3·r³) · [3·lz² + 2 - 2·lz² - 1]
+    # ΔT = (R³·χ·F) / (3·r³) · [lz² + 1]  # Incorrect simplification
+    #
+    # Correct formula:
+    # ΔT = (R³·χ·F) / (3·r³) · [2·lz² - lx² - ly²]
+    # Or equivalently:
+    # ΔT = (R³·χ·F) / (3·r³) · [2·lz² - (1 - lz²)]
+    # ΔT = (R³·χ·F) / (3·r³) · [3·lz² - 1]
+    
     analytical_tmi = []
     for z in observation_heights:
         r = abs(z - center[2])
         if r <= R:
-            # Inside sphere (not expected in this test set). Use outer formula at R for continuity.
             r = R
-        dT = (2.0 / 3.0) * (R ** 3) * chi * intensity_nT / (r ** 3)
+        
+        # Correct analytical formula for induced dipole on axis
+        # ΔT = (R³·χ·F) / (3·r³) · [3·lz² - 1]
+        dT = (R ** 3) * chi * intensity_nT / (3.0 * r ** 3) * (3 * lz**2 - 1)
         analytical_tmi.append(dT)
 
     analytical_tmi = np.array(analytical_tmi)
 
     # Report
     print("\n=== Magnetics Sphere Benchmark (Induced-only TMI) ===")
+    print(f"Field: I={inclination}°, D={declination}°, F={intensity_nT} nT")
+    print(f"Direction cosines: lx={lx:.4f}, ly={ly:.4f}, lz={lz:.4f}")
     print(f"Sphere: R={R}m, center={center.tolist()}, χ={chi}")
     print(f"Observation heights: {observation_heights}")
     print(f"Numerical ΔT (nT):  {numerical_tmi}")
     print(f"Analytical ΔT (nT): {analytical_tmi}")
     if np.all(analytical_tmi != 0):
-        print(
-            f"Relative error (%): {np.abs((numerical_tmi - analytical_tmi) / analytical_tmi) * 100}"
-        )
+        print(f"Relative error (%): {np.abs((numerical_tmi - analytical_tmi) / analytical_tmi) * 100}")
+    else:
+        print(f"Absolute difference (nT): {np.abs(numerical_tmi - analytical_tmi)}")
 
     # Tolerance varies by field geometry (vertical fields are most accurate)
     np.testing.assert_allclose(
         numerical_tmi,
         analytical_tmi,
         rtol=rtol,
+        atol=1.0,  # Add absolute tolerance for near-zero values
         err_msg=(
             f"Magnetic TMI calculation deviates significantly from analytical sphere solution "
             f"for I={inclination}°, D={declination}°"
@@ -216,7 +245,7 @@ def test_magnetics_sphere_analytical_benchmark_induced_only(inclination, declina
 
 
 @pytest.mark.parametrize("inclination,declination,intensity_nT", [
-    (90.0, 0.0, 50000.0),
+        (90.0, 0.0, 50000.0),
 ])
 def test_magnetics_line_profile_symmetry_induced_only(inclination, declination, intensity_nT):
     """
@@ -236,9 +265,9 @@ def test_magnetics_line_profile_symmetry_induced_only(inclination, declination,
     z_obs = 600.0
 
     centers = np.column_stack([
-        x_profile,
-        np.full_like(x_profile, y_center),
-        np.full_like(x_profile, z_obs),
+            x_profile,
+            np.full_like(x_profile, y_center),
+            np.full_like(x_profile, z_obs),
     ])
 
     geophysics_grid = CenteredGrid(
@@ -248,9 +277,9 @@ def test_magnetics_line_profile_symmetry_induced_only(inclination, declination,
     )
 
     igrf_params = {
-        "inclination": inclination,
-        "declination": declination,
-        "intensity": intensity_nT,
+            "inclination": inclination,
+            "declination": declination,
+            "intensity"  : intensity_nT,
     }
 
     try:
@@ -291,7 +320,7 @@ def test_magnetics_line_profile_symmetry_induced_only(inclination, declination,
     center_idx = len(tmi_profile) // 2
 
     left = tmi_profile[:center_idx]
-    right = tmi_profile[center_idx + 1 :][::-1]
+    right = tmi_profile[center_idx + 1:][::-1]
 
     print("\n=== Magnetics Line Profile Symmetry (Induced-only TMI) ===")
     print(f"x_profile: {x_profile}")
@@ -322,9 +351,9 @@ def test_multiple_devices(n_devices):
     # Create a grid of observation points
     x_coords = np.linspace(400, 600, n_devices)
     centers = np.column_stack([
-        x_coords,
-        np.full(n_devices, 500.0),
-        np.full(n_devices, 600.0),
+            x_coords,
+            np.full(n_devices, 500.0),
+            np.full(n_devices, 600.0),
     ])
 
     grid = CenteredGrid(
@@ -476,7 +505,8 @@ def test_kernel_decay_with_distance(distance_multiplier):
                                err_msg=f"Decay should be approximately 1/r³ (expected {expected_ratio:.2f}, got {actual_ratio:.2f})")
 
     print(f"✓ Decay follows 1/r³ law within tolerance")
-    
+
+
 # =============================================================================
 # Kernel Property Tests
 # =============================================================================
@@ -500,16 +530,16 @@ def test_kernel_symmetry():
     # For resolution [rx, ry, rz]: actual points are [rx+1, ry+1, rz+1]
     nx = 2 * (20 // 2) + 1  # 21
     ny = 2 * (20 // 2) + 1  # 21
-    nz = 20 + 1             # 21
-    
+    nz = 20 + 1  # 21
+
     # Verify expected shape
     expected_voxels = nx * ny * nz
     assert tmi_kernel.shape[0] == expected_voxels, \
         f"Expected {expected_voxels} voxels, got {tmi_kernel.shape[0]}"
 
     # Reshape kernel to 3D grid (z, y, x ordering for numpy convention)
     kernel_3d = tmi_kernel.reshape((nx, ny, nz))
-    kernel_3d[:,:,0]
+    kernel_3d[:, :, 0]
 
     # For vertical field, kernel should be symmetric about vertical axis
     # Check horizontal slices are approximately radially symmetric
@@ -518,8 +548,8 @@ def test_kernel_symmetry():
 
     # Check that corners are approximately equal (radial symmetry)
     corners = [
-        slice_mid[0, 0], slice_mid[0, -1], 
-        slice_mid[-1, 0], slice_mid[-1, -1]
+            slice_mid[0, 0], slice_mid[0, -1],
+            slice_mid[-1, 0], slice_mid[-1, -1]
     ]
 
     # All corners should be similar for radial symmetry
@@ -536,9 +566,9 @@ def test_v_components_reusability():
 
     # Different IGRF scenarios
     igrf_scenarios = [
-        {"inclination": 90.0, "declination": 0.0, "intensity": 50000.0},
-        {"inclination": 0.0, "declination": 0.0, "intensity": 45000.0},
-        {"inclination": 60.0, "declination": 30.0, "intensity": 48000.0},
+            {"inclination": 90.0, "declination": 0.0, "intensity": 50000.0},
+            {"inclination": 0.0, "declination": 0.0, "intensity": 45000.0},
+            {"inclination": 60.0, "declination": 30.0, "intensity": 48000.0},
     ]
 
     chi = np.full(V.shape[1], 0.01)