Use long double math for constructing ψ.

The extra precision is apparently needed for summing over the basis
functions.

1 files changed, 23 insertions, 23 deletions

diff --git a/src/brill.c b/src/brill.c
index 73b4cc6..be9c7f3 100644
--- a/src/brill.c
+++ b/src/brill.c
@@ -27,16 +27,16 @@
 /* a set of basis functions */
 typedef struct BasisSet {
     /* evaluate the idx-th basis function at the specified point*/
-    CCTK_REAL (*eval)      (CCTK_REAL coord, int idx);
+    long double (*eval)      (long double coord, int idx);
     /* evaluate the first derivative of the idx-th basis function at the specified point*/
-    CCTK_REAL (*eval_diff1)(CCTK_REAL coord, int idx);
+    long double (*eval_diff1)(long double coord, int idx);
     /* evaluate the second derivative of the idx-th basis function at the specified point*/
-    CCTK_REAL (*eval_diff2)(CCTK_REAL coord, int idx);
+    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)
      */
-    CCTK_REAL (*colloc_point)(int order, int idx);
+    long double (*colloc_point)(int order, int idx);
 } BasisSet;
 
 /*
@@ -47,40 +47,40 @@ typedef struct BasisSet {
 
 #define SCALE_FACTOR 1
 
-static CCTK_REAL sb_even_eval(CCTK_REAL coord, int idx)
+static long double sb_even_eval(long double coord, int idx)
 {
-    CCTK_REAL val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+    long double val = (coord == 0.0) ? M_PI_2 : atanl(SCALE_FACTOR / fabsl(coord));
 
     idx *= 2;   // even only
 
-    return sin((idx + 1) * val);
+    return sinl((idx + 1) * val);
 }
 
-static CCTK_REAL sb_even_eval_diff1(CCTK_REAL coord, int idx)
+static long double sb_even_eval_diff1(long double coord, int idx)
 {
-    CCTK_REAL val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / coord );
+    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) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+    return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cosl((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
 }
 
-static CCTK_REAL sb_even_eval_diff2(CCTK_REAL coord, int idx)
+static long double sb_even_eval_diff2(long double coord, int idx)
 {
-    CCTK_REAL val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / coord);
+    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 * coord * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+    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 CCTK_REAL sb_even_colloc_point(int order, int idx)
+static long double sb_even_colloc_point(int order, int idx)
 {
-    CCTK_REAL t;
+    long double t;
     order *= 2;
 
     t = (idx + 2) * M_PI / (order + 4);
-    return SCALE_FACTOR / tan(t);
+    return SCALE_FACTOR / tanl(t);
 }
 
 static const BasisSet sb_even_basis = {
@@ -241,7 +241,7 @@ void brill_data(CCTK_ARGUMENTS)
 
     BrillDataContext *bd;
 
-    CCTK_REAL *basis_val_r, *basis_val_z;
+    long double *basis_val_r, *basis_val_z;
 
     int64_t grid_size = CCTK_GFINDEX3D(cctkGH,
                                        cctk_lsh[0] - 1,
@@ -292,15 +292,15 @@ void brill_data(CCTK_ARGUMENTS)
         for (int j = 0; j < cctk_lsh[1]; j++)
             for (int i = 0; i < cctk_lsh[0]; i++) {
                 int    index = CCTK_GFINDEX3D(cctkGH, i, j, k);
-                CCTK_REAL xx = x[index], yy = y[index], zz = z[index];
+                long double xx = x[index], yy = y[index], zz = z[index];
 
-                CCTK_REAL r2 = SQR(xx) + SQR(yy);
-                CCTK_REAL r  = sqrt(r2);
+                long double r2 = SQR(xx) + SQR(yy);
+                long double r  = sqrtl(r2);
 
-                CCTK_REAL q   = bd->q_func(bd, r, zz);
-                CCTK_REAL e2q = exp(2 * q);
+                long double q   = bd->q_func(bd, r, zz);
+                long double e2q = expl(2 * q);
 
-                CCTK_REAL psi = 1.0, psi4;
+                long double psi = 1.0, psi4;
 
                 for (int n = 0; n < bd->nb_coeffs_z; n++)
                     for (int m = 0; m < bd->nb_coeffs_x; m++)