Implement evaluating derivatives of k_rtheta.

3 files changed, 23 insertions, 2 deletions

diff --git a/init.c b/init.c
index e5ad3d8..6b05048 100644
--- a/init.c
+++ b/init.c
@@ -830,9 +830,16 @@ int td_eval_krt(const TDContext *td,
     static const unsigned int nb_dummy_coords[2] = { 1, 1 };
 
     TDConstraintEvalContext *ce;
+    double (*eval)(TDConstraintEvalContext*, double, double);
     int ret;
 
-    if (diff_order[0] || diff_order[1])
+    if (diff_order[0] == 0 && diff_order[1] == 0)
+        eval = tdi_constraint_eval_k_rtheta;
+    else if (diff_order[0] == 1 && diff_order[1] == 0)
+        eval = tdi_constraint_eval_dk_rtheta_r;
+    else if (diff_order[0] == 0 && diff_order[1] == 1)
+        eval = tdi_constraint_eval_dk_rtheta_t;
+    else
         return -ENOSYS;
 
     ret = constraint_eval_alloc(td, nb_dummy_coords, dummy_coords,
@@ -844,7 +851,7 @@ int td_eval_krt(const TDContext *td,
         double theta_val = theta[i];
         double r_val     = r[i];
 
-        out[i] = tdi_constraint_eval_k_rtheta(ce, r_val, theta_val);
+        out[i] = eval(ce, r_val, theta_val);
     }
 
     tdi_constraint_eval_free(&ce);
diff --git a/td_constraints.c b/td_constraints.c
index dd62150..41c1c5b 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -218,6 +218,18 @@ double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, doub
     return fd->eval_krt(ctx, r, theta);
 }
 
+double tdi_constraint_eval_dk_rtheta_r(TDConstraintEvalContext *ctx, double r, double theta)
+{
+    const TDFamilyDef *fd = td_families[ctx->family];
+    return fd->eval_krt_dr(ctx, r, theta);
+}
+
+double tdi_constraint_eval_dk_rtheta_t(TDConstraintEvalContext *ctx, double r, double theta)
+{
+    const TDFamilyDef *fd = td_families[ctx->family];
+    return fd->eval_krt_dt(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 259cbfb..a2e43e2 100644
--- a/td_constraints.h
+++ b/td_constraints.h
@@ -69,5 +69,7 @@ int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
                         double *out);
 
 double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, double theta);
+double tdi_constraint_eval_dk_rtheta_r(TDConstraintEvalContext *ctx, double r, double theta);
+double tdi_constraint_eval_dk_rtheta_t(TDConstraintEvalContext *ctx, double r, double theta);
 
 #endif // TEUKOLSKY_DATA_CONSTRAINTS_H