A major rewrite.

Split the code into multiple files, drop the GSL dependency, introduce
configurable logging, random cleanups.

1 files changed, 57 insertions, 32 deletions

diff --git a/brill_data.h b/brill_data.h
index 37b338a..e928db9 100644
--- a/brill_data.h
+++ b/brill_data.h
@@ -18,8 +18,9 @@
 #ifndef BRILL_DATA_H
 #define BRILL_DATA_H
 
-#include <stdint.h>
+#include <stdarg.h>
 #include <stddef.h>
+#include <stdint.h>
 
 enum BDQFuncType {
     /**
@@ -32,9 +33,21 @@ enum BDQFuncType {
     BD_Q_FUNC_EPPLEY,
 };
 
+enum BDMetricComponent {
+    BD_METRIC_COMPONENT_RHORHO,
+    BD_METRIC_COMPONENT_RHOZ,
+    BD_METRIC_COMPONENT_RHOPHI,
+    BD_METRIC_COMPONENT_ZRHO,
+    BD_METRIC_COMPONENT_ZZ,
+    BD_METRIC_COMPONENT_ZPHI,
+    BD_METRIC_COMPONENT_PHIRHO,
+    BD_METRIC_COMPONENT_PHIZ,
+    BD_METRIC_COMPONENT_PHIPHI,
+};
+
 typedef struct BDContext {
     /**
-     *  private data
+     *  Solver internals, not to be accessed by the caller
      */
     void *priv;
 
@@ -86,26 +99,27 @@ typedef struct BDContext {
     int overdet_z;
 
     /**
-     * The difference between the index of the collocation grid used for the ρ
-     * direction and nb_coeffs_rho. The 'extra' collocation points furthest from
-     * the origin are discarded.
+     * The scaling factor used in the basis functions in the ρ direction.
      * Defaults to 3.
      */
-    unsigned int colloc_grid_offset_rho;
+    double basis_scale_factor_rho;
     /**
-     * Same as colloc_grid_offset_rho, but for the z direction.
+     * Same as basis_scale_factor_rho, but for the z direction
      */
-    unsigned int colloc_grid_offset_z;
+    double basis_scale_factor_z;
 
     /**
-     * The scaling factor used in the basis functions in the ρ direction.
-     * Defaults to 3.
+     * A callback that will be used to print diagnostic messages.
+     *
+     * Defaults to fprintf(stderr, ...)
      */
-    double basis_scale_factor_rho;
+    void (*log_callback)(const struct BDContext *bd, int level,
+                         const char *fmt, va_list);
+
     /**
-     * Same as basis_scale_factor_rho, but for the z direction
+     * Arbitrary user data, e.g. to be used by the log callback.
      */
-    double basis_scale_factor_z;
+    void *opaque;
 
     /**********
      * output *
@@ -121,7 +135,7 @@ typedef struct BDContext {
     /**
      * The number of array elements between two rows.
      */
-    ptrdiff_t stride;
+    unsigned int stride;
 } BDContext;
 
 /**
@@ -155,15 +169,17 @@ int bd_solve(BDContext *bd);
  *                   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.
+ *            I.e. psi[j * psi_stride + i] = ψ(rho[i], z[j]). The length of psi must be
+ *            at least psi_stride * nb_coords_z.
+ * @param psi_stride the distance (in double-sized elements) in psi between two elements corresponding
+ *                   to the same ρ but one step in z. Must be at least nb_coords_rho.
  *
  * @return >= 0 on success, a negative error code on failure.
  */
-int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho,
+int bd_eval_psi(const BDContext *bd, const double *rho, int nb_coords_rho,
                 const double *z, int nb_coords_z,
                 const unsigned int diff_order[2],
-                double *psi);
+                double *psi, unsigned int psi_stride);
 
 /**
  * Evaluate the 3-metric γ_ij at the specified rectangular grid (in cylindrical
@@ -174,21 +190,28 @@ int bd_eval_psi(BDContext *bd, const double *rho, int nb_coords_rho,
  * @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 comp the component of the metric to evaluate.
- * @param diff_order the order of the derivatives of the metric to evaluate. The first element
- *                   specifies the derivative wrt ρ, the second wrt z. I.e. diff_order = { 0, 0 }
- *                   evaluates the metric itself, diff_order = { 0, 1 } evaluates ∂γ/∂z etc.
- * @param out the array into which the values of the metric will be written. The metric
- *            is evaluated on the grid formed by the outer product of the rho
- *            and z vectors. I.e. out[j * nb_coords_rho + i] = γ_comp[0]comp[1](rho[i], z[j]).
- *            The length of out must be nb_coords_rho * nb_coords_z.
+ * @param nb_comp number of needed components of the metric
+ * @param comp a nb_comp-sized array specifying the components of the metric to evaluate
+ * @param diff_order a nb_comp-sized array specifying the order of the derivatives of the
+ *                   metric to evaluate. The first element specifies the derivative wrt ρ,
+ *                   the second wrt z. I.e. diff_order[i] = { 0, 0 } evaluates the metric
+ *                   itself, diff_order[i] = { 0, 1 } evaluates ∂γ/∂z etc.
+ * @param out a nb_comp-sized array of pointers to the arrays into which the values of the
+ *            metric will be written. Each requested component is evaluated on the grid
+ *            formed by the outer product of the rho and z vectors. I.e.
+ *            out[l][j * out_strides[l] + i] = γ_comp[l](rho[i], z[j]). The length of each
+ *            array in out must be nb_coords_rho * nb_coords_z.
+ * @param out_strides a nb_comp-sized array of distances (in double-sized elements), for each
+ *                    array in out, between two elements corresponding to the same ρ but one
+ *                    step in z. Each element in out_strides must be at least nb_coords_rho.
  *
  * @return >= 0 on success, a negative error code on failure.
  */
-int bd_eval_metric(BDContext *bd, const double *rho, int nb_coords_rho,
+int bd_eval_metric(const BDContext *bd, const double *rho, int nb_coords_rho,
                    const double *z, int nb_coords_z,
-                   const unsigned int comp[2], const unsigned int diff_order[2],
-                   double *out);
+                   int nb_comp, const enum BDMetricComponent *comp,
+                   const unsigned int (*diff_order)[2],
+                   double **out, unsigned int *out_strides);
 
 /**
  * Evaluate the q function at the specified rectangular grid (in cylindrical
@@ -204,13 +227,15 @@ int bd_eval_metric(BDContext *bd, const double *rho, int nb_coords_rho,
  *                   evaluates the q function itself, diff_order = { 0, 1 } evaluates ∂q/∂z etc.
  * @param out the array into which the values of the q function will be written. The q function
  *            is evaluated on the grid formed by the outer product of the rho
- *            and z vectors. I.e. out[j * nb_coords_rho + i] = γ_comp[0]comp[1](rho[i], z[j]).
+ *            and z vectors. I.e. out[j * out_stride + i] = γ_comp[0]comp[1](rho[i], z[j]).
  *            The length of out must be nb_coords_rho * nb_coords_z.
+ * @param out_stride the distance (in double-sized elements) in out between two elements corresponding
+ *                   to the same ρ but one step in z. Must be at least nb_coords_rho.
  *
  * @return >= 0 on success, a negative error code on failure.
  */
-int bd_eval_q(BDContext *bd, const double *rho, int nb_coords_rho,
+int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho,
               const double *z, int nb_coords_z, const unsigned int diff_order[2],
-              double *out);
+              double *out, unsigned int out_stride);
 
 #endif /* BRILL_DATA_H */