Use CCTK_REAL instead of double. This allows using higher precisions

that double.

Do not initialise the ghost zones; synchronise instead.  Since
interpolating the initial data is expensive this should save some
time.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@57 b2a53a04-0f4f-0410-87ed-f9f25ced00cf

8 files changed, 287 insertions, 278 deletions

diff --git a/src/CoordTransf.c b/src/CoordTransf.c
index 22291af..3a8adb7 100644
--- a/src/CoordTransf.c
+++ b/src/CoordTransf.c
@@ -12,7 +12,7 @@
 
 /*-----------------------------------------------------------*/
 void
-AB_To_XR (int nvar, double A, double B, double *X, double *R,
+AB_To_XR (int nvar, CCTK_REAL A, CCTK_REAL B, CCTK_REAL *X, CCTK_REAL *R,
 	  derivs U)
 /* On Entrance: U.d0[]=U[]; U.d1[] =U[]_A;  U.d2[] =U[]_B;  U.d3[] =U[]_3;  */
 /*                          U.d11[]=U[]_AA; U.d12[]=U[]_AB; U.d13[]=U[]_A3; */
@@ -22,7 +22,7 @@ AB_To_XR (int nvar, double A, double B, double *X, double *R,
 /*                          U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33; */
 {
   DECLARE_CCTK_PARAMETERS;
-  double At = 0.5 * (A + 1), A_X, A_XX, B_R, B_RR;
+  CCTK_REAL At = 0.5 * (A + 1), A_X, A_XX, B_R, B_RR;
   int ivar;
 
   *X = 2 * atanh (At);
@@ -47,7 +47,7 @@ AB_To_XR (int nvar, double A, double B, double *X, double *R,
 
 /*-----------------------------------------------------------*/
 void
-C_To_c (int nvar, double X, double R, double *x, double *r,
+C_To_c (int nvar, CCTK_REAL X, CCTK_REAL R, CCTK_REAL *x, CCTK_REAL *r,
 	derivs U)
 /* On Entrance: U.d0[]=U[]; U.d1[] =U[]_X;  U.d2[] =U[]_R;  U.d3[] =U[]_3;  */
 /*                          U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3; */
@@ -57,7 +57,7 @@ C_To_c (int nvar, double X, double R, double *x, double *r,
 /*                          U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; */
 {
   DECLARE_CCTK_PARAMETERS;
-  double C_c2, U_cb, U_CB;
+  CCTK_REAL C_c2, U_cb, U_CB;
   dcomplex C, C_c, C_cc, c, c_C, c_CC, U_c, U_cc, U_C, U_CC, One =
     Complex (1., 0.);
   int ivar;
@@ -115,8 +115,8 @@ C_To_c (int nvar, double X, double R, double *x, double *r,
 
 /*-----------------------------------------------------------*/
 void
-rx3_To_xyz (int nvar, double x, double r, double phi,
-	    double *y, double *z, derivs U)
+rx3_To_xyz (int nvar, CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi,
+	    CCTK_REAL *y, CCTK_REAL *z, derivs U)
 /* On Entrance: U.d0[]=U[]; U.d1[] =U[]_x;  U.d2[] =U[]_r;  U.d3[] =U[]_3;  */
 /*                          U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3; */
 /*                          U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33; */
@@ -125,7 +125,7 @@ rx3_To_xyz (int nvar, double x, double r, double phi,
 /*                          U.d22[]=U[]_yy; U.d2z[]=U[]_yz; U.dzz[]=U[]_zz; */
 {
   int jvar;
-  double
+  CCTK_REAL
     sin_phi = sin (phi),
     cos_phi = cos (phi),
     sin2_phi = sin_phi * sin_phi,
@@ -138,7 +138,7 @@ rx3_To_xyz (int nvar, double x, double r, double phi,
 
   for (jvar = 0; jvar < nvar; jvar++)
   {
-    double U_x = U.d1[jvar], U_r = U.d2[jvar], U_3 = U.d3[jvar],
+    CCTK_REAL U_x = U.d1[jvar], U_r = U.d2[jvar], U_3 = U.d3[jvar],
       U_xx = U.d11[jvar], U_xr = U.d12[jvar], U_x3 = U.d13[jvar],
       U_rr = U.d22[jvar], U_r3 = U.d23[jvar], U_33 = U.d33[jvar];
     U.d1[jvar] = U_x;		/* U_x*/
diff --git a/src/Equations.c b/src/Equations.c
index 42a2a08..1fd8876 100644
--- a/src/Equations.c
+++ b/src/Equations.c
@@ -20,12 +20,12 @@
 /* U.d23[ivar]  = U[ivar]_yz;*/
 /* U.d33[ivar]  = U[ivar]_zz;*/
 
-double
-BY_KKofxyz (double x, double y, double z)
+CCTK_REAL
+BY_KKofxyz (CCTK_REAL x, CCTK_REAL y, CCTK_REAL z)
 {
   DECLARE_CCTK_PARAMETERS;
   int i, j;
-  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
+  CCTK_REAL r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
     Aij, AijAij, n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
 
   r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
@@ -79,11 +79,11 @@ BY_KKofxyz (double x, double y, double z)
 }
 
 void
-BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
+BY_Aijofxyz (CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL Aij[3][3])
 {
   DECLARE_CCTK_PARAMETERS;
   int i, j;
-  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
+  CCTK_REAL r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
     n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
 
   r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
@@ -143,13 +143,13 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
 /*-----------------------------------------------------------*/
 void
 NonLinEquations (CCTK_REAL rho_adm,
-     double A, double B, double X, double R,
-		 double x, double r, double phi,
-		 double y, double z, derivs U, double *values)
+     CCTK_REAL A, CCTK_REAL B, CCTK_REAL X, CCTK_REAL R,
+		 CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi,
+		 CCTK_REAL y, CCTK_REAL z, derivs U, CCTK_REAL *values)
 {
   DECLARE_CCTK_PARAMETERS;
-  double r_plus, r_minus, psi, psi2, psi4, psi7;
-  double mu;
+  CCTK_REAL r_plus, r_minus, psi, psi2, psi4, psi7;
+  CCTK_REAL mu;
 
   r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
   r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
@@ -170,12 +170,12 @@ NonLinEquations (CCTK_REAL rho_adm,
 /********               Linear Equations                ***********/
 /*-----------------------------------------------------------*/
 void
-LinEquations (double A, double B, double X, double R,
-	      double x, double r, double phi,
-	      double y, double z, derivs dU, derivs U, double *values)
+LinEquations (CCTK_REAL A, CCTK_REAL B, CCTK_REAL X, CCTK_REAL R,
+	      CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi,
+	      CCTK_REAL y, CCTK_REAL z, derivs dU, derivs U, CCTK_REAL *values)
 {
   DECLARE_CCTK_PARAMETERS;
-  double r_plus, r_minus, psi, psi2, psi4, psi8;
+  CCTK_REAL r_plus, r_minus, psi, psi2, psi4, psi8;
 
   r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
   r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index d2c41a6..7a6e13b 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -16,7 +16,7 @@
 /*#define FAC sin(al)*sin(be)*sin(al)*sin(be)*/
 /*#define FAC 1*/
 
-static inline double min (double const x, double const y)
+static inline CCTK_REAL min (CCTK_REAL const x, CCTK_REAL const y)
 {
   return x<y ? x : y;
 }
@@ -84,7 +84,7 @@ void
 Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
 {
   int i, j, k, ivar, N, *indx;
-  double *p, *dp, *d2p, *q, *dq, *r, *dr;
+  CCTK_REAL *p, *dp, *d2p, *q, *dq, *r, *dr;
 
   N = maximum3 (n1, n2, n3);
   p = dvector (0, N);
@@ -190,7 +190,7 @@ Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
 /* --------------------------------------------------------------------------*/
 void
 F_of_v (CCTK_POINTER_TO_CONST cctkGH,
-        int nvar, int n1, int n2, int n3, derivs v, double *F,
+        int nvar, int n1, int n2, int n3, derivs v, CCTK_REAL *F,
         derivs u)
 {
   /*      Calculates the left hand sides of the non-linear equations F_m(v_n)=0*/
@@ -199,7 +199,7 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
   /*      at interior points and at the boundaries "+/-"*/
   DECLARE_CCTK_PARAMETERS;
   int i, j, k, ivar, indx;
-  double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+  CCTK_REAL al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
   derivs U;
   CCTK_REAL *sources;
 
@@ -264,7 +264,7 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
 
   Derivatives_AB3 (nvar, n1, n2, n3, v);
   FILE *debugfile;
-  double psi, psi2, psi4, psi7, r_plus, r_minus;
+  CCTK_REAL psi, psi2, psi4, psi7, r_plus, r_minus;
   if (do_residuum_debug_output)
   {
     debugfile=fopen("res.dat", "w");
@@ -338,18 +338,18 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
           psi7 = psi * psi2 * psi4;
           fprintf(debugfile,
                   "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
-             x, y, A, B,
-             (U.d11[0] +
-              U.d22[0] +
-              U.d33[0] +
-/*              0.125 * BY_KKofxyz (x, y, z) / psi7 +*/
-              2.0 * Pi / psi2/psi * sources[indx]) * FAC,
-             (U.d11[0] +
-              U.d22[0] +
-              U.d33[0])*FAC,
-             -(2.0 * Pi / psi2/psi * sources[indx]) * FAC,
-             sources[indx]
-             /*F[indx]*/
+             (double)x, (double)y, (double)A, (double)B,
+             (double)(U.d11[0] +
+                      U.d22[0] +
+                      U.d33[0] +
+/*                      0.125 * BY_KKofxyz (x, y, z) / psi7 +*/
+                      (2.0 * Pi / psi2/psi * sources[indx]) * FAC),
+             (double)((U.d11[0] +
+                       U.d22[0] +
+                       U.d33[0])*FAC),
+             (double)(-(2.0 * Pi / psi2/psi * sources[indx]) * FAC),
+             (double)sources[indx]
+             /*(double)F[indx]*/
              );
         }
       }
@@ -367,14 +367,14 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
 /* --------------------------------------------------------------------------*/
 void
 J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
-	    double *Jdv, derivs u)
+	    CCTK_REAL *Jdv, derivs u)
 {				/*      Calculates the left hand sides of the non-linear equations F_m(v_n)=0*/
   /*      and the function u (u.d0[]) as well as its derivatives*/
   /*      (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])*/
   /*      at interior points and at the boundaries "+/-"*/
   DECLARE_CCTK_PARAMETERS;
   int i, j, k, ivar, indx;
-  double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+  CCTK_REAL al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
   derivs dU, U;
 
 
@@ -448,17 +448,17 @@ J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
 /* --------------------------------------------------------------------------*/
 void
 JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2,
-	      int n3, derivs dv, derivs u, double *values)
+	      int n3, derivs dv, derivs u, CCTK_REAL *values)
 {				/* Calculates rows of the vector 'J(FD)*dv'.*/
   /* First row to be calculated: row = Index(0,      i, j, k; nvar, n1, n2, n3)*/
   /* Last  row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3)*/
   /* These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1].*/
   DECLARE_CCTK_PARAMETERS;
   int ivar, indx;
-  double al, be, A, B, X, R, x, r, phi, y, z, Am1;
-  double sin_al, sin_al_i1, sin_al_i2, sin_al_i3, cos_al;
-  double sin_be, sin_be_i1, sin_be_i2, sin_be_i3, cos_be;
-  double dV0, dV1, dV2, dV3, dV11, dV12, dV13, dV22, dV23, dV33,
+  CCTK_REAL al, be, A, B, X, R, x, r, phi, y, z, Am1;
+  CCTK_REAL sin_al, sin_al_i1, sin_al_i2, sin_al_i3, cos_al;
+  CCTK_REAL sin_be, sin_be_i1, sin_be_i2, sin_be_i3, cos_be;
+  CCTK_REAL dV0, dV1, dV2, dV3, dV11, dV12, dV13, dV22, dV23, dV33,
     ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
   derivs dU, U;
 
@@ -590,13 +590,13 @@ JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2,
 /* --------------------------------------------------------------------------*/
 void
 SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u,
-	       int *ncols, int **cols, double **Matrix)
+	       int *ncols, int **cols, CCTK_REAL **Matrix)
 {
   DECLARE_CCTK_PARAMETERS;
   int column, row, mcol;
   int i, i1, i_0, i_1, j, j1, j_0, j_1, k, k1, k_0, k_1, N1, N2, N3,
     ivar, ivar1, ntotal = nvar * n1 * n2 * n3;
-  double *values;
+  CCTK_REAL *values;
   derivs dv;
 
   values = dvector (0, nvar - 1);
@@ -679,12 +679,12 @@ SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u,
 
 /* --------------------------------------------------------------------------*/
 /* Calculates the value of v at an arbitrary position (A,B,phi)*/
-double
-PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
+CCTK_REAL
+PunctEvalAtArbitPosition (CCTK_REAL *v, CCTK_REAL A, CCTK_REAL B, CCTK_REAL phi,
 			  int n1, int n2, int n3)
 {
   int i, j, k, N;
-  double *p, *values1, **values2, result;
+  CCTK_REAL *p, *values1, **values2, result;
 
   N = maximum3 (n1, n2, n3);
   p = dvector (0, N);
@@ -725,7 +725,7 @@ void
 calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2,
 		  int n3, derivs v, derivs vv)
 {
-  double al = Pih * (2 * i + 1) / n1, be = Pih * (2 * j + 1) / n2,
+  CCTK_REAL al = Pih * (2 * i + 1) / n1, be = Pih * (2 * j + 1) / n2,
     sin_al = sin (al), sin2_al = sin_al * sin_al, cos_al = cos (al),
     sin_be = sin (be), sin2_be = sin_be * sin_be, cos_be = cos (be);
 
@@ -745,8 +745,8 @@ calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2,
 }
 
 /* --------------------------------------------------------------------------*/
-double
-interpol (double a, double b, double c, derivs v)
+CCTK_REAL
+interpol (CCTK_REAL a, CCTK_REAL b, CCTK_REAL c, derivs v)
 {
   return v.d0[0]
     + a * v.d1[0] + b * v.d2[0] + c * v.d3[0]
@@ -756,13 +756,13 @@ interpol (double a, double b, double c, derivs v)
 
 /* --------------------------------------------------------------------------*/
 /* Calculates the value of v at an arbitrary position (x,y,z)*/
-double
+CCTK_REAL
 PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
-                                  int n2, int n3, derivs v, double x, double y,
-                                  double z)
+                                  int n2, int n3, derivs v, CCTK_REAL x, CCTK_REAL y,
+                                  CCTK_REAL z)
 {
   DECLARE_CCTK_PARAMETERS;
-  double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
+  CCTK_REAL xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
     result, Ui;
   int i, j, k;
   derivs vv;
@@ -807,13 +807,13 @@ PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
 
 /* --------------------------------------------------------------------------*/
 /* Calculates the value of v at an arbitrary position (x,y,z)*/
-double
+CCTK_REAL
 PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
-			    int n2, int n3, derivs v, double x, double y,
-			    double z)
+			    int n2, int n3, derivs v, CCTK_REAL x, CCTK_REAL y,
+			    CCTK_REAL z)
 {
   DECLARE_CCTK_PARAMETERS;
-  double xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
+  CCTK_REAL xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
 
   xs = x / par_b;
   ys = y / par_b;
diff --git a/src/Newton.c b/src/Newton.c
index 7225c87..058a443 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -12,26 +12,26 @@
 static int
 bicgstab (CCTK_POINTER_TO_CONST cctkGH,
           int nvar, int n1, int n2, int n3, derivs v,
-          derivs dv, int output, int itmax, double tol, double *normres);
+          derivs dv, int output, int itmax, CCTK_REAL tol, CCTK_REAL *normres);
 static void
-relax (double *dv, int nvar, int n1, int n2, int n3,
-       double *rhs, int *ncols, int **cols, double **JFD);
+relax (CCTK_REAL *dv, int nvar, int n1, int n2, int n3,
+       CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
 static void
-resid (double *res, int ntotal, double *dv, double *rhs,
-       int *ncols, int **cols, double **JFD);
+resid (CCTK_REAL *res, int ntotal, CCTK_REAL *dv, CCTK_REAL *rhs,
+       int *ncols, int **cols, CCTK_REAL **JFD);
 static void
-LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3,
-	      double *rhs, int *ncols, int **cols, double **JFD);
+LineRelax_al (CCTK_REAL *dv, int j, int k, int nvar, int n1, int n2, int n3,
+	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
 static void
-LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
-	      double *rhs, int *ncols, int **cols, double **JFD);
+LineRelax_be (CCTK_REAL *dv, int i, int k, int nvar, int n1, int n2, int n3,
+	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
 /* --------------------------------------------------------------------------*/
 static void
-resid (double *res, int ntotal, double *dv, double *rhs,
-       int *ncols, int **cols, double **JFD)
+resid (CCTK_REAL *res, int ntotal, CCTK_REAL *dv, CCTK_REAL *rhs,
+       int *ncols, int **cols, CCTK_REAL **JFD)
 {
   int i, m;
-  double JFDdv_i;
+  CCTK_REAL JFDdv_i;
 
   for (i = 0; i < ntotal; i++)
   {
@@ -44,11 +44,11 @@ resid (double *res, int ntotal, double *dv, double *rhs,
 
 /* -------------------------------------------------------------------------*/
 static void
-LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3,
-	      double *rhs, int *ncols, int **cols, double **JFD)
+LineRelax_al (CCTK_REAL *dv, int j, int k, int nvar, int n1, int n2, int n3,
+	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
 {
   int i, m, Ic, Ip, Im, col, ivar;
-  double *a, *b, *c, *r, *u;
+  CCTK_REAL *a, *b, *c, *r, *u;
   a = dvector (0, n1 - 1);
   b = dvector (0, n1 - 1);
   c = dvector (0, n1 - 1);
@@ -98,11 +98,11 @@ LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3,
 
 /* --------------------------------------------------------------------------*/
 static void
-LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
-	      double *rhs, int *ncols, int **cols, double **JFD)
+LineRelax_be (CCTK_REAL *dv, int i, int k, int nvar, int n1, int n2, int n3,
+	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
 {
   int j, m, Ic, Ip, Im, col, ivar;
-  double *a, *b, *c, *r, *u;
+  CCTK_REAL *a, *b, *c, *r, *u;
   a = dvector (0, n2 - 1);
   b = dvector (0, n2 - 1);
   c = dvector (0, n2 - 1);
@@ -152,8 +152,8 @@ LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
 
 /* --------------------------------------------------------------------------*/
 static void
-relax (double *dv, int nvar, int n1, int n2, int n3,
-       double *rhs, int *ncols, int **cols, double **JFD)
+relax (CCTK_REAL *dv, int nvar, int n1, int n2, int n3,
+       CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
 {
   int i, j, k, n;
 
@@ -191,12 +191,12 @@ relax (double *dv, int nvar, int n1, int n2, int n3,
 void
 TestRelax (CCTK_POINTER_TO_CONST cctkGH,
      int nvar, int n1, int n2, int n3, derivs v,
-	   double *dv)
+	   CCTK_REAL *dv)
 {
   DECLARE_CCTK_PARAMETERS;
   int ntotal = n1 * n2 * n3 * nvar, **cols, *ncols, maxcol =
     StencilSize * nvar, j;
-  double *F, *res, **JFD;
+  CCTK_REAL *F, *res, **JFD;
   derivs u;
 
   F = dvector (0, ntotal - 1);
@@ -214,7 +214,7 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH,
   for (j = 0; j < ntotal; j++)
     dv[j] = 0;
   resid (res, ntotal, dv, F, ncols, cols, JFD);
-  printf ("Before: |F|=%20.15e\n", norm1 (res, ntotal));
+  printf ("Before: |F|=%20.15e\n", (double) norm1 (res, ntotal));
   fflush(stdout);
   for (j = 0; j < NRELAX; j++)
   {
@@ -222,13 +222,13 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH,
     if (j % Step_Relax == 0)
     {
       resid (res, ntotal, dv, F, ncols, cols, JFD);
-      printf ("j=%d\t |F|=%20.15e\n", j, norm1 (res, ntotal));
+      printf ("j=%d\t |F|=%20.15e\n", j, (double) norm1 (res, ntotal));
       fflush(stdout);
     }
   }
 
   resid (res, ntotal, dv, F, ncols, cols, JFD);
-  printf ("After: |F|=%20.15e\n", norm1 (res, ntotal));
+  printf ("After: |F|=%20.15e\n", (double) norm1 (res, ntotal));
   fflush(stdout);
 
   free_dvector (F, 0, ntotal - 1);
@@ -244,17 +244,17 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH,
 static int
 bicgstab (CCTK_POINTER_TO_CONST cctkGH,
           int nvar, int n1, int n2, int n3, derivs v,
-	        derivs dv, int output, int itmax, double tol, double *normres)
+	        derivs dv, int output, int itmax, CCTK_REAL tol, CCTK_REAL *normres)
 {
   DECLARE_CCTK_PARAMETERS;
   int ntotal = n1 * n2 * n3 * nvar, ii, j;
-  double alpha = 0, beta = 0;
-  double rho = 0, rho1 = 1, rhotol = 1e-50;
-  double omega = 0, omegatol = 1e-50;
-  double *p, *rt, *s, *t, *r, *vv;
-  double **JFD;
+  CCTK_REAL alpha = 0, beta = 0;
+  CCTK_REAL rho = 0, rho1 = 1, rhotol = 1e-50;
+  CCTK_REAL omega = 0, omegatol = 1e-50;
+  CCTK_REAL *p, *rt, *s, *t, *r, *vv;
+  CCTK_REAL **JFD;
   int **cols, *ncols, maxcol = StencilSize * nvar;
-  double *F;
+  CCTK_REAL *F;
   derivs u, ph, sh;
 
   F = dvector (0, ntotal - 1);
@@ -281,7 +281,7 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 
   /* check */
   if (output == 1) {
-    printf ("bicgstab:  itmax %d, tol %e\n", itmax, tol);
+    printf ("bicgstab:  itmax %d, tol %e\n", itmax, (double)tol);
     fflush(stdout);
   }
 
@@ -292,7 +292,7 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 
   *normres = norm2 (r, ntotal);
   if (output == 1) {
-    printf ("bicgstab: %5d  %10.3e\n", 0, *normres);
+    printf ("bicgstab: %5d  %10.3e\n", 0, (double) *normres);
     fflush(stdout);
   }
 
@@ -336,7 +336,7 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 	dv.d0[j] += alpha * ph.d0[j];
       if (output == 1) {
 	printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
-		ii + 1, *normres, alpha, beta, omega);
+		ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
         fflush(stdout);
       }
       break;
@@ -361,7 +361,7 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
     *normres = norm2 (r, ntotal);
     if (output == 1) {
       printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
-	      ii + 1, *normres, alpha, beta, omega);
+	      ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
       fflush(stdout);
     }
     if (*normres <= tol)
@@ -409,11 +409,11 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 void
 Newton (CCTK_POINTER_TO_CONST cctkGH,
   int nvar, int n1, int n2, int n3,
-	derivs v, double tol, int itmax)
+	derivs v, CCTK_REAL tol, int itmax)
 {
   DECLARE_CCTK_PARAMETERS;
   int ntotal = n1 * n2 * n3 * nvar, ii, j, it;
-  double *F, dmax, normres;
+  CCTK_REAL *F, dmax, normres;
   derivs u, dv;
 
   F = dvector (0, ntotal - 1);
@@ -435,7 +435,7 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
     }
     for (j = 0; j < ntotal; j++)
       dv.d0[j] = 0;
-    printf ("Newton: it=%d \t |F|=%e\n", it, dmax);
+    printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
     fflush(stdout);
     ii =
       bicgstab (cctkGH,
@@ -461,7 +461,7 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
         if (fabs (F[j]) > dmax)
           dmax = fabs (F[j]);
   }
-  printf ("Newton: it=%d \t |F|=%e \n", it, dmax);
+  printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
   fflush(stdout);
 
   free_dvector (F, 0, ntotal - 1);
diff --git a/src/TP_utilities.c b/src/TP_utilities.c
index 1736f7f..2d44511 100644
--- a/src/TP_utilities.c
+++ b/src/TP_utilities.c
@@ -31,13 +31,13 @@ ivector (long nl, long nh)
 }
 
 /*---------------------------------------------------------------------------*/
-double *
+CCTK_REAL *
 dvector (long nl, long nh)
-/* allocate a double vector with subscript range v[nl..nh] */
+/* allocate a CCTK_REAL vector with subscript range v[nl..nh] */
 {
-  double *v;
+  CCTK_REAL *v;
 
-  v = (double *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (double)));
+  v = (CCTK_REAL *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (CCTK_REAL)));
   if (!v)
     nrerror ("allocation failure in dvector()");
   return v - nl + NR_END;
@@ -74,15 +74,15 @@ imatrix (long nrl, long nrh, long ncl, long nch)
 }
 
 /*---------------------------------------------------------------------------*/
-double **
+CCTK_REAL **
 dmatrix (long nrl, long nrh, long ncl, long nch)
-/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+/* allocate a CCTK_REAL matrix with subscript range m[nrl..nrh][ncl..nch] */
 {
   long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
-  double **m;
+  CCTK_REAL **m;
 
   /* allocate pointers to rows */
-  m = (double **) malloc ((size_t) ((nrow + NR_END) * sizeof (double *)));
+  m = (CCTK_REAL **) malloc ((size_t) ((nrow + NR_END) * sizeof (CCTK_REAL *)));
   if (!m)
     nrerror ("allocation failure 1 in matrix()");
   m += NR_END;
@@ -90,7 +90,7 @@ dmatrix (long nrl, long nrh, long ncl, long nch)
 
   /* allocate rows and set pointers to them */
   m[nrl] =
-    (double *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (double)));
+    (CCTK_REAL *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (CCTK_REAL)));
   if (!m[nrl])
     nrerror ("allocation failure 2 in matrix()");
   m[nrl] += NR_END;
@@ -104,15 +104,15 @@ dmatrix (long nrl, long nrh, long ncl, long nch)
 }
 
 /*---------------------------------------------------------------------------*/
-double ***
+CCTK_REAL ***
 d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
-/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+/* allocate a CCTK_REAL 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
 {
   long i, j, nrow = nrh - nrl + 1, ncol = nch - ncl + 1, ndep = ndh - ndl + 1;
-  double ***t;
+  CCTK_REAL ***t;
 
   /* allocate pointers to pointers to rows */
-  t = (double ***) malloc ((size_t) ((nrow + NR_END) * sizeof (double **)));
+  t = (CCTK_REAL ***) malloc ((size_t) ((nrow + NR_END) * sizeof (CCTK_REAL **)));
   if (!t)
     nrerror ("allocation failure 1 in f3tensor()");
   t += NR_END;
@@ -120,8 +120,8 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
 
   /* allocate pointers to rows and set pointers to them */
   t[nrl] =
-    (double **)
-    malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (double *)));
+    (CCTK_REAL **)
+    malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (CCTK_REAL *)));
   if (!t[nrl])
     nrerror ("allocation failure 2 in f3tensor()");
   t[nrl] += NR_END;
@@ -129,8 +129,8 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
 
   /* allocate rows and set pointers to them */
   t[nrl][ncl] =
-    (double *)
-    malloc ((size_t) ((nrow * ncol * ndep + NR_END) * sizeof (double)));
+    (CCTK_REAL *)
+    malloc ((size_t) ((nrow * ncol * ndep + NR_END) * sizeof (CCTK_REAL)));
   if (!t[nrl][ncl])
     nrerror ("allocation failure 3 in f3tensor()");
   t[nrl][ncl] += NR_END;
@@ -160,8 +160,8 @@ free_ivector (int *v, long nl, long nh)
 
 /*--------------------------------------------------------------------------*/
 void
-free_dvector (double *v, long nl, long nh)
-/* free a double vector allocated with dvector() */
+free_dvector (CCTK_REAL *v, long nl, long nh)
+/* free a CCTK_REAL vector allocated with dvector() */
 {
   free ((FREE_ARG) (v + nl - NR_END));
 }
@@ -177,8 +177,8 @@ free_imatrix (int **m, long nrl, long nrh, long ncl, long nch)
 
 /*--------------------------------------------------------------------------*/
 void
-free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch)
-/* free a double matrix allocated by dmatrix() */
+free_dmatrix (CCTK_REAL **m, long nrl, long nrh, long ncl, long nch)
+/* free a CCTK_REAL matrix allocated by dmatrix() */
 {
   free ((FREE_ARG) (m[nrl] + ncl - NR_END));
   free ((FREE_ARG) (m + nrl - NR_END));
@@ -186,9 +186,9 @@ free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch)
 
 /*--------------------------------------------------------------------------*/
 void
-free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch,
+free_d3tensor (CCTK_REAL ***t, long nrl, long nrh, long ncl, long nch,
 	       long ndl, long ndh)
-/* free a double f3tensor allocated by f3tensor() */
+/* free a CCTK_REAL f3tensor allocated by f3tensor() */
 {
   free ((FREE_ARG) (t[nrl][ncl] + ndl - NR_END));
   free ((FREE_ARG) (t[nrl] + ncl - NR_END));
@@ -253,22 +253,22 @@ pow_int (int mantisse, int exponent)
 
 /*--------------------------------------------------------------------------*/
 #if 0
-double
-atanh (double x)
+CCTK_REAL
+atanh (CCTK_REAL x)
 {
   return 0.5 * log ((1 + x) / (1 - x));
 }
 
 /*--------------------------------------------------------------------------*/
-double
-asinh (double x)
+CCTK_REAL
+asinh (CCTK_REAL x)
 {
   return log (x + sqrt (1 + x * x));
 }
 
 /*--------------------------------------------------------------------------*/
-double
-acosh (double x)
+CCTK_REAL
+acosh (CCTK_REAL x)
 {
   return log (x + sqrt (x * x - 1));
 }
@@ -306,7 +306,7 @@ Cmul (dcomplex a, dcomplex b)
 
 /*--------------------------------------------------------------------------*/
 dcomplex
-RCmul (double x, dcomplex a)
+RCmul (CCTK_REAL x, dcomplex a)
 {
   dcomplex c;
   c.r = x * a.r;
@@ -319,7 +319,7 @@ dcomplex
 Cdiv (dcomplex a, dcomplex b)
 {
   dcomplex c;
-  double r, den;
+  CCTK_REAL r, den;
   if (fabs (b.r) >= fabs (b.i))
   {
     r = b.i / b.r;
@@ -339,7 +339,7 @@ Cdiv (dcomplex a, dcomplex b)
 
 /*--------------------------------------------------------------------------*/
 dcomplex
-Complex (double re, double im)
+Complex (CCTK_REAL re, CCTK_REAL im)
 {
   dcomplex c;
   c.r = re;
@@ -358,10 +358,10 @@ Conjg (dcomplex z)
 }
 
 /*--------------------------------------------------------------------------*/
-double
+CCTK_REAL
 Cabs (dcomplex z)
 {
-  double x, y, ans, temp;
+  CCTK_REAL x, y, ans, temp;
   x = fabs (z.r);
   y = fabs (z.i);
   if (x == 0.0)
@@ -386,7 +386,7 @@ dcomplex
 Csqrt (dcomplex z)
 {
   dcomplex c;
-  double x, y, w, r;
+  CCTK_REAL x, y, w, r;
   if ((z.r == 0.0) && (z.i == 0.0))
   {
     c.r = 0.0;
@@ -426,7 +426,7 @@ dcomplex
 Cexp (dcomplex z)
 {
   dcomplex c;
-  double exp_r = exp (z.r);
+  CCTK_REAL exp_r = exp (z.r);
 
   c.r = exp_r * cos (z.i);
   c.i = exp_r * sin (z.i);
@@ -524,10 +524,10 @@ Ccoth (dcomplex z)
 
 /*--------------------------------------------------------------------------*/
 void
-chebft_Zeros (double u[], int n, int inv)
+chebft_Zeros (CCTK_REAL u[], int n, int inv)
 {
   int k, j, isignum;
-  double fac, sum, Pion, *c;
+  CCTK_REAL fac, sum, Pion, *c;
 
   c = dvector (0, n);
   Pion = Pi / n;
@@ -571,10 +571,10 @@ chebft_Zeros (double u[], int n, int inv)
 /* --------------------------------------------------------------------------*/
 
 void
-chebft_Extremes (double u[], int n, int inv)
+chebft_Extremes (CCTK_REAL u[], int n, int inv)
 {
   int k, j, isignum, N = n - 1;
-  double fac, sum, PioN, *c;
+  CCTK_REAL fac, sum, PioN, *c;
 
   c = dvector (0, N);
   PioN = Pi / N;
@@ -614,7 +614,7 @@ chebft_Extremes (double u[], int n, int inv)
 /* --------------------------------------------------------------------------*/
 
 void
-chder (double *c, double *cder, int n)
+chder (CCTK_REAL *c, CCTK_REAL *cder, int n)
 {
   int j;
 
@@ -625,10 +625,10 @@ chder (double *c, double *cder, int n)
 }
 
 /* --------------------------------------------------------------------------*/
-double
-chebev (double a, double b, double c[], int m, double x)
+CCTK_REAL
+chebev (CCTK_REAL a, CCTK_REAL b, CCTK_REAL c[], int m, CCTK_REAL x)
 {
-  double d = 0.0, dd = 0.0, sv, y, y2;
+  CCTK_REAL d = 0.0, dd = 0.0, sv, y, y2;
   int j;
 
   y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
@@ -643,10 +643,10 @@ chebev (double a, double b, double c[], int m, double x)
 
 /* --------------------------------------------------------------------------*/
 void
-fourft (double *u, int N, int inv)
+fourft (CCTK_REAL *u, int N, int inv)
 {
   int l, k, iy, M;
-  double x, x1, fac, Pi_fac, *a, *b;
+  CCTK_REAL x, x1, fac, Pi_fac, *a, *b;
 
   M = N / 2;
   a = dvector (0, M);
@@ -705,7 +705,7 @@ fourft (double *u, int N, int inv)
 
 /* -----------------------------------------*/
 void
-fourder (double u[], double du[], int N)
+fourder (CCTK_REAL u[], CCTK_REAL du[], int N)
 {
   int l, M, lpM;
 
@@ -722,7 +722,7 @@ fourder (double u[], double du[], int N)
 
 /* -----------------------------------------*/
 void
-fourder2 (double u[], double d2u[], int N)
+fourder2 (CCTK_REAL u[], CCTK_REAL d2u[], int N)
 {
   int l, l2, M, lpM;
 
@@ -739,11 +739,11 @@ fourder2 (double u[], double d2u[], int N)
 }
 
 /* ----------------------------------------- */
-double
-fourev (double *u, int N, double x)
+CCTK_REAL
+fourev (CCTK_REAL *u, int N, CCTK_REAL x)
 {
   int l, M = N / 2;
-  double xl, result;
+  CCTK_REAL xl, result;
 
   result = 0.5 * (u[0] + u[M] * cos (x * M));
   for (l = 1; l < M; l++)
@@ -756,12 +756,12 @@ fourev (double *u, int N, double x)
 
 /* --------------------------------------------------------------------------*/
 void
-ludcmp (double **a, int n, int *indx, double *d)
+ludcmp (CCTK_REAL **a, int n, int *indx, CCTK_REAL *d)
 {				/* Version of 'ludcmp' of the numerical recipes for*/
   /* matrices a[0:n-1][0:n-1]*/
   int i, imax, j, k;
-  double big, dum, sum, temp;
-  double *vv;
+  CCTK_REAL big, dum, sum, temp;
+  CCTK_REAL *vv;
 
   vv = dvector (0, n - 1);
   *d = 1.0;
@@ -823,12 +823,12 @@ ludcmp (double **a, int n, int *indx, double *d)
 
 /* --------------------------------------------------------------------------*/
 void
-lubksb (double **a, int n, int *indx, double b[])
+lubksb (CCTK_REAL **a, int n, int *indx, CCTK_REAL b[])
 {				/* Version of 'lubksb' of the numerical recipes for*/
   /* matrices a[0:n-1][0:n-1] and vectors b[0:n-1]*/
 
   int i, ii = 0, ip, j;
-  double sum;
+  CCTK_REAL sum;
 
   for (i = 0; i < n; i++)
   {
@@ -853,11 +853,11 @@ lubksb (double **a, int n, int *indx, double b[])
 
 /* -------------------------------------------------------------------------*/
 void
-tridag (double a[], double b[], double c[], double r[], double u[], int n)
+tridag (CCTK_REAL a[], CCTK_REAL b[], CCTK_REAL c[], CCTK_REAL r[], CCTK_REAL u[], int n)
 {				/* Version of 'tridag' of the numerical recipes for*/
   /* vectors a, b, c, r, u with indices in the range [0:n-1]*/
   int j;
-  double bet, *gam;
+  CCTK_REAL bet, *gam;
 
   gam = dvector (0, n - 1);
   if (b[0] == 0.0)
@@ -877,11 +877,11 @@ tridag (double a[], double b[], double c[], double r[], double u[], int n)
 }
 
 /* ------------------------------------------------------------------------*/
-double
-norm1 (double *v, int n)
+CCTK_REAL
+norm1 (CCTK_REAL *v, int n)
 {
   int i;
-  double result = -1;
+  CCTK_REAL result = -1;
 
   for (i = 0; i < n; i++)
     if (fabs (v[i]) > result)
@@ -891,11 +891,11 @@ norm1 (double *v, int n)
 }
 
 /* -------------------------------------------------------------------------*/
-double
-norm2 (double *v, int n)
+CCTK_REAL
+norm2 (CCTK_REAL *v, int n)
 {
   int i;
-  double result = 0;
+  CCTK_REAL result = 0;
 
   for (i = 0; i < n; i++)
     result += v[i] * v[i];
@@ -904,11 +904,11 @@ norm2 (double *v, int n)
 }
 
 /* -------------------------------------------------------------------------*/
-double
-scalarproduct (double *v, double *w, int n)
+CCTK_REAL
+scalarproduct (CCTK_REAL *v, CCTK_REAL *w, int n)
 {
   int i;
-  double result = 0;
+  CCTK_REAL result = 0;
 
   for (i = 0; i < n; i++)
     result += v[i] * w[i];
@@ -917,11 +917,11 @@ scalarproduct (double *v, double *w, int n)
 }
 
 /* -------------------------------------------------------------------------*/
-double
-plgndr (int l, int m, double x)
+CCTK_REAL
+plgndr (int l, int m, CCTK_REAL x)
 {
   void nrerror (char error_text[]);
-  double fact, pll, pmm, pmmp1, somx2;
+  CCTK_REAL fact, pll, pmm, pmmp1, somx2;
   int i, ll;
 
   if (m < 0 || m > l || fabs (x) > 1.0)
diff --git a/src/TP_utilities.h b/src/TP_utilities.h
index fb05cb2..95c2045 100644
--- a/src/TP_utilities.h
+++ b/src/TP_utilities.h
@@ -2,6 +2,8 @@
 
 #include <math.h>
 
+#include "cctk.h"
+
 #define Pi  3.14159265358979323846264338328
 #define Pih 1.57079632679489661923132169164	/* Pi/2*/
 #define Piq 0.78539816339744830961566084582	/* Pi/4*/
@@ -13,7 +15,7 @@
 
 typedef struct DCOMPLEX
 {
-  double r, i;
+  CCTK_REAL r, i;
 } dcomplex;
 
 #define nrerror TP_nrerror
@@ -30,16 +32,16 @@ typedef struct DCOMPLEX
 
 void nrerror (char error_text[]);
 int *ivector (long nl, long nh);
-double *dvector (long nl, long nh);
+CCTK_REAL *dvector (long nl, long nh);
 int **imatrix (long nrl, long nrh, long ncl, long nch);
-double **dmatrix (long nrl, long nrh, long ncl, long nch);
-double ***d3tensor (long nrl, long nrh, long ncl, long nch, long ndl,
+CCTK_REAL **dmatrix (long nrl, long nrh, long ncl, long nch);
+CCTK_REAL ***d3tensor (long nrl, long nrh, long ncl, long nch, long ndl,
 		    long ndh);
 void free_ivector (int *v, long nl, long nh);
-void free_dvector (double *v, long nl, long nh);
+void free_dvector (CCTK_REAL *v, long nl, long nh);
 void free_imatrix (int **m, long nrl, long nrh, long ncl, long nch);
-void free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch);
-void free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch,
+void free_dmatrix (CCTK_REAL **m, long nrl, long nrh, long ncl, long nch);
+void free_d3tensor (CCTK_REAL ***t, long nrl, long nrh, long ncl, long nch,
 		    long ndl, long ndh);
 
 int minimum2 (int i, int j);
@@ -48,19 +50,19 @@ int maximum2 (int i, int j);
 int maximum3 (int i, int j, int k);
 int pow_int (int mantisse, int exponent);
 #if 0
-double atanh (double x);
-double asinh (double x);
-double acosh (double x);
+CCTK_REAL atanh (CCTK_REAL x);
+CCTK_REAL asinh (CCTK_REAL x);
+CCTK_REAL acosh (CCTK_REAL x);
 #endif
 
 dcomplex Cadd (dcomplex a, dcomplex b);
 dcomplex Csub (dcomplex a, dcomplex b);
 dcomplex Cmul (dcomplex a, dcomplex b);
-dcomplex RCmul (double x, dcomplex a);
+dcomplex RCmul (CCTK_REAL x, dcomplex a);
 dcomplex Cdiv (dcomplex a, dcomplex b);
-dcomplex Complex (double re, double im);
+dcomplex Complex (CCTK_REAL re, CCTK_REAL im);
 dcomplex Conjg (dcomplex z);
-double Cabs (dcomplex z);
+CCTK_REAL Cabs (dcomplex z);
 
 dcomplex Csqrt (dcomplex z);
 dcomplex Cexp (dcomplex z);
@@ -74,22 +76,22 @@ dcomplex Ccosh (dcomplex z);
 dcomplex Ctanh (dcomplex z);
 dcomplex Ccoth (dcomplex z);
 
-void chebft_Zeros (double u[], int n, int inv);
-void chebft_Extremes (double u[], int n, int inv);
-void chder (double *c, double *cder, int n);
-double chebev (double a, double b, double c[], int m, double x);
-void fourft (double *u, int N, int inv);
-void fourder (double u[], double du[], int N);
-void fourder2 (double u[], double d2u[], int N);
-double fourev (double *u, int N, double x);
+void chebft_Zeros (CCTK_REAL u[], int n, int inv);
+void chebft_Extremes (CCTK_REAL u[], int n, int inv);
+void chder (CCTK_REAL *c, CCTK_REAL *cder, int n);
+CCTK_REAL chebev (CCTK_REAL a, CCTK_REAL b, CCTK_REAL c[], int m, CCTK_REAL x);
+void fourft (CCTK_REAL *u, int N, int inv);
+void fourder (CCTK_REAL u[], CCTK_REAL du[], int N);
+void fourder2 (CCTK_REAL u[], CCTK_REAL d2u[], int N);
+CCTK_REAL fourev (CCTK_REAL *u, int N, CCTK_REAL x);
 
 
-void ludcmp (double **a, int n, int *indx, double *d);
-void lubksb (double **a, int n, int *indx, double b[]);
-void tridag (double a[], double b[], double c[], double r[], double u[],
+void ludcmp (CCTK_REAL **a, int n, int *indx, CCTK_REAL *d);
+void lubksb (CCTK_REAL **a, int n, int *indx, CCTK_REAL b[]);
+void tridag (CCTK_REAL a[], CCTK_REAL b[], CCTK_REAL c[], CCTK_REAL r[], CCTK_REAL u[],
 	     int n);
-double norm1 (double *v, int n);
-double norm2 (double *v, int n);
-double scalarproduct (double *v, double *w, int n);
+CCTK_REAL norm1 (CCTK_REAL *v, int n);
+CCTK_REAL norm2 (CCTK_REAL *v, int n);
+CCTK_REAL scalarproduct (CCTK_REAL *v, CCTK_REAL *w, int n);
 
-double plgndr (int l, int m, double x);
+CCTK_REAL plgndr (int l, int m, CCTK_REAL x);
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index a48cf79..cfde44d 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -88,39 +88,39 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH,
         fprintf(debug_file,
                 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g "
                 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
-                s_x[indx], s_y[indx],
-                A,B,
-                U.d0[0],
-                (-cos(Pih * (2 * i + 1) / n1)-1.0),
-                U.d1[0],
-                U.d2[0],
-                U.d3[0],
-                U.d11[0],
-                U.d22[0],
-                U.d33[0],
-                v.d0[indx],
-                v.d1[indx],
-                v.d2[indx],
-                v.d3[indx],
-                v.d11[indx],
-                v.d22[indx],
-                v.d33[indx]
+                (double)s_x[indx], (double)s_y[indx],
+                (double)A,(double)B,
+                (double)U.d0[0],
+                (double)(-cos(Pih * (2 * i + 1) / n1)-1.0),
+                (double)U.d1[0],
+                (double)U.d2[0],
+                (double)U.d3[0],
+                (double)U.d11[0],
+                (double)U.d22[0],
+                (double)U.d33[0],
+                (double)v.d0[indx],
+                (double)v.d1[indx],
+                (double)v.d2[indx],
+                (double)v.d3[indx],
+                (double)v.d11[indx],
+                (double)v.d22[indx],
+                (double)v.d33[indx]
                 );
       }
     fprintf(debug_file, "\n\n");
     for (i=n2-10; i<n2; i++)
     {
-      double d;
+      CCTK_REAL d;
       indx = Index(0,0,i,0,1,n1,n2,n3);
       d = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,
               s_x[indx], 0.0, 0.0);
       fprintf(debug_file, "%.16g %.16g\n",
-                s_x[indx], d);
+                (double)s_x[indx], (double)d);
     }
     fprintf(debug_file, "\n\n");
     for (i=n2-10; i<n2-1; i++)
     {
-      double d;
+      CCTK_REAL d;
       int ip= Index(0,0,i+1,0,1,n1,n2,n3);
       indx = Index(0,0,i,0,1,n1,n2,n3);
       for (j=-10; j<10; j++)
@@ -129,7 +129,7 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH,
                 s_x[indx]+(s_x[ip]-s_x[indx])*j/10,
                 0.0, 0.0);
         fprintf(debug_file, "%.16g %.16g\n",
-                s_x[indx]+(s_x[ip]-s_x[indx])*j/10, d);
+                (double)(s_x[indx]+(s_x[ip]-s_x[indx])*j/10), (double)d);
       }
     }
     fprintf(debug_file, "\n\n");
@@ -152,13 +152,13 @@ void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH,
           C_To_c (nvar, X, R, &(s_x[indx]), &r, U);
           fprintf(debug_file,
                   "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
-                  s_x[indx], r, X, R, U.d0[0],
-                  U.d1[0],
-                  U.d2[0],
-                  U.d3[0],
-                  U.d11[0],
-                  U.d22[0],
-                  U.d33[0]);
+                  (double)s_x[indx], (double)r, (double)X, (double)R, (double)U.d0[0],
+                  (double)U.d1[0],
+                  (double)U.d2[0],
+                  (double)U.d3[0],
+                  (double)U.d11[0],
+                  (double)U.d22[0],
+                  (double)U.d33[0]);
         }
       }
     fclose(debug_file);
@@ -184,11 +184,12 @@ TwoPunctures (CCTK_ARGUMENTS)
 
   int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
+  int imin[3], imax[3], d;
   int i, j, k, ntotal = n1 * n2 * n3 * nvar;
   int percent10 = 0;
-  static double *F = NULL;
+  static CCTK_REAL *F = NULL;
   static derivs u, v;
-  double admMass;
+  CCTK_REAL admMass;
 
   if (! F) {
     /* Solve only when called for the first time */
@@ -232,7 +233,7 @@ TwoPunctures (CCTK_ARGUMENTS)
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (par_m_plus + par_m_minus
                - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", admMass);
+    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass);
   }
 
   if (CCTK_EQUALS(grid_setup_method, "Taylor expansion"))
@@ -252,8 +253,8 @@ TwoPunctures (CCTK_ARGUMENTS)
   averaged_lapse = CCTK_EQUALS(initial_lapse, "twopunctures-averaged");
 	pmn_lapse = CCTK_EQUALS(initial_lapse, "psi^n");
 	if (pmn_lapse)
-		CCTK_VInfo(CCTK_THORNSTRING, "Setting initial lapse to psi^%f profile.", 
-               initial_lapse_psi_exponent);
+		CCTK_VInfo(CCTK_THORNSTRING, "Setting initial lapse to psi^%f profile.",
+               (double)initial_lapse_psi_exponent);
 
   CCTK_INFO ("Interpolating result");
   if (CCTK_EQUALS(metric_type, "static conformal")) {
@@ -268,11 +269,17 @@ TwoPunctures (CCTK_ARGUMENTS)
     *conformal_state = 0;
   }
 
-  for (k = 0; k < cctk_lsh[2]; ++k)
+  for (d = 0; d < 3; ++ d)
   {
-    for (j = 0; j < cctk_lsh[1]; ++j)
+    imin[d] = 0           + (cctk_bbox[2*d  ] ? 0 : cctk_nghostzones[2]);
+    imax[d] = cctk_lsh[d] - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[2]);
+  }
+
+  for (k = imin[2]; k < imax[2]; ++k)
+  {
+    for (j = imin[1]; j < imax[1]; ++j)
     {
-      for (i = 0; i < cctk_lsh[0]; ++i)
+      for (i = imin[0]; i < imax[0]; ++i)
       {
         if (percent10 != 10*(i+j*cctk_lsh[0]+k*cctk_lsh[0]*cctk_lsh[1]) /
                             (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]))
@@ -284,12 +291,12 @@ TwoPunctures (CCTK_ARGUMENTS)
 
 	const int ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
         
-        double r_plus
+        CCTK_REAL r_plus
           = sqrt(pow(x[ind] - par_b, 2) + pow(y[ind], 2) + pow(z[ind], 2));
-        double r_minus
+        CCTK_REAL r_minus
           = sqrt(pow(x[ind] + par_b, 2) + pow(y[ind], 2) + pow(z[ind], 2));
         
-        double U;
+        CCTK_REAL U;
         switch (gsm)
         {
         case GSM_Taylor_expansion:
@@ -309,7 +316,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             r_plus = TP_Tiny;
         if (r_minus < TP_Tiny)
             r_minus = TP_Tiny;
-        double psi1 = 1
+        CCTK_REAL psi1 = 1
           + 0.5 * par_m_plus / r_plus
           + 0.5 * par_m_minus / r_minus + U;
 #define EXTEND(M,r) \
@@ -326,16 +333,16 @@ TwoPunctures (CCTK_ARGUMENTS)
              + 0.5 * EXTEND(par_m_minus,r_minus)
              + 0.5 * par_m_plus / r_plus + U;
         }
-        double static_psi = 1;
+        CCTK_REAL static_psi = 1;
         
-        double Aij[3][3];
+        CCTK_REAL Aij[3][3];
         BY_Aijofxyz (x[ind], y[ind], z[ind], Aij);
 
         if ((*conformal_state > 0) || (pmn_lapse)) {
 
-          double xp, yp, zp, rp, ir;
-          double s1, s3, s5;
-          double p, px, py, pz, pxx, pxy, pxz, pyy, pyz, pzz;
+          CCTK_REAL xp, yp, zp, rp, ir;
+          CCTK_REAL s1, s3, s5;
+          CCTK_REAL p, px, py, pz, pxx, pxy, pxz, pyy, pyz, pzz;
           p = 1.0;
           px = py = pz = 0.0;
           pxx = pxy = pxz = 0.0;
@@ -443,9 +450,9 @@ TwoPunctures (CCTK_ARGUMENTS)
         kzz[ind] = Aij[2][2] / pow(psi1, 2);
 
         if (antisymmetric_lapse || averaged_lapse) {
-/*           const double alp1 = ((1.0 - 0.5 * par_m_plus / r_plus) */
+/*           const CCTK_REAL alp1 = ((1.0 - 0.5 * par_m_plus / r_plus) */
 /*                                / (1.0 + 0.5 * par_m_plus / r_plus)); */
-/*           const double alp2 = ((1.0 - 0.5 * par_m_minus / r_minus) */
+/*           const CCTK_REAL alp2 = ((1.0 - 0.5 * par_m_minus / r_minus) */
 /*                                / (1.0 + 0.5 * par_m_minus / r_minus)); */
 /*           alp[ind] = alp1 * alp2; */
           alp[ind] =
diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h
index 14637a1..4c4399f 100644
--- a/src/TwoPunctures.h
+++ b/src/TwoPunctures.h
@@ -7,7 +7,7 @@
 
 typedef struct DERIVS
 {
-  double *d0, *d1, *d2, *d3, *d11, *d12, *d13, *d22, *d23, *d33;
+  CCTK_REAL *d0, *d1, *d2, *d3, *d11, *d12, *d13, *d22, *d23, *d33;
 } derivs;
 
 /*
@@ -22,8 +22,8 @@ Files of "TwoPunctures":
 */
 
 /* Routines in  "TwoPunctures.c"*/
-double TestSolution (double A, double B, double X, double R, double phi);
-void TestVector_w (double *par, int nvar, int n1, int n2, int n3, double *w);
+CCTK_REAL TestSolution (CCTK_REAL A, CCTK_REAL B, CCTK_REAL X, CCTK_REAL R, CCTK_REAL phi);
+void TestVector_w (CCTK_REAL *par, int nvar, int n1, int n2, int n3, CCTK_REAL *w);
 
 /* Routines in  "FuncAndJacobian.c"*/
 int Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3);
@@ -32,50 +32,50 @@ void free_derivs (derivs * v, int n);
 void Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v);
 void F_of_v (CCTK_POINTER_TO_CONST cctkGH,
        int nvar, int n1, int n2, int n3, derivs v,
-	     double *F, derivs u);
+	     CCTK_REAL *F, derivs u);
 void J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
-		 double *Jdv, derivs u);
+		 CCTK_REAL *Jdv, derivs u);
 void JFD_times_dv (int i, int j, int k, int nvar, int n1,
-		   int n2, int n3, derivs dv, derivs u, double *values);
+		   int n2, int n3, derivs dv, derivs u, CCTK_REAL *values);
 void SetMatrix_JFD (int nvar, int n1, int n2, int n3,
-		    derivs u, int *ncols, int **cols, double **Matrix);
-double PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
+		    derivs u, int *ncols, int **cols, CCTK_REAL **Matrix);
+CCTK_REAL PunctEvalAtArbitPosition (CCTK_REAL *v, CCTK_REAL A, CCTK_REAL B, CCTK_REAL phi,
 				 int n1, int n2, int n3);
 void calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1,
 		       int n2, int n3, derivs v, derivs vv);
-double interpol (double a, double b, double c, derivs v);
-double PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
-                                         int n2, int n3, derivs v, double x,
-                                         double y, double z);
-double PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
-				   int n2, int n3, derivs v, double x,
-				   double y, double z);
+CCTK_REAL interpol (CCTK_REAL a, CCTK_REAL b, CCTK_REAL c, derivs v);
+CCTK_REAL PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
+                                         int n2, int n3, derivs v, CCTK_REAL x,
+                                         CCTK_REAL y, CCTK_REAL z);
+CCTK_REAL PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+				   int n2, int n3, derivs v, CCTK_REAL x,
+				   CCTK_REAL y, CCTK_REAL z);
 
 /* Routines in  "CoordTransf.c"*/
-void AB_To_XR (int nvar, double A, double B, double *X,
-	       double *R, derivs U);
-void C_To_c (int nvar, double X, double R, double *x,
-	     double *r, derivs U);
-void rx3_To_xyz (int nvar, double x, double r, double phi, double *y,
-		 double *z, derivs U);
+void AB_To_XR (int nvar, CCTK_REAL A, CCTK_REAL B, CCTK_REAL *X,
+	       CCTK_REAL *R, derivs U);
+void C_To_c (int nvar, CCTK_REAL X, CCTK_REAL R, CCTK_REAL *x,
+	     CCTK_REAL *r, derivs U);
+void rx3_To_xyz (int nvar, CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi, CCTK_REAL *y,
+		 CCTK_REAL *z, derivs U);
 
 /* Routines in  "Equations.c"*/
-double BY_KKofxyz (double x, double y, double z);
-void BY_Aijofxyz (double x, double y, double z, double Aij[3][3]);
+CCTK_REAL BY_KKofxyz (CCTK_REAL x, CCTK_REAL y, CCTK_REAL z);
+void BY_Aijofxyz (CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL Aij[3][3]);
 void NonLinEquations (CCTK_REAL rho_adm,
-          double A, double B, double X, double R,
-		      double x, double r, double phi,
-		      double y, double z, derivs U, double *values);
-void LinEquations (double A, double B, double X, double R,
-		   double x, double r, double phi,
-		   double y, double z, derivs dU, derivs U, double *values);
+          CCTK_REAL A, CCTK_REAL B, CCTK_REAL X, CCTK_REAL R,
+		      CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi,
+		      CCTK_REAL y, CCTK_REAL z, derivs U, CCTK_REAL *values);
+void LinEquations (CCTK_REAL A, CCTK_REAL B, CCTK_REAL X, CCTK_REAL R,
+		   CCTK_REAL x, CCTK_REAL r, CCTK_REAL phi,
+		   CCTK_REAL y, CCTK_REAL z, derivs dU, derivs U, CCTK_REAL *values);
 
 /* Routines in  "Newton.c"*/
 void TestRelax (CCTK_POINTER_TO_CONST cctkGH,
-                int nvar, int n1, int n2, int n3, derivs v, double *dv);
+                int nvar, int n1, int n2, int n3, derivs v, CCTK_REAL *dv);
 void Newton (CCTK_POINTER_TO_CONST cctkGH,
              int nvar, int n1, int n2, int n3, derivs v,
-	           double tol, int itmax);
+	           CCTK_REAL tol, int itmax);
 
 
 /*