diff options
Diffstat (limited to 'brill_data.c')
-rw-r--r-- | brill_data.c | 36 |
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; +} |