trumpet: macrofication

1 files changed, 28 insertions, 19 deletions

diff --git a/src/trumpet.c b/src/trumpet.c
index 1598a0c..7b65822 100644
--- a/src/trumpet.c
+++ b/src/trumpet.c
@@ -9,9 +9,21 @@
 #include "cctk_Parameters.h"
 
 #define MAX(x,y) (x) > (y) ? (x) : (y)
-#define SQR(x) (x)*(x)
-#define C 3*sqrt(3)/4
-#define EPS 1E-10
+#define SQR(x) ((x)*(x))
+
+#define MASS 1
+/*
+ * the constant C from e.g. 0908.1063v4
+ */
+#define TRUMPET_CONST (3*sqrt(3)*SQR(MASS)/4)
+/*
+ * small number to avoid r=0 singularities
+ */
+#define EPS 1E-08
+/*
+ * isotropic/coordinate radius
+ */
+#define ISO_R(index) (sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS))
 
 static inline double r_areal_to_isotropic(double r, double m)
 {
@@ -40,7 +52,7 @@ static CCTK_REAL get_max_r(CCTK_INT dim)
 
 static void get_r_tables(double **pr_iso, double **pr_areal, CCTK_INT *size, CCTK_INT dim)
 {
-#define STEP 0.01
+#define STEP 0.0001
     double *r_iso, *r_areal, max = get_max_r(dim);
     int count = max*2/STEP, i = 0;
 
@@ -90,18 +102,16 @@ void trumpet_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 r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]));
+                CCTK_REAL r = ISO_R(index);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
-                if (r < EPS)
-                    r = EPS;
                 CCTK_REAL psi2 = R/r, psi4 = psi2*psi2;
-                CCTK_REAL k_fact = C/(r*r*r*r*R);
+                CCTK_REAL k_fact = TRUMPET_CONST/(r*r*r*r*R);
 
                 gxx[index] = gyy[index] = gzz[index] = psi4;
 
-                kxx[index] = -k_fact*(3*SQR(x[index]) - r*r);
-                kyy[index] = -k_fact*(3*SQR(y[index]) - r*r);
-                kzz[index] = -k_fact*(3*SQR(z[index]) - r*r);
+                kxx[index] = -k_fact*(3*SQR(x[index]) - SQR(r));
+                kyy[index] = -k_fact*(3*SQR(y[index]) - SQR(r));
+                kzz[index] = -k_fact*(3*SQR(z[index]) - SQR(r));
 
                 kxy[index] = -k_fact*3*x[index]*y[index];
                 kxz[index] = -k_fact*3*x[index]*z[index];
@@ -133,10 +143,8 @@ void trumpet_lapse(CCTK_ARGUMENTS)
                 int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
                 CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
-                if (r < EPS)
-                    r = EPS;
 
-                alp[index] = sqrt(1 - 2*1/R + C*C/(R*R*R*R));
+                alp[index] = sqrt(1 - 2*MASS/R + TRUMPET_CONST*TRUMPET_CONST/(R*R*R*R));
             }
 }
 
@@ -164,11 +172,12 @@ void trumpet_shift(CCTK_ARGUMENTS)
                 int   index = CCTK_GFINDEX3D(cctkGH, i, j, k);
                 CCTK_REAL r = sqrt(SQR(x[index]) + SQR(y[index]) + SQR(z[index]) + EPS);
                 CCTK_REAL R = gsl_spline_eval(spline, r, acc);
-                if (r < EPS)
-                    r = EPS;
 
-                betax[index] = C*x[index]/(R*R*R);
-                betay[index] = C*y[index]/(R*R*R);
-                betaz[index] = C*z[index]/(R*R*R);
+                betax[index] = TRUMPET_CONST*x[index]/(R*R*R);
+                betay[index] = TRUMPET_CONST*y[index]/(R*R*R);
+                betaz[index] = TRUMPET_CONST*z[index]/(R*R*R);
+            }
+}
+
             }
 }