Add the Eppley(1977) version of q.

2 files changed, 29 insertions, 2 deletions

diff --git a/param.ccl b/param.ccl
index 9b214da..1f80535 100644
--- a/param.ccl
+++ b/param.ccl
@@ -13,6 +13,11 @@ CCTK_REAL amplitude "Wave amplitude A."
     0: :: ""
 } 1.0
 
+CCTK_INT eppley_n ""
+{
+    1: :: ""
+} 5
+
 CCTK_INT basis_order_r "Number of the basis functions in the radial direction"
 {
     1: :: ""
diff --git a/src/brill.c b/src/brill.c
index cc45587..a2cfbd5 100644
--- a/src/brill.c
+++ b/src/brill.c
@@ -100,6 +100,7 @@ static const BasisSet sb_even_basis = {
 typedef struct BrillDataContext {
     /* options */
     CCTK_REAL amplitude;
+    int       eppley_n;
 
     double (*q_func)(struct BrillDataContext *bd, double rho, double z);
     double (*laplace_q_func)(struct BrillDataContext *bd, double rho, double z);
@@ -337,7 +338,24 @@ static double laplace_q_gundlach(BrillDataContext *bd, double rho, double z)
     return 2 * bd->amplitude * e * (1.0 + 2 * r2 * (r2 + z2 - 3));
 }
 
-static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude,
+static double q_eppley(BrillDataContext *bd, double rho, double z)
+{
+    double r2 = SQR(rho) + SQR(z);
+    return bd->amplitude * SQR(rho) / (1.0 + pow(r2, bd->eppley_n / 2.0));
+}
+
+static double laplace_q_eppley(BrillDataContext *bd, double rho, double z)
+{
+    const int    n   = bd->eppley_n;
+    const double r2  = SQR(rho) + SQR(z);
+    const double r2n = pow(r2, n);
+    const double rn  = sqrt(r2n);
+    const double U   = 1.0 / (1 + rn);
+
+    return bd->amplitude * (2 * U - n * (n + 4) * SQR(rho) * SQR(U) * rn / r2 + 2 * SQR(n) * SQR(rho) * SQR(U) * U * r2n /r2);
+}
+
+static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude, int eppley_n,
                                  int basis_order_r, int basis_order_z,
                                  int colloc_offset_r, int colloc_offset_z,
                                  double sf, int overdet)
@@ -351,7 +369,11 @@ static BrillDataContext *init_bd(cGH *cctkGH, CCTK_REAL amplitude,
 
     bd->q_func         = q_gundlach;
     bd->laplace_q_func = laplace_q_gundlach;
+    //bd->q_func           = q_eppley;
+    //bd->laplace_q_func   = laplace_q_eppley;
+
     bd->amplitude      = amplitude;
+    bd->eppley_n       = eppley_n;
 
     bd->basis      = &sb_even_basis;
 
@@ -391,7 +413,7 @@ void brill_data(CCTK_ARGUMENTS)
 
     /* on the first run, solve the equation for ψ */
     if (!prev_bd) {
-        bd = init_bd(cctkGH, amplitude, basis_order_r, basis_order_z,
+        bd = init_bd(cctkGH, amplitude, eppley_n, basis_order_r, basis_order_z,
                      colloc_grid_offset_r, colloc_grid_offset_z,
                      scale_factor, overdet);