From e2745670894cac8ca73c6bee5c1b4cc2c8c94a3e Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Thu, 30 Jan 2020 13:36:41 +0100
Subject: teukolsky_data: change TD_FAMILY_AE_TIME_ANTISYM

This is obviously not the data that A&E were using in their paper
(though the reason is unknown, as the free functions should be the
same), so stop referring to it as such. Also drop the constant
multiplier of 60 * \sqrt{2/pi}, since it has no other purpose than being
consistent with their data.
Name the resulting family TD_FAMILY_TIME_ANTISYM_CUBIC.

This is an API break.
---
 init.c            |  2 +-
 td_constraints.c  | 30 +++++++++++++++---------------
 teukolsky_data.h  | 12 +++++++-----
 teukolsky_data.py |  2 +-
 4 files changed, 24 insertions(+), 22 deletions(-)

diff --git a/init.c b/init.c
index 6b05048..4df0f8f 100644
--- a/init.c
+++ b/init.c
@@ -453,7 +453,7 @@ TDContext *td_context_alloc(void)
 
     td->nb_threads = 1;
 
-    td->family          = TD_FAMILY_SIMPLE_TIME_ANTISYM;
+    td->family          = TD_FAMILY_TIME_ANTISYM_CUBIC;
     td->solution_branch = 0;
     td->amplitude       = 1.0;
 
diff --git a/td_constraints.c b/td_constraints.c
index 5618cb3..73d48a7 100644
--- a/td_constraints.c
+++ b/td_constraints.c
@@ -158,36 +158,36 @@ static const ConstraintEvalCB constraint_funcs_confflat[TD_CONSTRAINT_EQ_NB] = {
     [TD_CONSTRAINT_EQ_MOM_1] = constraint_eval_mom_t_confflat,
 };
 
-static double k_rtheta_eval_ae_time_antisym(TDConstraintEvalContext *ctx,
-                                            double r, double theta)
+static double k_rtheta_eval_time_antisym_cubic(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);
+    return (3.0 * r - 2.0 * r3) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude;
 }
 
-static double k_rtheta_eval_dr_ae_time_antisym(TDConstraintEvalContext *ctx,
+static double k_rtheta_eval_dr_time_antisym_cubic(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);
+    return (3.0 - 12.0 * r2 + 4.0 * r4) * exp(-r2) * sin(2.0 * theta) * ctx->amplitude;
 }
 
-static double k_rtheta_eval_dt_ae_time_antisym(TDConstraintEvalContext *ctx,
-                                               double r, double theta)
+static double k_rtheta_eval_dt_time_antisym_cubic(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);
+    return 2.0 * (3.0 * r - 2.0 * r3) * exp(-r2) * cos(2.0 * theta) * ctx->amplitude;
 }
 
-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 const TDFamilyDef time_antisym_cubic = {
+    .eval_krt    = k_rtheta_eval_time_antisym_cubic,
+    .eval_krt_dr = k_rtheta_eval_dr_time_antisym_cubic,
+    .eval_krt_dt = k_rtheta_eval_dt_time_antisym_cubic,
+    .a_converge  = 0.440809326171875,
+    .a_diverge   = 0.4408099365234375,
     .constraint_eval = constraint_funcs_confflat,
 };
 
@@ -223,7 +223,7 @@ static const TDFamilyDef simple_time_antisym = {
 };
 
 static const TDFamilyDef *td_families[] = {
-    [TD_FAMILY_AE_TIME_ANTISYM] = &ae_time_antisym,
+    [TD_FAMILY_TIME_ANTISYM_CUBIC] = &time_antisym_cubic,
     [TD_FAMILY_SIMPLE_TIME_ANTISYM] = &simple_time_antisym,
 };
 
diff --git a/teukolsky_data.h b/teukolsky_data.h
index 4967f8f..943407b 100644
--- a/teukolsky_data.h
+++ b/teukolsky_data.h
@@ -41,13 +41,15 @@
  */
 enum TDFamily {
     /**
-     * The time-antisymmetric initial data used in Abrahams&Evans PhysRevD v49,n8 (1994).
+     * Time-antisymmetric initial data with a cubic radial term multiplied by an
+     * exponential and a quadrupole angular term.
+     * Similar to that 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     /
+     *  r       /   x         x  3 \        /    x  2  \
+     * K  =  a | 3 ---  - 2 (---)   | * exp| - (---)    | sin(2θ)
+     *  θ       \   L         L    /        \    L     /
      */
-    TD_FAMILY_AE_TIME_ANTISYM = 0,
+    TD_FAMILY_TIME_ANTISYM_CUBIC = 0,
     /**
      * Simpler time-antisymmetric initial data.
      * Conformally flat spatial metric.
diff --git a/teukolsky_data.py b/teukolsky_data.py
index bb6ba92..2b6198e 100644
--- a/teukolsky_data.py
+++ b/teukolsky_data.py
@@ -21,7 +21,7 @@ import numpy as np
 
 class TeukolskyData(object):
 
-    TD_FAMILY_AE_TIME_ANTISYM     = 0
+    TD_FAMILY_TIME_ANTISYM_CUBIC  = 0
     TD_FAMILY_SIMPLE_TIME_ANTISYM = 1
 
     coeffs = None
-- 
cgit v1.2.3