From 33eab328cdd9be877f9a3767f8d776194b9dab36 Mon Sep 17 00:00:00 2001
From: Hagen Wierstorf <hagen.wierstorf@posteo.de>
Date: Mon, 18 Mar 2019 10:55:17 +0100
Subject: [PATCH] Hide warnings for a few undefined points

---
 sfs/fd/source.py |  9 +++++++--
 sfs/td/source.py | 10 +++++++---
 2 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/sfs/fd/source.py b/sfs/fd/source.py
index 993bcba..3a32b86 100644
--- a/sfs/fd/source.py
+++ b/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
 
diff --git a/sfs/td/source.py b/sfs/td/source.py
index bf9e9aa..017aeb4 100644
--- a/sfs/td/source.py
+++ b/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,