From 1f443d3fdb93e71a70244b7521ae0f4c40601351 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 27 Aug 2021 13:35:52 +0200 Subject: Add public API for evaluating the residual at an arbitrary location. --- teukolsky_data.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) (limited to 'teukolsky_data.py') diff --git a/teukolsky_data.py b/teukolsky_data.py index 00ec5c4..a9bd490 100644 --- a/teukolsky_data.py +++ b/teukolsky_data.py @@ -170,6 +170,27 @@ class TeukolskyData(object): return metric, curv + def eval_residual(self, r, theta): + if r.ndim != 1 or theta.ndim != 1: + raise TypeError('r and theta must both be 1-dimensional NumPy arrays') + + out = np.empty((3, r.shape[0] * theta.shape[0])) + + c_out = (ctypes.POINTER(ctypes.c_double) * 3)() + for i in range(3): + c_out[i] = ctypes.cast(np.ctypeslib.as_ctypes(out[i]), ctypes.POINTER(ctypes.c_double)) + + c_r = np.ctypeslib.as_ctypes(r) + c_theta = np.ctypeslib.as_ctypes(theta) + + ret = self._libtd.td_eval_residual(self._tdctx, r.shape[0], c_r, + theta.shape[0], c_theta, c_out) + if ret < 0: + raise RuntimeError('Error evaluating the variable: %d' % ret) + + out.shape = (3, theta.shape[0], r.shape[0]) + return out + @property def amplitude(self): return self._tdctx.contents.__getattribute__('amplitude') -- cgit v1.2.3