Commit 7e19751e authored by Kloss's avatar Kloss

fix problem when roots are exactly on knots

parent 8edf45e7
......@@ -938,6 +938,9 @@ class BandSketching:
if derivative_order == 0: # use tabulated values, faster
y = self.y[:, band]
dy = self.dy[:, band]
elif derivative_order == 1:
y = self.dy[:, band]
dy = self._func(self.x, derivative_order + 1)[:, band]
else:
y = self._func(self.x, derivative_order)[:, band]
dy = self._func(self.x, derivative_order + 1)[:, band]
......@@ -950,7 +953,22 @@ class BandSketching:
y[np.abs(y - y_mean) < ytol] = y_mean
dy[np.abs(dy) < 100 * ytol] = 0 # derivative dy less acurate than y
g = y - f # roots of g give the intersection points we are seeking
roots = self.interpolation(self.x, y - f, dy).roots()
# the interpolation/rootfinding can fail if roots are exactly on the
# knots (the gridpoints self.x). In that case we check if g(i) == 0
# and also if the sign on g(i - 1) and g(i + 1) is different
# (necessary condition for g(i) == 0). This condition is important
# to not find spurious zeros of g when it is actually flat.
g_signchange = [np.sign(g[i - 1] * g[i + 1]) == -1
for i in range(1, len(g) - 1)]
# at the endpoints, check only that g is not flat
g_signchange.append(np.sign(g[-2]) != 0)
g_signchange.insert(0, np.sign(g[1]) != 0)
roots_on_knots = self.x[np.logical_and(np.abs(g) < ytol, g_signchange)]
roots = np.append(roots, roots_on_knots)
roots = remove_nan(roots)
if self.period:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment