Merge pull request #81 from jcapriot/circle_fix

Fix test for small rho

.github/workflows/release.yml

@@ -64,6 +64,7 @@ jobs:
         with:
           path: ./dist/geoana*.whl
 
+
   distribute:
     name: Distribute documentation
     runs-on: ubuntu-latest

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] * (

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)
+    npt.assert_allclose(test, other, rtol=1E-3, atol=1E-20)