Rewrite constraint evaluation code.

Make it much more efficient and easier to add other seed functions.

1 files changed, 46 insertions, 5 deletions

diff --git a/td_constraints.h b/td_constraints.h
index 44f5257..f0a892c 100644
--- a/td_constraints.h
+++ b/td_constraints.h
@@ -20,6 +20,20 @@
 
 #include "pssolve.h"
 
+/**
+ * Identifies the specific initial data family to use.
+ */
+enum TDFamily {
+    /**
+     * The time-antisymmetric initial data used in Abrahams&Evans PhysRevD v49,n8 (1994).
+     * Conformally flat spatial metric.
+     *  r                / x        x  3 \        /    x  2  \
+     * K  = -60√(2/π) a | ---  -  (---)   | * exp| - (---)    | sin(2θ)
+     *  θ                \ L        L    /        \    L     /
+     */
+    TD_FAMILY_AE_TIME_ANTISYM = 0,
+};
+
 enum TDConstraintEq {
     TD_CONSTRAINT_EQ_HAM = 0,
     TD_CONSTRAINT_EQ_MOM_0,
@@ -34,10 +48,37 @@ enum TDConstraintVar {
     TD_CONSTRAINT_VAR_NB,
 };
 
-int tdi_constraints_eval(enum TDConstraintEq eq,
-                         const double a,
-                         unsigned int nb_coords[2], double *coords[2],
-                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
-                         double *out);
+
+typedef struct TDConstraintEvalPriv TDConstraintEvalPriv;
+
+typedef struct TDConstraintEvalContext {
+    /**
+     * Private data, not to be touched by the caller.
+     */
+    TDConstraintEvalPriv *priv;
+
+    /**
+     * The logging context.
+     * Set by the caller before tdi_constraint_eval_init().
+     */
+    TDLogger logger;
+
+    unsigned int nb_coords[2];
+    double *coords[2];
+
+    enum TDFamily family;
+
+    double amplitude;
+    double a_converge;
+    double a_diverge;
+} TDConstraintEvalContext;
+
+int tdi_constraint_eval_alloc(TDConstraintEvalContext **ctx, enum TDFamily family);
+int tdi_constraint_eval_init (TDConstraintEvalContext  *ctx);
+void tdi_constraint_eval_free(TDConstraintEvalContext **ctx);
+
+int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
+                        const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                        double *out);
 
 #endif // TEUKOLSKY_DATA_CONSTRAINTS_H