Split the basis sets definitions into a separate file.

4 files changed, 220 insertions, 208 deletions

diff --git a/src/basis.c b/src/basis.c
new file mode 100644
index 0000000..5ea043e
--- /dev/null
+++ b/src/basis.c
@@ -0,0 +1,189 @@
+
+#include <math.h>
+
+#include "maximal_slicing_axi.h"
+
+static double sb_even_eval(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx *= 2;   // even only
+
+    return sin((idx + 1) * val);
+}
+
+static double sb_even_eval_diff1(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx *= 2;   // even only
+
+    return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double sb_even_eval_diff2(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx *= 2;   // even only
+
+    return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabs(coord) * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double sb_even_colloc_point(int order, int idx)
+{
+    double t;
+
+    idx = order - idx - 1;
+    //order *= 2;
+
+    //t = (idx + 2) * M_PI / (order + 4);
+    t = (idx + 2) * M_PI / (2 * order + 2);
+    return SCALE_FACTOR / tan(t);
+}
+
+const BasisSet msa_sb_even_basis = {
+    .eval         = sb_even_eval,
+    .eval_diff1   = sb_even_eval_diff1,
+    .eval_diff2   = sb_even_eval_diff2,
+    .colloc_point = sb_even_colloc_point,
+};
+
+static double tb_even_eval(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx++;
+    idx *= 2;   // even only
+
+    return cos(idx * val) - 1.0;
+}
+
+static double tb_even_eval_diff1(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx++;
+    idx *= 2;   // even only
+
+    return SCALE_FACTOR * idx * SGN(coord) * sin(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double tb_even_eval_diff2(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+
+    idx++;
+    idx *= 2;   // even only
+
+    return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabs(coord) * sin(idx * val) + SCALE_FACTOR * idx * cos(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+}
+
+static double tb_even_colloc_point(int order, int idx)
+{
+    double t;
+
+    idx = order - idx - 1;
+    //order *= 2;
+
+    //t = (idx + 2) * M_PI / (order + 4);
+    t = (idx + 3) * M_PI / (2 * (order + 2));
+    return SCALE_FACTOR / tan(t);
+}
+
+const BasisSet msa_tb_even_basis = {
+    .eval         = tb_even_eval,
+    .eval_diff1   = tb_even_eval_diff1,
+    .eval_diff2   = tb_even_eval_diff2,
+    .colloc_point = tb_even_colloc_point,
+};
+
+static double full_eval(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+    int flag = idx & 1;
+
+    idx /= 2;
+
+    if (flag) return sin((2 * idx + 1) * 4 * val);
+    else      return cos((2 * idx + 2) * 4 * val) - 1;
+}
+
+static double full_eval_diff1(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+    int flag = idx & 1;
+
+    idx /= 2;
+
+    if (flag) {
+        idx = 2 * idx + 1;
+        return -4 * idx * SCALE_FACTOR * cos(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+    } else {
+        idx = 2 * (idx + 1);
+        return 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+    }
+}
+
+static double full_eval_diff2(double coord, int idx)
+{
+    double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+    int flag = idx & 1;
+
+    idx /= 2;
+
+    if (flag) {
+        idx = 2 * idx + 1;
+        return (16 * SQR(idx * SCALE_FACTOR) * cos(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sin(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+    } else {
+        idx = 2 * (idx + 1);
+        return (-16 * SQR(idx * SCALE_FACTOR) * sin(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cos(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+    }
+}
+
+static double full_colloc_point(int order, int idx)
+{
+    double t;
+
+    idx = order - idx - 1;
+
+    t = (idx + 0.5) * M_PI / (2 * order);
+
+    return SCALE_FACTOR / tan(t);
+
+}
+
+const BasisSet msa_full_basis = {
+    .eval         = full_eval,
+    .eval_diff1   = full_eval_diff1,
+    .eval_diff2   = full_eval_diff2,
+    .colloc_point = full_colloc_point,
+};
+
+static double cheb_eval(double coord, int idx)
+{
+    return cos(2 * idx * acos(coord / SCALE_FACTOR));
+}
+
+static double cheb_eval_diff1(double coord, int idx)
+{
+    return 2 * idx * sin(2 * idx * acos(coord / SCALE_FACTOR)) / sqrt(SQR(SCALE_FACTOR) - SQR(coord));
+}
+
+static double cheb_eval_diff2(double coord, int idx)
+{
+    double t = acos(coord / SCALE_FACTOR);
+    return 2 * idx * (cos(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sin(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5));
+}
+
+static double cheb_colloc_point(int order, int idx)
+{
+    return SCALE_FACTOR * cos((idx + 0.01) * M_PI / (4 * order + 4));
+}
+
+const BasisSet msa_cheb_basis = {
+    .eval       = cheb_eval,
+    .eval_diff1 = cheb_eval_diff1,
+    .eval_diff2 = cheb_eval_diff2,
+    .colloc_point = cheb_colloc_point,
+};
diff --git a/src/make.code.defn b/src/make.code.defn
index 27625ce..a5091a4 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn MaximalSlicingAxi
 
 # Source files in this directory
-SRCS = maximal_slicing_axi.c
+SRCS = basis.c maximal_slicing_axi.c
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/maximal_slicing_axi.c b/src/maximal_slicing_axi.c
index 7d3e6f4..ee39b3f 100644
--- a/src/maximal_slicing_axi.c
+++ b/src/maximal_slicing_axi.c
@@ -21,15 +21,12 @@
 #include "util_Table.h"
 #include "Slicing.h"
 
-double ms_scalarproduct_metric_avx(int offset_r, int offset_z, const double *coeffs,
-                                   const double *basis_r, const double *basis_z);
+#include "maximal_slicing_axi.h"
 
 #define ACC_TEST 0
 
 #define MAX(x, y) ((x) > (y) ? (x) : (y))
 #define MIN(x, y) ((x) > (y) ? (y) : (x))
-#define SQR(x) ((x) * (x))
-#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
 #define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr))
 
 /*
@@ -45,215 +42,14 @@ int64_t gettime(void)
     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 }
 
-/* a set of basis functions */
-typedef struct BasisSet {
-    /* evaluate the idx-th basis function at the specified point*/
-    long double (*eval)      (long double coord, int idx);
-    /* evaluate the first derivative of the idx-th basis function at the specified point*/
-    long double (*eval_diff1)(long double coord, int idx);
-    /* evaluate the second derivative of the idx-th basis function at the specified point*/
-    long double (*eval_diff2)(long double coord, int idx);
-    /**
-     * Get the idx-th collocation point for the specified order.
-     * idx runs from 0 to order - 1 (inclusive)
-     */
-    long double (*colloc_point)(int order, int idx);
-} BasisSet;
-
 /*
  * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
  * SB(x, n) = sin((n + 1) arccot(|x| / L))
  * They are symmetric wrt origin and decay as 1/x in infinity.
  */
 
-static CCTK_REAL scale_factor;
-
-#define SCALE_FACTOR scale_factor
-
-static long double sb_even_eval(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return sinl((idx + 1) * val);
-}
-
-static long double sb_even_eval_diff1(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cosl((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double sb_even_eval_diff2(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx *= 2;   // even only
-
-    return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabsl(coord) * cosl((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sinl((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double sb_even_colloc_point(int order, int idx)
-{
-    long double t;
-
-    idx = order - idx - 1;
-    //order *= 2;
-
-    //t = (idx + 2) * M_PI / (order + 4);
-    t = (idx + 2) * M_PI / (2 * order + 2);
-    return SCALE_FACTOR / tanl(t);
-}
+double scale_factor;
 
-static const BasisSet sb_even_basis = {
-    .eval         = sb_even_eval,
-    .eval_diff1   = sb_even_eval_diff1,
-    .eval_diff2   = sb_even_eval_diff2,
-    .colloc_point = sb_even_colloc_point,
-};
-
-static long double tb_even_eval(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx++;
-    idx *= 2;   // even only
-
-    return cosl(idx * val) - 1.0;
-}
-
-static long double tb_even_eval_diff1(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx++;
-    idx *= 2;   // even only
-
-    return SCALE_FACTOR * idx * SGN(coord) * sinl(idx * val) / (SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double tb_even_eval_diff2(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-
-    idx++;
-    idx *= 2;   // even only
-
-    return -SCALE_FACTOR * idx * SGN(coord) * (2 * fabsl(coord) * sinl(idx * val) + SCALE_FACTOR * idx * cosl(idx * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
-}
-
-static long double tb_even_colloc_point(int order, int idx)
-{
-    long double t;
-
-    idx = order - idx - 1;
-    //order *= 2;
-
-    //t = (idx + 2) * M_PI / (order + 4);
-    t = (idx + 3) * M_PI / (2 * (order + 2));
-    return SCALE_FACTOR / tanl(t);
-}
-
-static const BasisSet tb_even_basis = {
-    .eval         = tb_even_eval,
-    .eval_diff1   = tb_even_eval_diff1,
-    .eval_diff2   = tb_even_eval_diff2,
-    .colloc_point = tb_even_colloc_point,
-};
-
-static long double full_eval(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-    int flag = idx & 1;
-
-    idx /= 2;
-
-    if (flag) return sinl((2 * idx + 1) * 4 * val);
-    else      return cosl((2 * idx + 2) * 4 * val) - 1;
-}
-
-static long double full_eval_diff1(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-    int flag = idx & 1;
-
-    idx /= 2;
-
-    if (flag) {
-        idx = 2 * idx + 1;
-        return -4 * idx * SCALE_FACTOR * cosl(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
-    } else {
-        idx = 2 * (idx + 1);
-        return 4 * idx * SCALE_FACTOR * sinl(idx * 4 * val) / (SQR(SCALE_FACTOR) + SQR(coord));
-    }
-}
-
-static long double full_eval_diff2(long double coord, int idx)
-{
-    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
-    int flag = idx & 1;
-
-    idx /= 2;
-
-    if (flag) {
-        idx = 2 * idx + 1;
-        return (16 * SQR(idx * SCALE_FACTOR) * cosl(idx * 4 * val) - 4 * idx * SCALE_FACTOR * sinl(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
-    } else {
-        idx = 2 * (idx + 1);
-        return (-16 * SQR(idx * SCALE_FACTOR) * sinl(idx * 4 * val) - 4 * idx * SCALE_FACTOR * cosl(idx * 4 * val) * 2 * coord) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
-    }
-}
-
-static long double full_colloc_point(int order, int idx)
-{
-    long double t;
-
-    idx = order - idx - 1;
-
-    t = (idx + 0.5) * M_PI / (2 * order);
-
-    return SCALE_FACTOR / tanl(t);
-
-}
-
-static const BasisSet full_basis = {
-    .eval         = full_eval,
-    .eval_diff1   = full_eval_diff1,
-    .eval_diff2   = full_eval_diff2,
-    .colloc_point = full_colloc_point,
-};
-
-static long double cheb_eval(long double coord, int idx)
-{
-    return cosl(2 * idx * acosl(coord / SCALE_FACTOR));
-}
-
-static long double cheb_eval_diff1(long double coord, int idx)
-{
-    return 2 * idx * sinl(2 * idx * acosl(coord / SCALE_FACTOR)) / sqrtl(SQR(SCALE_FACTOR) - SQR(coord));
-}
-
-static long double cheb_eval_diff2(long double coord, int idx)
-{
-    long double t = acosl(coord / SCALE_FACTOR);
-    return 2 * idx * (cosl(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sinl(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5));
-}
-
-static long double cheb_colloc_point(int order, int idx)
-{
-    return SCALE_FACTOR * cosl((idx + 0.01) * M_PI / (4 * order + 4));
-}
-
-static const BasisSet cheb_basis = {
-    .eval       = cheb_eval,
-    .eval_diff1 = cheb_eval_diff1,
-    .eval_diff2 = cheb_eval_diff2,
-    .colloc_point = cheb_colloc_point,
-};
 
 /* indices (in our code, not cactus structs) of the grid functions which we'll need to
  * interpolate on the pseudospectral grid */
@@ -1163,7 +959,7 @@ static MaximalSlicingContext *init_ms(cGH *cctkGH,
 
     ms->gh = cctkGH;
 
-    ms->basis      = &sb_even_basis;
+    ms->basis      = &msa_sb_even_basis;
     //ms->basis      = &full_basis;
     //ms->basis      = &cheb_basis;
 
diff --git a/src/maximal_slicing_axi.h b/src/maximal_slicing_axi.h
new file mode 100644
index 0000000..b15420c
--- /dev/null
+++ b/src/maximal_slicing_axi.h
@@ -0,0 +1,27 @@
+
+#define SQR(x) ((x) * (x))
+#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
+
+#define SCALE_FACTOR scale_factor
+
+/* a set of basis functions */
+typedef struct BasisSet {
+    /* evaluate the idx-th basis function at the specified point*/
+    double (*eval)      (double coord, int idx);
+    /* evaluate the first derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff1)(double coord, int idx);
+    /* evaluate the second derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff2)(double coord, int idx);
+    /**
+     * Get the idx-th collocation point for the specified order.
+     * idx runs from 0 to order - 1 (inclusive)
+     */
+    double (*colloc_point)(int order, int idx);
+} BasisSet;
+
+extern const BasisSet msa_cheb_basis;
+extern const BasisSet msa_full_basis;
+extern const BasisSet msa_tb_even_basis;
+extern const BasisSet msa_sb_even_basis;
+
+extern double scale_factor;