Add a function for evaluating q.

1 files changed, 36 insertions, 0 deletions

diff --git a/brill_data.c b/brill_data.c
index 1f6378f..96b2d3a 100644
--- a/brill_data.c
+++ b/brill_data.c
@@ -803,3 +803,39 @@ fail:
 
     return ret;
 }
+
+int bd_eval_q(BDContext *bd, const double *rho, int nb_coords_rho,
+              const double *z, int nb_coords_z, const unsigned int diff_order[2],
+              double *out)
+{
+    BDPriv *s = bd->priv;
+    double (*eval)(struct BDContext *bd, double rho, double z);
+
+    if (diff_order[0] > 2 || diff_order[1] > 2) {
+        fprintf(stderr, "Derivatives of higher order than 2 are not supported.\n");
+        return -ENOSYS;
+    }
+
+    switch (diff_order[0]) {
+    case 0:
+        switch (diff_order[1]) {
+        case 0: eval = s->q_func->q;     break;
+        case 1: eval = s->q_func->dq_z;  break;
+        case 2: eval = s->q_func->d2q_z; break;
+        }
+        break;
+    case 1:
+        switch (diff_order[1]) {
+        case 0: eval = s->q_func->dq_rho;    break;
+        case 1: eval = s->q_func->d2q_rho_z; break;
+        }
+        break;
+    case 2:     eval = s->q_func->d2q_rho;   break;
+    }
+
+    for (int i = 0; i < nb_coords_z; i++)
+        for (int j = 0; j < nb_coords_rho; j++)
+            out[i * nb_coords_rho + j] = eval(bd, rho[j], z[i]);
+
+    return 0;
+}