Add public API for evaluating the residual at an arbitrary location.

1 files changed, 21 insertions, 0 deletions

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')