From 2ee593d56a83f0264cc41620c2feb3e03a13974f Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Sun, 22 Apr 2018 17:44:04 +0200
Subject: Add a function for evaluating K_rtheta

---
 init.c            | 27 +++++++++++++++++++++++++++
 td_constraints.c  |  6 ++++++
 td_constraints.h  |  2 ++
 teukolsky_data.h  |  4 ++++
 teukolsky_data.py |  2 ++
 5 files changed, 41 insertions(+)

diff --git a/init.c b/init.c
index 5050f5e..a7d54a8 100644
--- a/init.c
+++ b/init.c
@@ -583,3 +583,30 @@ int td_eval_kpp(const TDContext *td,
 {
     return eval_var(td, 2, nb_coords, r, theta, diff_order, out);
 }
+
+int td_eval_krt(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out)
+{
+    TDConstraintEvalContext *ce;
+    int ret;
+
+    if (diff_order[0] || diff_order[1])
+        return -ENOSYS;
+
+    ret = constraint_eval_alloc(td, td->amplitude, &ce);
+    if (ret < 0)
+        return ret;
+
+    for (int i = 0; i < nb_coords; i++) {
+        double theta_val = theta[i];
+        double r_val     = r[i];
+
+        out[i] = tdi_constraint_eval_k_rtheta(ce, r_val, theta_val);
+    }
+
+    tdi_constraint_eval_free(&ce);
+
+    return 0;
+}
diff --git a/td_constraints.c b/td_constraints.c
index fe6bc18..663b83a 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -192,6 +192,12 @@ constraint_eval_mom_t_confflat(TDConstraintEvalContext *ctx,
     }
 }
 
+double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, double theta)
+{
+    const TDFamilyDef *fd = td_families[ctx->family];
+    return fd->eval_krt(ctx, r, theta);
+}
+
 int tdi_constraint_eval_init(TDConstraintEvalContext *ctx)
 {
     TDConstraintEvalPriv *priv = ctx->priv;
diff --git a/td_constraints.h b/td_constraints.h
index f0a892c..6746c3e 100644
--- a/td_constraints.h
+++ b/td_constraints.h
@@ -81,4 +81,6 @@ int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
                         const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
                         double *out);
 
+double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, double theta);
+
 #endif // TEUKOLSKY_DATA_CONSTRAINTS_H
diff --git a/teukolsky_data.h b/teukolsky_data.h
index 1668bff..e0877cd 100644
--- a/teukolsky_data.h
+++ b/teukolsky_data.h
@@ -163,5 +163,9 @@ int td_eval_kpp(const TDContext *td,
                 size_t nb_coords, const double *r, const double *theta,
                 const unsigned int diff_order[2],
                 double *out);
+int td_eval_krt(const TDContext *td,
+                size_t nb_coords, const double *r, const double *theta,
+                const unsigned int diff_order[2],
+                double *out);
 
 #endif /* TEUKOLSKY_DATA_H */
diff --git a/teukolsky_data.py b/teukolsky_data.py
index c9150de..f6c78c5 100644
--- a/teukolsky_data.py
+++ b/teukolsky_data.py
@@ -122,6 +122,8 @@ class TeukolskyData(object):
         return self._eval_var(self._libtd.td_eval_krr, r, theta, diff_order)
     def eval_kpp(self, r, theta, diff_order = None):
         return self._eval_var(self._libtd.td_eval_kpp, r, theta, diff_order)
+    def eval_krt(self, r, theta, diff_order = None):
+        return self._eval_var(self._libtd.td_eval_krt, r, theta, diff_order)
 
 
     @property
-- 
cgit v1.2.3