td_constraints: tie the constraint equations to the family

Preparation for the ingoing pulse family, which has a slightly different
metric expression and so different form of the constraint equations.

1 files changed, 83 insertions, 76 deletions

diff --git a/td_constraints.c b/td_constraints.c
index 41c1c5b..5618cb3 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -27,79 +27,19 @@ struct TDConstraintEvalPriv {
     double *dk_rtheta[2];
 };
 
+typedef void (*ConstraintEvalCB)(TDConstraintEvalContext *ctx,
+                                 double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
+                                 double *out);
+
 typedef struct TDFamilyDef {
     double (*eval_krt)   (TDConstraintEvalContext *ctx, double r, double theta);
     double (*eval_krt_dr)(TDConstraintEvalContext *ctx, double r, double theta);
     double (*eval_krt_dt)(TDConstraintEvalContext *ctx, double r, double theta);
     double a_converge;
     double a_diverge;
-} TDFamilyDef;
-
-static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx,
-                                            double r, double theta)
-{
-    const double r2 = SQR(r);
-    const double r3 = r * r2;
-    return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
-}
-
-static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx,
-                                               double r, double theta)
-{
-    const double r2 = SQR(r);
-    const double r4 = SQR(r2);
-    return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
-}
-
-static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx,
-                                               double r, double theta)
-{
-    const double r2 = SQR(r);
-    const double r3 = r * r2;
-    return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
-}
-
-static const TDFamilyDef ae_time_antisym = {
-    .eval_krt    = k_rtheta_eval_ae_time_antisym,
-    .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym,
-    .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym,
-    .a_converge  = 0.0092078837420,
-    .a_diverge   = 0.0092078837427,
-};
-
-static double k_rtheta_eval_simple_time_antisym(TDConstraintEvalContext *ctx,
-                                                double r, double theta)
-{
-    const double r2 = SQR(r);
-    return ctx->amplitude * r * exp(-r2) * sin(2.0 * theta);
-}
-
-static double k_rtheta_eval_dr_simple_time_antisym(TDConstraintEvalContext *ctx,
-                                                   double r, double theta)
-{
-    const double r2 = SQR(r);
-    return ctx->amplitude * (1.0 - 2.0 * r2) * exp(-r2) * sin(2.0 * theta);
-}
-
-static double k_rtheta_eval_dt_simple_time_antisym(TDConstraintEvalContext *ctx,
-                                                   double r, double theta)
-{
-    const double r2 = SQR(r);
-    return 2.0 * ctx->amplitude * r * exp(-r2) * cos(2.0 * theta);
-}
 
-static const TDFamilyDef simple_time_antisym = {
-    .eval_krt    = k_rtheta_eval_simple_time_antisym,
-    .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym,
-    .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym,
-    .a_converge  = 1.018918628476058,
-    .a_diverge   = 1.018918628476968,
-};
-
-static const TDFamilyDef *td_families[] = {
-    [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym,
-    [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym,
-};
+    const ConstraintEvalCB *constraint_eval;
+} TDFamilyDef;
 
 static void
 constraint_eval_ham_confflat(TDConstraintEvalContext *ctx,
@@ -212,6 +152,81 @@ constraint_eval_mom_t_confflat(TDConstraintEvalContext *ctx,
     }
 }
 
+static const ConstraintEvalCB constraint_funcs_confflat[TD_CONSTRAINT_EQ_NB] = {
+    [TD_CONSTRAINT_EQ_HAM]   = constraint_eval_ham_confflat,
+    [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat,
+    [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat,
+};
+
+static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                            double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r3 = r * r2;
+    return -60.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                               double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r4 = SQR(r2);
+    return -60.0 * sqrt(2) * (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx,
+                                               double r, double theta)
+{
+    const double r2 = SQR(r);
+    const double r3 = r * r2;
+    return -120.0 * sqrt(2) * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude / sqrt(M_PI);
+}
+
+static const TDFamilyDef ae_time_antisym = {
+    .eval_krt    = k_rtheta_eval_ae_time_antisym,
+    .eval_krt_dr = k_rtheta_eval_dr_ae_time_antisym,
+    .eval_krt_dt = k_rtheta_eval_dt_ae_time_antisym,
+    .a_converge  = 0.0092078837420,
+    .a_diverge   = 0.0092078837427,
+    .constraint_eval = constraint_funcs_confflat,
+};
+
+static double k_rtheta_eval_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                double r, double theta)
+{
+    const double r2 = SQR(r);
+    return ctx->amplitude * r * exp(-r2) * sin(2.0 * theta);
+}
+
+static double k_rtheta_eval_dr_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                   double r, double theta)
+{
+    const double r2 = SQR(r);
+    return ctx->amplitude * (1.0 - 2.0 * r2) * exp(-r2) * sin(2.0 * theta);
+}
+
+static double k_rtheta_eval_dt_simple_time_antisym(TDConstraintEvalContext *ctx,
+                                                   double r, double theta)
+{
+    const double r2 = SQR(r);
+    return 2.0 * ctx->amplitude * r * exp(-r2) * cos(2.0 * theta);
+}
+
+static const TDFamilyDef simple_time_antisym = {
+    .eval_krt    = k_rtheta_eval_simple_time_antisym,
+    .eval_krt_dr = k_rtheta_eval_dr_simple_time_antisym,
+    .eval_krt_dt = k_rtheta_eval_dt_simple_time_antisym,
+    //.a_converge  = 1.018918628476058,
+    .a_converge  = 1.0189186,
+    .a_diverge   = 1.018918628476968,
+    .constraint_eval = constraint_funcs_confflat,
+};
+
+static const TDFamilyDef *td_families[] = {
+    [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym,
+    [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym,
+};
+
 double tdi_constraint_eval_k_rtheta(TDConstraintEvalContext *ctx, double r, double theta)
 {
     const TDFamilyDef *fd = td_families[ctx->family];
@@ -310,19 +325,11 @@ void tdi_constraint_eval_free(TDConstraintEvalContext **pctx)
     *pctx = NULL;
 }
 
-static const void
-(*constraint_funcs[TD_CONSTRAINT_EQ_NB])(TDConstraintEvalContext *ctx,
-                                         double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
-                                         double *out) = {
-    [TD_CONSTRAINT_EQ_HAM]   = constraint_eval_ham_confflat,
-    [TD_CONSTRAINT_EQ_MOM_0] = constraint_eval_mom_r_confflat,
-    [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat,
-};
-
 int tdi_constraint_eval(TDConstraintEvalContext *ctx, enum TDConstraintEq eq,
                         const double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB],
                         double *out)
 {
-    constraint_funcs[eq](ctx, vars, out);
+    const TDFamilyDef *fd = td_families[ctx->family];
+    fd->constraint_eval[eq](ctx, vars, out);
     return 0;
 }