qfunc: disentangle the API from the global solver API

1 files changed, 103 insertions, 43 deletions

diff --git a/qfunc.c b/qfunc.c
index e7827df..82557ba 100644
--- a/qfunc.c
+++ b/qfunc.c
@@ -20,63 +20,74 @@
  * definitions of the q functions in the exponential in the Brill data
  */
 
+#include <errno.h>
 #include <math.h>
+#include <stdlib.h>
 
 #include "brill_data.h"
 #include "internal.h"
 
-static double q_gundlach(const BDContext *bd, double rho, double z)
+typedef double (*QFunc)(const struct QFuncContext *s, double rho, double z);
+
+struct QFuncContext {
+    const QFunc *qfunc;
+
+    double amplitude;
+    int n;
+};
+
+static double q_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * SQR(rho) * exp(- (SQR(rho) + SQR(z)));
+    return s->amplitude * SQR(rho) * exp(-(SQR(rho) + SQR(z)));
 }
 
-static double dq_rho_gundlach(const BDContext *bd, double rho, double z)
+static double dq_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return s->amplitude * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
 }
 
-static double d2q_rho_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_rho_gundlach(const QFuncContext *s, double rho, double z)
 {
     double rho2 = SQR(rho);
-    return bd->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
+    return s->amplitude * 2 * exp(-rho2 - SQR(z)) * ((1 - rho2) * (1 - 2 * rho2) - 2 * rho2);
 }
 
-static double d2q_rho_z_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_rho_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return -bd->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
+    return -s->amplitude * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (1 - SQR(rho));
 }
 
-static double dq_z_gundlach(const BDContext *bd, double rho, double z)
+static double dq_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return - bd->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
+    return -s->amplitude * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z));
 }
 
-static double d2q_z_gundlach(const BDContext *bd, double rho, double z)
+static double d2q_z_gundlach(const QFuncContext *s, double rho, double z)
 {
-    return bd->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
+    return s->amplitude * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1);
 }
 
 // q function from from PHYSICAL REVIEW D 88, 103009 (2013)
 // with σ and ρ_0 hardcoded to 1 for now
-const QFunc bdi_q_func_gundlach = {
-    .q         = q_gundlach,
-    .dq_rho    = dq_rho_gundlach,
-    .dq_z      = dq_z_gundlach,
-    .d2q_rho   = d2q_rho_gundlach,
-    .d2q_z     = d2q_z_gundlach,
-    .d2q_rho_z = d2q_rho_z_gundlach,
+static const QFunc q_func_gundlach[9] = {
+    [0] = q_gundlach,
+    [1] = dq_rho_gundlach,
+    [3] = dq_z_gundlach,
+    [2] = d2q_rho_gundlach,
+    [6] = d2q_z_gundlach,
+    [4] = d2q_rho_z_gundlach,
 };
 
-static double q_eppley(const BDContext *bd, double rho, double z)
+static double q_eppley(const QFuncContext *s, double rho, double z)
 {
     double r2 = SQR(rho) + SQR(z);
-    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
+    return s->amplitude * SQR(rho) / (1.0 + pow(r2, s->n / 2.0));
 }
 
-static double dq_rho_eppley(const BDContext *bd, double rho, double z)
+static double dq_rho_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -87,10 +98,10 @@ static double dq_rho_eppley(const BDContext *bd, double rho, double z)
     return A * rho * (2 * (1 + rnm2 * r2) - n * rho2 * rnm2) / SQR(1 + rnm2 * r2);
 }
 
-static double dq_z_eppley(const BDContext *bd, double rho, double z)
+static double dq_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -101,10 +112,10 @@ static double dq_z_eppley(const BDContext *bd, double rho, double z)
     return - A * n * rho2 * z * rnm2 / SQR(1 + rnm2 * r2);
 }
 
-static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
+static double d2q_rho_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -117,10 +128,10 @@ static double d2q_rho_eppley(const BDContext *bd, double rho, double z)
     return A * (SQR(n) * SQR(rho2) * rnm4 * (rn - 1) - n * rho2 * rnm4 * (3 * rho2 + 5 * z2) * rnp1 + 2 * SQR(rnp1)) / (rnp1 * SQR(rnp1));
 }
 
-static double d2q_z_eppley(const BDContext *bd, double rho, double z)
+static double d2q_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -133,10 +144,10 @@ static double d2q_z_eppley(const BDContext *bd, double rho, double z)
     return A * n * rho2 * rnm4 * (z2 - rho2 - n * z2 + rn * ((1 + n) * z2 - rho2)) / (rnp1 * SQR(rnp1));
 }
 
-static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
+static double d2q_rho_z_eppley(const QFuncContext *s, double rho, double z)
 {
-    double A = bd->amplitude;
-    double n = bd->eppley_n;
+    double A = s->amplitude;
+    double n = s->n;
     double rho2 = SQR(rho);
     double z2   = SQR(z);
 
@@ -151,11 +162,60 @@ static double d2q_rho_z_eppley(const BDContext *bd, double rho, double z)
 
 // q function from from PHYSICAL REVIEW D 16, 1609 (1977)
 // with λ hardcoded to 1 for now
-const QFunc bdi_q_func_eppley = {
-    .q         = q_eppley,
-    .dq_rho    = dq_rho_eppley,
-    .dq_z      = dq_z_eppley,
-    .d2q_rho   = d2q_rho_eppley,
-    .d2q_z     = d2q_z_eppley,
-    .d2q_rho_z = d2q_rho_z_eppley,
+static const QFunc q_func_eppley[9] = {
+    [0] = q_eppley,
+    [1] = dq_rho_eppley,
+    [3] = dq_z_eppley,
+    [2] = d2q_rho_eppley,
+    [6] = d2q_z_eppley,
+    [4] = d2q_rho_z_eppley,
 };
+
+double bdi_qfunc_eval(QFuncContext *s, double rho, double z, int diff_idx)
+{
+    return s->qfunc[diff_idx](s, rho, z);
+}
+
+void bdi_qfunc_free(QFuncContext **ps)
+{
+    QFuncContext *s = *ps;
+    free(s);
+    *ps = NULL;
+}
+
+int bdi_qfunc_init(const BDContext *bd, QFuncContext **out, enum BDQFuncType type,
+                   double amplitude, int n)
+{
+    QFuncContext *s = calloc(1, sizeof(*s));
+    int ret = 0;
+
+    if (!s)
+        return -ENOMEM;
+
+    switch (type) {
+    case BD_Q_FUNC_GUNDLACH:
+        s->qfunc = q_func_gundlach;
+        break;
+    case BD_Q_FUNC_EPPLEY:
+        if (n < 4) {
+            bdi_log(bd, 0, "Invalid n: %d < 4\n", n);
+            ret = -EINVAL;
+            goto fail;
+        }
+        s->qfunc = q_func_eppley;
+        break;
+    default:
+        bdi_log(bd, 0, "Unknown q function type: %d\n", type);
+        ret = -EINVAL;
+        goto fail;
+    }
+
+    s->amplitude = amplitude;
+    s->n         = n;
+
+    *out = s;
+    return 0;
+fail:
+    free(s);
+    return ret;
+}