Add a function for evaluating ψ.

2 files changed, 89 insertions, 0 deletions

diff --git a/brill_data.c b/brill_data.c
index 6eeadb6..9902385 100644
--- a/brill_data.c
+++ b/brill_data.c
@@ -463,6 +463,72 @@ void bd_context_free(BDContext **pbd)
     *pbd = NULL;
 }
 
+int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho,
+                const double *z, int nb_coords_z,
+                const unsigned int diff_order[2],
+                double *psi)
+{
+    BDPriv *s = bd->priv;
+    const long double sfr = bd->basis_scale_factor_rho, sfz = bd->basis_scale_factor_z;
+    long double (*eval)(long double coord, int idx, long double sf);
+
+    long double *basis_val_rho, *basis_val_z;
+
+    long double add = 0.0;
+
+    if (diff_order[0] > 2 || diff_order[1] > 2) {
+        fprintf(stderr, "Derivatives of higher order than 2 are not supported.\n");
+        return -ENOSYS;
+    }
+
+    if (diff_order[0] == 0 && diff_order[1] == 0)
+        add = 1.0;
+
+    /* precompute the basis values on the grid points */
+    basis_val_rho = malloc(sizeof(*basis_val_rho) * bd->nb_coeffs_rho * nb_coords_rho);
+    basis_val_z   = malloc(sizeof(*basis_val_z)   * bd->nb_coeffs_z   * nb_coords_z);
+
+    switch (diff_order[0]) {
+    case 0: eval = s->basis->eval;       break;
+    case 1: eval = s->basis->eval_diff1; break;
+    case 2: eval = s->basis->eval_diff2; break;
+    }
+
+    for (int i = 0; i < nb_coords_rho; i++) {
+        double rrho = rho[i];
+        for (int j = 0; j < bd->nb_coeffs_rho; j++)
+            basis_val_rho[i * bd->nb_coeffs_rho + j] = eval(rrho, j, sfr);
+    }
+
+    switch (diff_order[1]) {
+    case 0: eval = s->basis->eval;       break;
+    case 1: eval = s->basis->eval_diff1; break;
+    case 2: eval = s->basis->eval_diff2; break;
+    }
+    for (int i = 0; i < nb_coords_z; i++) {
+        double zz   = z[i];
+        for (int j = 0; j < bd->nb_coeffs_z; j++)
+            basis_val_z[i * bd->nb_coeffs_z + j] = eval(zz, j, sfz);
+    }
+
+    for (int i = 0; i < nb_coords_z; i++)
+        for (int j = 0; j < nb_coords_rho; j++) {
+            long double ppsi = add;
+
+            for (int m = 0; m < bd->nb_coeffs_z; m++)
+                for (int n = 0; n < bd->nb_coeffs_rho; n++)
+                    ppsi += s->coeffs->data[m * bd->nb_coeffs_rho + n] * basis_val_rho[j * bd->nb_coeffs_rho + n] * basis_val_z[i * bd->nb_coeffs_z + m];
+
+            psi[i * nb_coords_rho + j] = ppsi;
+        }
+
+    free(basis_val_rho);
+    free(basis_val_z);
+
+    return 0;
+
+}
+
 #if 0
 void brill_data(CCTK_ARGUMENTS)
 {
diff --git a/brill_data.h b/brill_data.h
index b8b63e8..85e1aac 100644
--- a/brill_data.h
+++ b/brill_data.h
@@ -138,3 +138,26 @@ void bd_context_free(BDContext **bd);
  * @return >= 0 on success, a negative error code on failure
  */
 int bd_solve(BDContext *bd);
+
+/**
+ * Evaluate the conformal factor ψ or its derivatives at the specified rectangular grid.
+ *
+ * @param bd the solver context
+ * @param rho the array of ρ coordinates.
+ * @param nb_coords_rho the number of elements in rho.
+ * @param z   the array of z coordinates.
+ * @param nb_coords_z the number of elements in z.
+ * @param diff_order the order of the derivatives of ψ to evaluate. The first element specifies the
+ *                   derivative wrt ρ, the second wrt z. I.e. diff_order = { 0, 0 } evaluates ψ
+ *                   itself, diff_order = { 0, 1 } evaluates ∂ψ/∂z etc.
+ * @param psi the array into which the values of ψ will be written. ψ is evaluated on the grid
+ *            formed by the outer product of the rho and z vectors.
+ *            I.e. psi[j * nb_coords_rho + i] = ψ(rho[i], z[j]). The length of psi must be
+ *            nb_coords_rho * nb_coords_z.
+ *
+ * @return >= 0 on success, a negative error code on failure.
+ */
+int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho,
+                const double *z, int nb_coords_z,
+                const unsigned int diff_order[2],
+                double *psi);