qfunc: disentangle the API from the global solver API

1 files changed, 11 insertions, 25 deletions

diff --git a/eval.c b/eval.c
index 6c70790..98fd0a7 100644
--- a/eval.c
+++ b/eval.c
@@ -1,5 +1,5 @@
 /*
- * Copyright 2014-2015 Anton Khirnov <anton@khirnov.net>
+ * Copyright 2014-2016 Anton Khirnov <anton@khirnov.net>
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -27,6 +27,7 @@
 
 #include "brill_data.h"
 #include "internal.h"
+#include "qfunc.h"
 
 int bd_eval_psi(const BDContext *bd,
                 const double *rho, int nb_coords_rho,
@@ -172,14 +173,14 @@ do {                                                                \
  * derivatives as needed */
 #define DO_EVAL_Q_WRAPPER(DO_EVAL, eval_drho, eval_dz, eval_d2rho, eval_d2z, eval_d2rhoz)                                                    \
 do {                                                                                                                                         \
-    const double base = psi4 * exp(2 * s->q_func->q(bd, rrho, zz));                                                                          \
+    const double base = psi4 * exp(2 * bdi_qfunc_eval(s->qfunc, rrho, zz, 0));                                                               \
     double drho, dz, d2rho, d2z, d2rho_z;                                                                                                    \
                                                                                                                                              \
-    if (eval_drho)   drho    = s->q_func->dq_rho(bd, rrho, zz)    + 2 * GRID(dpsi_rho)    / ppsi;                                            \
-    if (eval_dz)     dz      = s->q_func->dq_z(bd, rrho, zz)      + 2 * GRID(dpsi_z)      / ppsi;                                            \
-    if (eval_d2rho)  d2rho   = s->q_func->d2q_rho(bd, rrho, zz)   + 2 * GRID(d2psi_rho)   / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi);           \
-    if (eval_d2z)    d2z     = s->q_func->d2q_z(bd, rrho, zz)     + 2 * GRID(d2psi_z)     / ppsi - 2 * SQR(GRID(dpsi_z)   / ppsi);           \
-    if (eval_d2rhoz) d2rho_z = s->q_func->d2q_rho_z(bd, rrho, zz) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \
+    if (eval_drho)   drho    = bdi_qfunc_eval(s->qfunc, rrho, zz, 1) + 2 * GRID(dpsi_rho)    / ppsi;                                            \
+    if (eval_dz)     dz      = bdi_qfunc_eval(s->qfunc, rrho, zz, 3) + 2 * GRID(dpsi_z)      / ppsi;                                            \
+    if (eval_d2rho)  d2rho   = bdi_qfunc_eval(s->qfunc, rrho, zz, 2) + 2 * GRID(d2psi_rho)   / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi);           \
+    if (eval_d2z)    d2z     = bdi_qfunc_eval(s->qfunc, rrho, zz, 6) + 2 * GRID(d2psi_z)     / ppsi - 2 * SQR(GRID(dpsi_z)   / ppsi);           \
+    if (eval_d2rhoz) d2rho_z = bdi_qfunc_eval(s->qfunc, rrho, zz, 4) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \
                                                                                                                                              \
     DO_EVAL;                                                                                                                                 \
 } while (0)
@@ -286,34 +287,19 @@ int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho,
               double *out, unsigned int out_stride)
 {
     BDPriv *s = bd->priv;
-    double (*eval)(const struct BDContext *bd, double rho, double z);
+    int diff_idx;
 
     if (diff_order[0] > 2 || diff_order[1] > 2) {
         bdi_log(bd, 0, "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;
-    }
+    diff_idx = diff_order[1] * 3 + diff_order[0];
 
     for (int i = 0; i < nb_coords_z; i++) {
         double *dst = out + i * out_stride;
         for (int j = 0; j < nb_coords_rho; j++)
-            *dst++ = eval(bd, rho[j], z[i]);
+            *dst++ = bdi_qfunc_eval(s->qfunc, rho[j], z[i], diff_idx);
     }
 
     return 0;