From 4cb4a1720e4ac6bd96d6b33400dfc7a8f23aed2d Mon Sep 17 00:00:00 2001
From: Joseph Capriotti <josephrcapriotti@gmail.com>
Date: Thu, 24 Oct 2024 17:31:10 -0600
Subject: [PATCH] fix test for small rho

Do a check for rho < sqrt(a**2 + z**2)
---
 .github/workflows/release.yml        |  1 +
 geoana/em/static/wholespace.py       | 12 ++++--------
 tests/em/static/test_current_loop.py | 13 +++++--------
 3 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 64fa0946..31ac9703 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -64,6 +64,7 @@ jobs:
         with:
           path: ./dist/geoana*.whl
 
+
   distribute:
     name: Distribute documentation
     runs-on: ubuntu-latest
diff --git a/geoana/em/static/wholespace.py b/geoana/em/static/wholespace.py
index 4d9e828f..32612d56 100644
--- a/geoana/em/static/wholespace.py
+++ b/geoana/em/static/wholespace.py
@@ -597,10 +597,8 @@ def vector_potential(self, xyz, coordinates="cartesian"):
 
         A_cyl = np.zeros_like(r_vec)
 
-        # when rho is small relative to the radius and z, k_sq -> 0
-        # this is unstable so use a small argument approximation.
-        # check k_sq for the relative sizes
-        small_rho = k_sq < 1E-3
+        # when rho is small relative to the radius and z
+        small_rho = rho**2/(a**2 + z**2) < 1E-6
 
         temp = np.sqrt(a**2 + z[small_rho]**2)
         A_cyl[small_rho, 1] = np.pi * C * (
@@ -731,10 +729,8 @@ def magnetic_flux_density(self, xyz, coordinates="cartesian"):
 
         B_cyl = np.zeros_like(r_cyl)
 
-        # when rho is small relative to the radius and z, k_sq -> 0
-        # this is unstable so use a small argument approximation.
-        # check k_sq for the relative sizes
-        small_rho = k_sq < 1E-3
+        # when rho is small relative to the radius and z
+        small_rho = rho**2/(a**2 + z**2) < 1E-6
 
         temp = np.sqrt(a**2 + z[small_rho]**2)
         B_cyl[small_rho, 0] = 3 * C * np.pi * a**2 * z[small_rho] * (
diff --git a/tests/em/static/test_current_loop.py b/tests/em/static/test_current_loop.py
index b9bb08e0..403444d7 100644
--- a/tests/em/static/test_current_loop.py
+++ b/tests/em/static/test_current_loop.py
@@ -24,10 +24,10 @@ def circle_loop():
 
 @pytest.fixture()
 def xyz():
-    nx, ny, nz = (11, 13, 9)
-    x = np.linspace(-100, 100, nx)
-    y = np.linspace(-90, 90, ny)
-    z = np.linspace(-80, 80, nz)
+    nx, ny, nz = (32, 32, 32)
+    x = np.linspace(-50, 50, nx)
+    y = np.linspace(-50, 50, ny)
+    z = np.linspace(-50, 50, nz)
     xyz = np.meshgrid(x, y, z)
     return xyz
 
@@ -87,8 +87,5 @@ def test_circular_loop(circle_loop, orient, method, xyz):
     test = getattr(circle_loop, method)(xyz)
     other = getattr(loop_approx, method)(xyz)
 
-    atol = 1E-16
-    if method == 'magnetic_field':
-        atol /= mu_0
 
-    npt.assert_allclose(test, other, rtol=1E-3, atol=atol)
\ No newline at end of file
+    npt.assert_allclose(test, other, rtol=1E-3, atol=1E-20)
\ No newline at end of file