Hide warnings for a few undefined points

Author: hagenw

sfs/fd/source.py

@@ -82,7 +82,10 @@ def point(omega, x0, grid, *, c=None):
     grid = util.as_xyz_components(grid)
 
     r = np.linalg.norm(grid - x0)
-    return 1 / (4*np.pi) * np.exp(-1j * k * r) / r
+    # If r is 0, the sound pressure is complex infinity
+    numerator = np.exp(-1j * k * r) / (4*np.pi)
+    with np.errstate(invalid='ignore', divide='ignore'):
+        return numerator / r
 
 
 def point_velocity(omega, x0, grid, *, c=None, rho0=None):
@@ -395,7 +398,9 @@ def point_image_sources(omega, x0, grid, L, *, max_order, coeffs=None, c=None):
     p = 0
     for position, strength in zip(xs, source_strengths):
         if strength != 0:
-            p += strength * point(omega, position, grid, c=c)
+            # point can be complex infinity
+            with np.errstate(invalid='ignore'):
+                p += strength * point(omega, position, grid, c=c)
 
     return p
 

sfs/td/source.py

@@ -79,13 +79,17 @@ def point(xs, signal, observation_time, grid, c=None):
     if c is None:
         c = default.c
     r = np.linalg.norm(grid - xs)
-    # evaluate g over grid
-    weights = 1 / (4 * np.pi * r)
+    # If r is +-0, the sound pressure is +-infinity
+    with np.errstate(divide='ignore'):
+        weights = 1 / (4 * np.pi * r)
     delays = r / c
     base_time = observation_time - signal_offset
-    return weights * np.interp(base_time - delays,
+    points_at_time = np.interp(base_time - delays,
                                np.arange(len(data)) / samplerate,
                                data, left=0, right=0)
+    # weights can be +-infinity
+    with np.errstate(invalid='ignore'):
+        return weights * points_at_time
 
 
 def point_image_sources(x0, signal, observation_time, grid, L, max_order,