Tom's C version of IDAnalyticBH is now the default. Passes testsuites

for Schwarzschild and 2 misner holes, need to test multiple misner
and BL still


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAnalyticBH/trunk@40 6a3ddf76-46e1-4315-99d9-bc56cac1ef84

15 files changed, 1059 insertions, 725 deletions

diff --git a/interface.ccl b/interface.ccl
index 498195c..6465590 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -6,10 +6,10 @@ inherits: einstein
 
 private:
 
-real misner_workspace TYPE = GF
-{
-  csch,
-  coth,
-  r1,
-  r2
-} "Workspace arrays for Misner_standard"
+#real misner_workspace TYPE = GF
+#{
+#  csch,
+#  coth,
+#  r1,
+#  r2
+#} "Workspace arrays for Misner_standard"
diff --git a/param.ccl b/param.ccl
index 87fb919..a2c5f14 100644
--- a/param.ccl
+++ b/param.ccl
@@ -113,7 +113,6 @@ EXTENDS KEYWORD initial_data ""
 EXTENDS KEYWORD initial_lapse ""
 {
   "schwarz" :: "Set lapse to Schwarzschild"
-  "cadez"   :: "Set initial lapse to Cadez"
 }
 
 EXTENDS BOOLEAN use_conformal ""
diff --git a/schedule.ccl b/schedule.ccl
index 269f497..ad0c94f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -8,7 +8,7 @@ if (CCTK_Equals(initial_data,"schwarzschild") ||
 {
   schedule ParamChecker at CCTK_PARAMCHECK
   {
-    LANG: Fortran
+    LANG: C
   } "Construct parameters for analytic black hole solutions"
 }
 
@@ -17,7 +17,7 @@ if (CCTK_Equals(initial_data,"schwarzschild"))
 	
    schedule Schwarzschild at CCTK_INITIAL after InitialEinstein
    {
-     LANG: Fortran
+     LANG: C
    } "Construct initial data for a single Schwarzschild black hole"
 
 }
@@ -27,17 +27,17 @@ if (CCTK_Equals(initial_data,"bl_bh"))
    STORAGE: confac
    schedule BrillLindquist at CCTK_INITIAL
    {
-     LANG: Fortran
+     LANG: C
    } "Construct initial data for Brill Lindquist black holes"
 
 }
 
 if (CCTK_Equals(initial_data,"misner_bh")) 
 { 
-   STORAGE: confac, misner_workspace
+   STORAGE: confac
    schedule Misner_standard at CCTK_INITIAL
    {
-     LANG: Fortran
+     LANG: C
    } "Construct initial data for multiple Misner black holes"
 
 }
@@ -47,7 +47,7 @@ if (CCTK_Equals(initial_data,"multiple_misner_bh"))
    STORAGE: confac
    schedule Misner_multiple at CCTK_INITIAL
    {
-     LANG: Fortran
+     LANG: C
    } "Construct initial data for two Misner black holes"
 
 }
diff --git a/src/BrillLindquist.F b/src/BrillLindquist.F
deleted file mode 100644
index d052668..0000000
--- a/src/BrillLindquist.F
+++ /dev/null
@@ -1,191 +0,0 @@
- /*@@
-   @file      BrillLindquist.F
-   @date      
-   @author    
-   @desc 
-      Set up initial data for Brill Lindquist black holes
-   @enddesc 
- @@*/
-
-#include "cctk.h"
-#include "cctk_arguments.h"
-#include "cctk_parameters.h"
-
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
- /*@@
-   @routine    BrillLindquist
-   @date       
-   @author     
-   @desc 
-      Set up initial data for Brill Lindquist black holes
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-      subroutine BrillLindquist(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer i
-      CCTK_REAL, dimension(4) :: hole_mass, hole_x0, hole_y0, hole_z0
-
-
-c     Put parameters into arrays for following calculations
-c     -----------------------------------------------------
-      hole_x0(1)   = bl_x0_1
-      hole_y0(1)   = bl_y0_1
-      hole_z0(1)   = bl_z0_1
-      hole_mass(1) = bl_M_1
-
-      hole_x0(2)   = bl_x0_2
-      hole_y0(2)   = bl_y0_2
-      hole_z0(2)   = bl_z0_2
-      hole_mass(2) = bl_M_2
-
-      hole_x0(3)   = bl_x0_3
-      hole_y0(3)   = bl_y0_3
-      hole_z0(3)   = bl_z0_3
-      hole_mass(3) = bl_M_3
-
-      hole_x0(4)   = bl_x0_4
-      hole_y0(4)   = bl_y0_4
-      hole_z0(4)   = bl_z0_4
-      hole_mass(4) = bl_M_4
-
-c     Do not ask...
-c     ----------
-      hole_x0 = -hole_x0
-      hole_y0 = -hole_y0
-      hole_z0 = -hole_z0
-
-c     Initialize to zero and then use +=
-c     ----------------------------------
-      psi = 1.0D0
-
-      if (use_conformal_derivs == 1) then
-
-        psix = 0.0D0
-        psiy = 0.0D0
-        psiz = 0.0D0
-        psixx = 0.0D0
-        psixy = 0.0D0
-        psixz = 0.0D0
-        psiyy = 0.0D0
-        psiyz = 0.0D0
-        psizz = 0.0D0
-
-      end if
-
-      do i = 1, bl_nbh
-
-c     Maple Output
-c     ------------
-      Psi = Psi+hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2.D
-     &0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)+h
-     &ole_z0(i)**2.D0)/2.D0
-
-      if (use_conformal_derivs == 1) then
-
-      Psix = Psix-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2
-     &.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)
-     &+hole_z0(i)**2.D0)**3.D0*(2.D0*x+2.D0*hole_x0(i))/4.D0
-      Psiy = Psiy-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2
-     &.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)
-     &+hole_z0(i)**2.D0)**3.D0*(2.D0*y+2.D0*hole_y0(i))/4.D0
-      Psiz = Psiz-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2
-     &.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)
-     &+hole_z0(i)**2.D0)**3.D0*(2.D0*z+2.D0*hole_z0(i))/4.D0
-      Psixx = Psixx+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*x+2.D0*hole_x0(i))**2.D
-     &0-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2.D0+y**2+2
-     &.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)+hole_z0(i
-     &)**2.D0)**3.D0/2.D0
-      Psixy = Psixy+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*x+2.D0*hole_x0(i))*(2.D
-     &0*y+2.D0*hole_y0(i))
-      Psixz = Psixz+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*x+2.D0*hole_x0(i))*(2.D
-     &0*z+2.D0*hole_z0(i))
-      Psiyy = Psiyy+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*y+2.D0*hole_y0(i))**2.D
-     &0-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2.D0+y**2+2
-     &.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)+hole_z0(i
-     &)**2.D0)**3.D0/2.D0
-      Psiyz = Psiyz+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*y+2.D0*hole_y0(i))*(2.D
-     &0*z+2.D0*hole_z0(i))
-      Psizz = Psizz+3.D0/8.D0*hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+h
-     &ole_x0(i)**2.D0+y**2+2.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*
-     &z*hole_z0(i)+hole_z0(i)**2.D0)**5.D0*(2.D0*z+2.D0*hole_z0(i))**2.D
-     &0-hole_mass(i)/sqrt(x**2+2.D0*x*hole_x0(i)+hole_x0(i)**2.D0+y**2+2
-     &.D0*y*hole_y0(i)+hole_y0(i)**2.D0+z**2+2.D0*z*hole_z0(i)+hole_z0(i
-     &)**2.D0)**3.D0/2.D0
-
-      end if
-
-      end do
-
-c     Cactus conventions
-c     ------------------
-      if (use_conformal_derivs == 1) then
-
-        psix = psix / psi
-        psiy = psiy / psi
-        psiz = psiz / psi
-        psixx = psixx / psi
-        psixy = psixy / psi
-        psixz = psixz / psi
-        psiyy = psiyy / psi
-        psiyz = psiyz / psi
-        psizz = psizz / psi
-
-      end if
-
-c     Metric depends on conformal choice
-c     ----------------------------------
-      if (conformal_state == CONFORMAL_METRIC) then
-
-        gxx = 1d0
-        gyy = 1d0
-        gzz = 1d0
-        gxy = 0d0 
-        gxz = 0d0  
-        gyz = 0d0   
-
-      else 
-
-        gxx = psi**4
-        gyy = gxx
-        gzz = gxx
-        gxy = 0d0
-        gyz = 0d0
-        gzz = 0d0
-
-      end if
-
-c     Time-symmetric data
-c     -------------------
-      kxx = 0d0
-      kyy = 0d0
-      kzz = 0d0
-      kxy = 0d0 
-      kxz = 0d0  
-      kyz = 0d0   
-
-      return
-      end
-
diff --git a/src/BrillLindquist.c b/src/BrillLindquist.c
new file mode 100644
index 0000000..bc890c0
--- /dev/null
+++ b/src/BrillLindquist.c
@@ -0,0 +1,249 @@
+ /*@@
+   @file      BrillLindquist.F
+   @date      
+   @author    
+   @desc 
+      Set up initial data for Brill Lindquist black holes
+   @enddesc 
+ @@*/
+
+#include <math.h>
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+#define SQR(a) ((a)*(a))
+
+#define MAX_HOLES 4
+ /*@@
+   @routine    BrillLindquist
+   @date       
+   @author     
+   @desc 
+      Set up initial data for Brill Lindquist black holes
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void BrillLindquist(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  int n;
+  CCTK_REAL hole_mass[MAX_HOLES], hole_x0[MAX_HOLES], hole_y0[MAX_HOLES], hole_z0[MAX_HOLES];
+  CCTK_REAL tmp1, tmp2, tmp3;
+  CCTK_REAL xval, yval, zval;
+  CCTK_REAL x_2, y_2, z_2;
+  int i,j,k;
+  int nx,ny,nz;
+  int index;
+
+  nx = cctk_lsh[0];
+  ny = cctk_lsh[1];
+  nz = cctk_lsh[2];
+
+  /*     Put parameters into arrays for following calculations
+   *     -----------------------------------------------------
+   */
+  hole_x0[0]   = -bl_x0_1;
+  hole_y0[0]   = -bl_y0_1;
+  hole_z0[0]   = -bl_z0_1;
+  hole_mass[0] = bl_M_1;
+
+  hole_x0[1]   = -bl_x0_2;
+  hole_y0[1]   = -bl_y0_2;
+  hole_z0[1]   = -bl_z0_2;
+  hole_mass[1] = bl_M_2;
+
+  hole_x0[2]   = -bl_x0_3;
+  hole_y0[2]   = -bl_y0_3;
+  hole_z0[2]   = -bl_z0_3;
+  hole_mass[2] = bl_M_3;
+
+  hole_x0[3]   = -bl_x0_4;
+  hole_y0[3]   = -bl_y0_4;
+  hole_z0[3]   = -bl_z0_4;
+  hole_mass[3] = bl_M_4;
+
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+
+        /*     Initialize to zero and then use +=
+         *     ----------------------------------
+         */
+
+        psi[index] = 1.0;
+
+        if (use_conformal_derivs == 1)
+        {
+          psix[index]  = 0.0;
+          psiy[index]  = 0.0;
+          psiz[index]  = 0.0;
+          psixx[index] = 0.0;
+          psixy[index] = 0.0;
+          psixz[index] = 0.0;
+          psiyy[index] = 0.0;
+          psiyz[index] = 0.0;
+          psizz[index] = 0.0;
+        }
+
+        x_2 = SQR(x[index]);
+        y_2 = SQR(y[index]);
+
+        xval = x[index];
+        yval = y[index];
+        zval = z[index];
+
+        for(n = 0;  n < bl_nbh; i++)
+        {
+
+          /*     Maple Output
+           *     ------------
+           */
+          
+          tmp1 = sqrt(x_2+2.0*xval*hole_x0[n]+SQR(hole_x0[n])
+                     +y_2+2.0*yval*hole_y0[n]+SQR(hole_y0[n])
+                     +z_2+2.0*zval*hole_z0[n]+SQR(hole_z0[n]));
+
+          psi[index] += hole_mass[n]/tmp1/2.0;
+
+          if (use_conformal_derivs == 1)
+          {
+            tmp2 = pow(tmp1, 3);
+
+            psix[index] +=  -hole_mass[n]/tmp2*(2.0*xval+2.0*hole_x0[n])/4.0;
+            psiy[index] +=  -hole_mass[n]/tmp2*(2.0*yval+2.0*hole_y0[n])/4.0;
+            psiz[index] +=  -hole_mass[n]/tmp2*(2.0*zval+2.0*hole_z0[n])/4.0;
+
+            tmp2 = pow(tmp1, 5.0);
+            tmp3 = pow(tmp1, 3.0);
+
+            psixx[index] += 3.0/8.0*hole_mass[n]/tmp2*SQR(2.0*xval+2.0*hole_x0[n])-hole_mass[n]/tmp3/2.0;
+        
+            psixy[index] += 3.0/8.0*hole_mass[n]/tmp2*(2.0*xval+2.0*hole_x0[n])*(2.0*yval+2.0*hole_y0[n]);
+            psixz[index] += 3.0/8.0*hole_mass[n]/tmp2*(2.0*xval+2.0*hole_x0[n])*(2.0*zval+2.0*hole_z0[n]);
+
+            psiyy[index] += 3.0/8.0*hole_mass[n]/tmp2*SQR(2.0*yval+2.0*hole_y0[n])-hole_mass[n]/tmp3/2.0;
+            psiyz[index] += 3.0/8.0*hole_mass[n]/tmp2*(2.0*yval+2.0*hole_y0[n])*(2.0*zval+2.0*hole_z0[n]);
+
+            psizz[index] += 3.0/8.0*hole_mass[n]/tmp2*SQR(2.0*zval+2.0*hole_z0[n])-hole_mass[n]/tmp3/2.0;
+          }
+        }
+      }
+    }
+  }
+
+  /*     Cactus conventions
+   *     ------------------
+   */
+
+  if (use_conformal_derivs == 1)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          psix[index]  /= psi[index];
+          psiy[index]  /= psi[index];
+          psiz[index]  /= psi[index];
+          psixx[index] /= psi[index];
+          psixy[index] /= psi[index];
+          psixz[index] /= psi[index];
+          psiyy[index] /= psi[index];
+          psiyz[index] /= psi[index];
+          psizz[index] /= psi[index];
+        }
+      }
+    }
+  }
+
+  /*     Metric depends on conformal state
+   *     ---------------------------------
+   */
+
+  if (*conformal_state == CONFORMAL_METRIC)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+
+          gxx[index] = 1.0;
+          gyy[index] = 1.0;
+          gzz[index] = 1.0;
+          gxy[index] = 0.0;
+          gxz[index] = 0.0;
+          gyz[index] = 0.0;
+        }
+      }
+    }
+  }
+  else 
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = pow(psi[index],4);
+          gyy[index] = gxx[index];
+          gzz[index] = gxx[index];
+          gxy[index] = 0.0;
+          gxz[index] = 0.0;
+          gyz[index] = 0.0;
+        }
+      }
+    }
+  }
+
+  /*     Time-symmetric data
+   *     -------------------
+   */
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+        kxx[index] = 0.0;
+        kyy[index] = 0.0;
+        kzz[index] = 0.0;
+        kxy[index] = 0.0;
+        kxz[index] = 0.0;
+        kyz[index] = 0.0;
+      }
+    }
+  }
+
+  return;
+}
diff --git a/src/Misner_multiple.F b/src/Misner_multiple.F
deleted file mode 100644
index d5b729e..0000000
--- a/src/Misner_multiple.F
+++ /dev/null
@@ -1,153 +0,0 @@
- /*@@
-   @file      Misner_multiple.F
-   @date      
-   @author    Carsten Gundlach
-   @desc 
-      Set up initial data for multiple Misner black holes
-   @enddesc 
- @@*/
-
-#include "cctk.h"
-#include "cctk_arguments.h"
-#include "cctk_parameters.h"
-
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
- /*@@
-   @routine    Misner_multiple
-   @date       
-   @author     Carsten Gundlach
-   @desc 
-      Set up initial data for multiple Misner black holes
-   @enddesc 
-   @calls     
-   @history 
- 
-   @endhistory 
-
-@@*/
-      subroutine Misner_multiple(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer :: i, j, k
-      CCTK_REAL :: xval, yval, zval
-      CCTK_REAL :: one, zero, tmp0, tmp1, tmp2, tmp3, tmp4
-      CCTK_REAL, parameter :: nm_eps = 1.d-6   ! finite differencing step
-
-      one  = 1.0D0
-      zero = 0.0D0
-
-c     Initialize C global variables
-c     -----------------------------
-      call Misner_init(misner_nbh,mu,nmax)
-
-C     Get value of psi pointwise from a C function
-c     --------------------------------------------
-      do k=1,cctk_lsh(1)
-         do j=1,cctk_lsh(2)
-            do i=1,cctk_lsh(3)
- 
-               xval = x(i,j,k)
-               yval = y(i,j,k)
-               zval = z(i,j,k)
-
-               call MisnerEvalPsi(xval,yval,zval,tmp0)
-               psi(i,j,k) = tmp0
-               
-c              Only calculate derivatives of psi if required
-c              ---------------------------------------------
-               if (use_conformal_derivs == 1) then
-
-                 call MisnerEvalPsi(xval+nm_eps,yval,zval,tmp1)
-                 call MisnerEvalPsi(xval-nm_eps,yval,zval,tmp2)
-                 psix(i,j,k) = 0.5d0*(tmp1-tmp2)/nm_eps
-                 psixx(i,j,k) = (tmp1+tmp2-2.d0*tmp0)/nm_eps**2
-
-                 call MisnerEvalPsi(xval,yval+nm_eps,zval,tmp1)
-                 call MisnerEvalPsi(xval,yval-nm_eps,zval,tmp2)
-                 psiy(i,j,k) = 0.5d0*(tmp1-tmp2)/nm_eps
-                 psiyy(i,j,k) = (tmp1+tmp2-2.d0*tmp0)/nm_eps**2
-
-                 call MisnerEvalPsi(xval,yval,zval+nm_eps,tmp1)
-                 call MisnerEvalPsi(xval,yval,zval-nm_eps,tmp2)
-                 psiz(i,j,k) = 0.5d0*(tmp1-tmp2)/nm_eps
-                 psizz(i,j,k) = (tmp1+tmp2-2.d0*tmp0)/nm_eps**2
-
-                 call MisnerEvalPsi(xval+nm_eps,yval+nm_eps,zval,tmp1)
-                 call MisnerEvalPsi(xval+nm_eps,yval-nm_eps,zval,tmp2)
-                 call MisnerEvalPsi(xval-nm_eps,yval+nm_eps,zval,tmp3)
-                 call MisnerEvalPsi(xval-nm_eps,yval-nm_eps,zval,tmp4)
-                 psixy(i,j,k) = 0.25d0*(tmp1-tmp2-tmp3+tmp4)/nm_eps**2
-
-                 call MisnerEvalPsi(xval,yval+nm_eps,zval+nm_eps,tmp1)
-                 call MisnerEvalPsi(xval,yval-nm_eps,zval+nm_eps,tmp2)
-                 call MisnerEvalPsi(xval,yval+nm_eps,zval-nm_eps,tmp3)
-                 call MisnerEvalPsi(xval,yval-nm_eps,zval-nm_eps,tmp4)
-                 psiyz(i,j,k) = 0.25d0*(tmp1-tmp2-tmp3+tmp4)/nm_eps**2
-
-                 call MisnerEvalPsi(xval+nm_eps,yval,zval+nm_eps,tmp1)
-                 call MisnerEvalPsi(xval+nm_eps,yval,zval-nm_eps,tmp2)
-                 call MisnerEvalPsi(xval-nm_eps,yval,zval+nm_eps,tmp3)
-                 call MisnerEvalPsi(xval-nm_eps,yval,zval-nm_eps,tmp4)
-                 psixz(i,j,k) = 0.25d0*(tmp1-tmp2-tmp3+tmp4)/nm_eps**2
-
-              end if
-
-            end do
-         end do
-      end do
-
-c     Cactus conventions
-c     ------------------
-      if (use_conformal_derivs == 1) then
-
-        psix = psix/psi
-        psiy = psiy/psi
-        psiz = psiz/psi
-        psixx = psixx/psi
-        psixy = psixy/psi
-        psixz = psixz/psi
-        psiyy = psiyy/psi
-        psiyz = psiyz/psi
-        psizz = psizz/psi
-
-      end if
-
-c     Metric depends on conformal state
-c     ---------------------------------
-      if (conformal_state == CONFORMAL_METRIC) then
-
-        gxx = one
-        gyy = one
-        gzz = one
-        gxy = zero 
-        gxz = zero  
-        gyz = zero   
-
-      else 
-
-        gxx = psi**4
-        gyy = gxx
-        gzz = gxx
-        gxy = zero
-        gxz = zero
-        gyz = zero
-
-      end if
-
-c     Time-symmetric data
-c     -------------------
-      kxx = zero
-      kyy = zero
-      kzz = zero
-      kxy = zero 
-      kxz = zero  
-      kyz = zero   
-
-      return
-      end
-
diff --git a/src/Misner_multiple.c b/src/Misner_multiple.c
new file mode 100644
index 0000000..2ba3a1d
--- /dev/null
+++ b/src/Misner_multiple.c
@@ -0,0 +1,219 @@
+ /*@@
+   @file      Misner_multiple.F
+   @date      
+   @author    Carsten Gundlach
+   @desc 
+      Set up initial data for multiple Misner black holes
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+#define SQR(a) ((a)*(a))
+
+void Misner_init(int n, CCTK_REAL mu, int terms);
+void MisnerEvalPsi(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL *res); 
+
+ /*@@
+   @routine    Misner_multiple
+   @date       
+   @author     Carsten Gundlach
+   @desc 
+      Set up initial data for multiple Misner black holes
+   @enddesc 
+   @calls     
+   @history 
+ 
+   @endhistory 
+
+@@*/
+void Misner_multiple(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  int i, j, k;
+  CCTK_REAL xval, yval, zval;
+  CCTK_REAL one, zero, tmp0, tmp1, tmp2, tmp3, tmp4;
+  CCTK_REAL nm_eps = 0.000001;   /* finite differencing step*/
+
+  int nx,ny,nz;
+  int index;
+
+  nx = cctk_lsh[0];
+  ny = cctk_lsh[1];
+  nz = cctk_lsh[2];
+
+  one  = 1.0;
+  zero = 0.0;
+
+  /*     Initialize C global variables
+   *     -----------------------------
+   */
+
+  Misner_init(misner_nbh,mu,nmax);
+
+  /*     Get value of psi pointwise from a C function
+   *     --------------------------------------------
+   */
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+
+        xval = x[index];
+        yval = y[index];
+        zval = z[index];
+
+        MisnerEvalPsi(xval,yval,zval,&tmp0);
+        psi[index] = tmp0;
+               
+        /*              Only calculate derivatives of psi if required
+         *              ---------------------------------------------
+         */
+        if (use_conformal_derivs == 1)
+        {
+	  printf("Hello\n");
+          MisnerEvalPsi(xval+nm_eps,yval,zval,&tmp1);
+          MisnerEvalPsi(xval-nm_eps,yval,zval,&tmp2);
+          psix[index] = 0.5*(tmp1-tmp2)/nm_eps;
+          psixx[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval+nm_eps,zval,&tmp1);
+          MisnerEvalPsi(xval,yval-nm_eps,zval,&tmp2);
+          psiy[index]  = 0.5*(tmp1-tmp2)/nm_eps;
+          psiyy[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval,yval,zval-nm_eps,&tmp2);
+          psiz[index] = 0.5*(tmp1-tmp2)/nm_eps;
+          psizz[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval+nm_eps,yval+nm_eps,zval,&tmp1);
+          MisnerEvalPsi(xval+nm_eps,yval-nm_eps,zval,&tmp2);
+          MisnerEvalPsi(xval-nm_eps,yval+nm_eps,zval,&tmp3);
+          MisnerEvalPsi(xval-nm_eps,yval-nm_eps,zval,&tmp4);
+          psixy[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval+nm_eps,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval,yval-nm_eps,zval+nm_eps,&tmp2);
+          MisnerEvalPsi(xval,yval+nm_eps,zval-nm_eps,&tmp3);
+          MisnerEvalPsi(xval,yval-nm_eps,zval-nm_eps,&tmp4);
+          psiyz[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval+nm_eps,yval,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval+nm_eps,yval,zval-nm_eps,&tmp2);
+          MisnerEvalPsi(xval-nm_eps,yval,zval+nm_eps,&tmp3);
+          MisnerEvalPsi(xval-nm_eps,yval,zval-nm_eps,&tmp4);
+          psixz[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+        }
+
+      }
+    }
+  }
+
+  /*     Cactus conventions
+   *     ------------------
+   */
+
+  if (use_conformal_derivs == 1)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          psix[index]  /= psi[index];
+          psiy[index]  /= psi[index];
+          psiz[index]  /= psi[index];
+          psixx[index] /= psi[index];
+          psixy[index] /= psi[index];
+          psixz[index] /= psi[index];
+          psiyy[index] /= psi[index];
+          psiyz[index] /= psi[index];
+          psizz[index] /= psi[index];
+        }
+      }
+    }
+  }
+
+  /*     Metric depends on conformal state
+   *     ---------------------------------
+   */
+
+  if (*conformal_state == CONFORMAL_METRIC)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = one;
+          gyy[index] = one;
+          gzz[index] = one;
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+  else 
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = pow(psi[index],4);
+          gyy[index] = gxx[index];
+          gzz[index] = gxx[index];
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+
+  /*     Time-symmetric data
+   *     -------------------
+   */
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+        kxx[index] = zero;
+        kyy[index] = zero;
+        kzz[index] = zero;
+        kxy[index] = zero;
+        kxz[index] = zero;
+        kyz[index] = zero;
+      }
+    }
+  }
+
+  return;
+}
diff --git a/src/Misner_points.c b/src/Misner_points.c
index 2eb5541..5df8fb7 100644
--- a/src/Misner_points.c
+++ b/src/Misner_points.c
@@ -38,10 +38,10 @@ struct bhole {
 /* The seed black holes. */
 struct bhole bholes[MAXBHOLES];
 
-CCTK_REAL csch(CCTK_REAL mu) {
+static CCTK_REAL csch(CCTK_REAL mu) {
   return 1.0/sinh(mu); 
 }
-CCTK_REAL coth(CCTK_REAL mu) {
+static CCTK_REAL coth(CCTK_REAL mu) {
   return cosh(mu)/sinh(mu); 
 }
 
@@ -59,7 +59,7 @@ CCTK_REAL coth(CCTK_REAL mu) {
 
 @@*/
 
-void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3) 
+static void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3) 
 {
   CCTK_REAL rad,radtwo;
   radtwo=(
@@ -87,7 +87,7 @@ void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3)
 
 @@*/
 
-void fill_iso(struct bhole *b, int n) 
+static void fill_iso(struct bhole *b, int n) 
 {
   int i,j;
   if(n==0) 
@@ -122,29 +122,29 @@ void fill_iso(struct bhole *b, int n)
 
 @@*/
 
-void FMODIFIER FORTRAN_NAME(Misner_init)(int *n, CCTK_REAL *mu, int *terms) 
+void Misner_init(int n, CCTK_REAL mu, int terms) 
 {
 
   int i;
   CCTK_REAL pi,ang;
 
-  assert((nbholes=*n) < MAXBHOLES);
+  assert((nbholes=n) < MAXBHOLES);
 
   pi = 4.0*atan(1.);
 
-  ang = 2.*pi/(*n);
+  ang = 2.*pi/(n);
 
-  for(i=0;i<*n;i++) 
+  for(i=0;i<n;i++) 
   {
-    bholes[i].x = coth(*mu)*cos(ang*i);
-    bholes[i].y = coth(*mu)*sin(ang*i);
-    bholes[i].mass = csch(*mu);
+    bholes[i].x = coth(mu)*cos(ang*i);
+    bholes[i].y = coth(mu)*sin(ang*i);
+    bholes[i].mass = csch(mu);
     bholes[i].i = i;
     bholes[i].isos = 0;
   }
 
-  for(i=0;i<*n;i++)
-    fill_iso(&bholes[i],*terms);
+  for(i=0;i<n;i++)
+    fill_iso(&bholes[i],terms);
 
 }
 
@@ -163,7 +163,7 @@ void FMODIFIER FORTRAN_NAME(Misner_init)(int *n, CCTK_REAL *mu, int *terms)
 
 @@*/
 
-CCTK_REAL eval_bh_psi(struct bhole *b, CCTK_REAL x, CCTK_REAL y, CCTK_REAL z) 
+static CCTK_REAL eval_bh_psi(struct bhole *b, CCTK_REAL x, CCTK_REAL y, CCTK_REAL z) 
 {
   int i;
   CCTK_REAL res;
@@ -198,10 +198,10 @@ CCTK_REAL eval_bh_psi(struct bhole *b, CCTK_REAL x, CCTK_REAL y, CCTK_REAL z)
 
 @@*/
 
-void FMODIFIER FORTRAN_NAME(MisnerEvalPsi)(CCTK_REAL *x, CCTK_REAL *y, CCTK_REAL *z, CCTK_REAL *res) 
+void MisnerEvalPsi(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL *res) 
 {
   int i;
   *res = 1;
   for(i=0;i<nbholes;i++)
-    *res += eval_bh_psi(&bholes[i],*x,*y,*z);
+    *res += eval_bh_psi(&bholes[i],x,y,z);
 }
diff --git a/src/Misner_standard.F b/src/Misner_standard.F
deleted file mode 100644
index d43f30f..0000000
--- a/src/Misner_standard.F
+++ /dev/null
@@ -1,182 +0,0 @@
- /*@@
-   @file      Misner_standard.F
-   @date      March 1997
-   @author    Joan Masso
-   @desc 
-      Set up initial data for two Misner black holes
-   @enddesc 
- @@*/
-
-#include "cctk.h"
-#include "cctk_arguments.h"
-#include "cctk_parameters.h"
-
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
- /*@@
-   @routine    Misner_standard
-   @date       
-   @author     Joan Masso, Ed Seidel
-   @desc 
-            Initialize the metric with a time symmetrical 
-            black hole spacetime containing
-            two axially symmetric misner black holes with a 
-            mass/length parameter mu. The mass is computed.
-            The spacetime line element has the form:
-                 $$ ds^2 = -dt^2 + \Psi^4 (dx^2+dy^2+dz^2) $$
-            and only $\Psi$ differs.
-            (Conformal factor from Karen Camarda)
-   @enddesc 
-   @calls     
-   @history 
- 
-   @endhistory 
-
-   @par      mu
-   @pdesc    Misner parameter. 
-   @ptype    real
-   @pcomment Values less than 1.8 do not really correspond to two 
-            black holes, as there is an initial single event horizon 
-            surrounding the throats. So, with low values of mu we also
-            have distorted single black holes.
-   @endpar
- 
-   @par      nmax
-   @pdesc    Summation limit for the misner series in the 'twobh' case.
-   @ptype    integer
-   @endpar
-
-@@*/
-      subroutine Misner_standard(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer :: n
-      INTEGER :: CCTK_Equals
-      CCTK_REAL :: zero, one, three
-
-      CHARACTER*200 :: infoline
-
-      zero  = 0.d0
-      one   = 1.d0
-      three = 3.d0
-
-c     Initialize so we can accumulate
-c     -------------------------------
-      psi = one
-      if (use_conformal_derivs == 1) then
-         psix = zero
-         psiy = zero
-         psiz = zero
-         psixx = zero
-         psixy = zero
-         psixz = zero
-         psiyy = zero
-         psiyz = zero
-         psizz = zero
-      end if
-
-      do n = nmax,1,-1
-         csch = one/sinh(mu*n)
-         coth = one/tanh(mu*n)
-         r1 = sqrt(x**2+y**2+(z+coth)**2)
-         r2 = sqrt(x**2+y**2+(z-coth)**2)
-         psi = psi + csch*(one/r1 + one/r2)
-
-         if (use_conformal_derivs == 1) then
-            psix = psix + (-(x/r2**3)-x/r1**3)*csch
-            psiy = psiy + (-(y/r2**3)-y/r1**3)*csch
-            psiz = psiz + (-((z-coth)/r2**3)-(z + coth)/r1**3)*csch
-            psixx = psixx +
-     &           ((three*x**2)/r2**5-one/r2**3+(three*x**2)/r1**5-one/r1**3)*csch
-            psixy = psixy + ((three*x*y)/r2**5+(three*x*y)/r1**5)*csch
-            psixz = psixz + ((three*x*(z-coth))/r2**5
-     &           +(three*x*(z+coth))/r1**5)*csch
-            psiyy = psiyy +
-     &           ((three*y**2)/r2**5-one/r2**3
-     &           +(three*y**2)/r1**5-one/r1**3)*csch
-            psiyz = psiyz + ((three*y*(z-coth))/r2**5
-     &           +(three*y*(z+coth))/r1**5)*csch
-            psizz = psizz + (-one/r2**3+(three*(z-coth)**2)/r2**5+
-     &           (three*(z+coth)**2)/r1**5-one/r1**3)*csch
-
-         end if
-
-      end do
-
-c     Cactus convention
-c     -----------------
-      if (use_conformal_derivs == 1) then
-         psix = psix/psi
-         psiy = psiy/psi
-         psiz = psiz/psi
-         psixx = psixx/psi
-         psixy = psixy/psi
-         psixz = psixz/psi
-         psiyy = psiyy/psi
-         psiyz = psiyz/psi
-         psizz = psizz/psi
-      end if
-
-c     compute the ADM mass
-c     -------------------- 
-      mass = zero
-      do n = 1,nmax
-         mass = mass + 4./sinh(n*mu)
-      enddo
-      call CCTK_INFO("Two misner black holes")
-      write(infoline,'(A17,G12.7)') 
-     &     'ADM mass of each ',mass/2d0
-      call CCTK_INFO(infoline)
-
-c     Should initialize lapse to Cadez value if possible
-c     --------------------------------------------------
-      if (CCTK_Equals(initial_lapse,"cadez") == 1) then 
-         write (*,*)"Initialise with cadez lapse"
-         alp = one
-         do n = 1,nmax
-            coth = one/tanh(mu*n)
-            r1 = sqrt(x**2+y**2+(z+coth)**2)
-            r2 = sqrt(x**2+y**2+(z-coth)**2)
-            alp = alp + (-1.)**n * one/sinh(mu*n)*(one/r1 + one/r2)
-         enddo
-         alp = alp/psi
-      endif
-         
-      
-c     Metric depends on conformal state
-c     ---------------------------------
-      if (conformal_state == CONFORMAL_METRIC) then
-
-        gxx = one
-        gyy = one
-        gzz = one
-        gxy = zero 
-        gxz = zero  
-        gyz = zero   
-
-      else 
-
-        gxx = psi**4
-        gyy = gxx
-        gzz = gxx
-        gxy = zero
-        gxz = zero
-        gyz = zero
-
-      end if
-
-c     Time-symmetric data
-c     -------------------
-      kxx = zero
-      kyy = zero
-      kzz = zero
-      kxy = zero 
-      kxz = zero  
-      kyz = zero   
-
-      return
-      end
diff --git a/src/Misner_standard.c b/src/Misner_standard.c
new file mode 100644
index 0000000..628eee6
--- /dev/null
+++ b/src/Misner_standard.c
@@ -0,0 +1,321 @@
+ /*@@
+   @file      Misner_standard.c
+   @date      March 1997
+   @author    Joan Masso
+   @desc 
+      Set up initial data for two Misner black holes
+   @enddesc 
+   @history
+   @hdate Sun Oct 17 11:05:48 1999 @hauthor Tom Goodale
+   @hdesc Converted to C
+   @endhistory
+ @@*/
+
+#include <stdio.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+#define SQR(a) ((a)*(a))
+
+ /*@@
+   @routine    Misner_standard
+   @date       
+   @author     Joan Masso, Ed Seidel
+   @desc 
+            Initialize the metric with a time symmetrical 
+            black hole spacetime containing
+            two axially symmetric misner black holes with a 
+            mass/length parameter mu. The mass is computed.
+            The spacetime line element has the form:
+                 $$ ds^2 = -dt^2 + \Psi^4 (dx^2+dy^2+dz^2) $$
+            and only $\Psi$ differs.
+            (Conformal factor from Karen Camarda)
+   @enddesc 
+   @calls     
+   @history 
+ 
+   @endhistory 
+
+   @par      mu
+   @pdesc    Misner parameter. 
+   @ptype    real
+   @pcomment Values less than 1.8 do not really correspond to two 
+            black holes, as there is an initial single event horizon 
+            surrounding the throats. So, with low values of mu we also
+            have distorted single black holes.
+   @endpar
+ 
+   @par      nmax
+   @pdesc    Summation limit for the misner series in the 'twobh' case.
+   @ptype    integer
+   @endpar
+
+@@*/
+void Misner_standard(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  int i,j,k;
+  int n;
+  CCTK_REAL zero, one, three;
+
+  CCTK_REAL csch, coth, r1, r2;
+
+  CCTK_REAL x_squared, y_squared;
+  CCTK_REAL r1_cubed, r2_cubed;
+  CCTK_REAL r1_5, r2_5;
+  CCTK_INT powfac;
+
+  int nx,ny,nz;
+  int index;
+  nx = cctk_lsh[0];
+  ny = cctk_lsh[1];
+  nz = cctk_lsh[2];
+
+
+  zero  = 0.0;
+  one   = 1.0;
+  three = 3.0;
+
+  /*     Initialize so we can accumulate
+   *     ------------------------------- 
+   */
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+        
+        psi[index] = one;
+      }
+    }
+  }
+
+  if (use_conformal_derivs == 1)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          psix[index]  = zero;
+          psiy[index]  = zero;
+          psiz[index]  = zero;
+          psixx[index] = zero;
+          psixy[index] = zero;
+          psixz[index] = zero;
+          psiyy[index] = zero;
+          psiyz[index] = zero;
+          psizz[index] = zero;
+        }
+      }
+    }
+  }
+
+  for(n = nmax; n >= 1; n--)
+  {
+    csch = one/sinh(mu*n);
+    coth = one/tanh(mu*n);
+    
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          x_squared = SQR(x[index]);
+          y_squared = SQR(y[index]);
+          r1 = sqrt(x_squared+y_squared+SQR(z[index]+coth));
+          r2 = sqrt(x_squared+y_squared+SQR(z[index]-coth));
+          
+          psi[index] += csch*(one/r1 + one/r2);
+          
+          if (use_conformal_derivs == 1)
+          {
+            r1_cubed = r1*r1*r1;
+            r2_cubed = r2*r2*r2;
+            r1_5     = pow(r1, 5);
+            r2_5     = pow(r2, 5);
+            psix[index]  +=  (-(x[index]/r2_cubed)-x[index]/r1_cubed)*csch;
+            psiy[index]  +=  (-(y[index]/r2_cubed)-y[index]/r1_cubed)*csch;
+            psiz[index]  +=  (-((z[index]-coth)/r2_cubed)-(z[index] + coth)/r1_cubed)*csch;
+            psixx[index] +=  ((three*x_squared)/r2_5-one/r2_cubed
+                              +(three*x_squared)/r1_5-one/r1_cubed)*csch;
+            psixy[index] +=  ((three*x[index]*y[index])/r2_5
+                              +(three*x[index]*y[index])/r1_5)*csch;
+            psixz[index] +=  ((three*x[index]*(z[index]-coth))/r2_5 
+                              +(three*x[index]*(z[index]+coth))/r1_5)*csch;
+            psiyy[index] +=  ((three*y_squared)/r2_5-one/r2_cubed 
+                              +(three*y_squared)/r1_5-one/r1_cubed)*csch;
+            psiyz[index] +=  ((three*y[index]*(z[index]-coth))/r2_5
+                              +(three*y[index]*(z[index]+coth))/r1_5)*csch;
+            psizz[index] += (-one/r2_cubed+(three*SQR(z[index]-coth))/r2_5+
+                             (three*SQR(z[index]+coth))/r1_5-one/r1_cubed)*csch;
+          }
+        }
+      }
+    }
+  }
+
+  /*     Cactus convention
+   *     -----------------
+   */
+
+  if (use_conformal_derivs == 1)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          psix[index]  /= psi[index];
+          psiy[index]  /= psi[index];
+          psiz[index]  /= psi[index];
+          psixx[index] /= psi[index];
+          psixy[index] /= psi[index];
+          psixz[index] /= psi[index];
+          psiyy[index] /= psi[index];
+          psiyz[index] /= psi[index];
+          psizz[index] /= psi[index];
+        }
+      }
+    }
+  }
+  
+  /*     compute the ADM mass
+   *     -------------------- 
+   */
+  mass = zero;
+  
+  for(n = 1; n <= nmax; n++)
+  {
+    mass += 4./sinh(n*mu);
+  }
+
+  printf("   >>ADM mass  %f\n",mass);
+  
+  /*     Should initialize lapse to Cadez value if possible
+   *     --------------------------------------------------
+   */
+
+  if (CCTK_Equals(initial_lapse,"cadez") )
+  { 
+    CCTK_INFO("Initialise with cadez lapse");
+    
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          x_squared = SQR(x[index]);
+          y_squared = SQR(y[index]);
+            
+          alp[index] = one;
+
+          powfac = 1;
+
+          for(n = 1; n <= nmax; n++)
+          {
+            coth = one/tanh(mu*n);
+            r1 = sqrt(x_squared+y_squared+SQR(z[index]+coth));
+            r2 = sqrt(x_squared+y_squared+SQR(z[index]-coth));
+            powfac *= -1;
+
+            alp[index] += powfac * one/sinh(mu*n)*(one/r1 + one/r2);
+          }
+          
+          alp[index] /= psi[index];
+        }
+      }
+    }
+  }
+  
+  /*     Metric depends on conformal state
+   *     ---------------------------------
+   */
+
+  if (*conformal_state == CONFORMAL_METRIC)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+
+          gxx[index] = one;
+          gyy[index] = one;
+          gzz[index] = one;
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+  else 
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = pow(psi[index],4);
+          gyy[index] = gxx[index];
+          gzz[index] = gxx[index];
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+
+  /*     Time-symmetric data
+   *     -------------------
+   */
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+        kxx[index] = zero;
+        kyy[index] = zero;
+        kzz[index] = zero;
+        kxy[index] = zero;
+        kxz[index] = zero;
+        kyz[index] = zero;
+      }
+    }
+  }
+
+  return;
+}
diff --git a/src/ParamChecker.F b/src/ParamChecker.F
deleted file mode 100644
index 3ed7497..0000000
--- a/src/ParamChecker.F
+++ /dev/null
@@ -1,79 +0,0 @@
- /*@@
-   @file      ParamChecker.F
-   @date      March 1999
-   @author    Gabrielle Allen
-   @desc 
-      Check parameters for black hole initial data and give some
-      information
-   @enddesc 
- @@*/
-
-#include "cctk.h"
-#include "cctk_arguments.h"
-#include "cctk_parameters.h"
-
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
- /*@@
-   @routine    ParamChecker
-   @date       March 1999
-   @author     Gabrielle Allen
-   @desc 
-      Check parameters for black hole initial data and give some
-      information
-   @enddesc 
-   @calls     
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-      subroutine ParamChecker(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer CCTK_Equals
-
-      if (CCTK_Equals(initial_data,"schwarzschild") == 1) then
-
-         write(*,*) 'One Black Hole:  Throat at :',mass/2d0
-
-      else if (CCTK_Equals(initial_data,"bl") == 1) then
-
-         write(*,*) "Setting up Brill Lindquist Solution"
-         write(*,*) "  Number of black holes (bl_nbh) :",bl_nbh
-         write(*,*) "  Black hole masses (bl_M_?) :",bl_M_1,bl_M_2,
-     &        bl_M_3,bl_M_4
-
-      else if (CCTK_Equals(initial_data,"multiple_misner") == 1) then
-
-         write(*,*) "Setting up Misner solution for multiple holes"
-
-      else if (CCTK_Equals(initial_data,"misner")==1) then
-
-         write(*,*) "Setting up Misner solution for two holes"
-         write(*,*) "  mu = ",mu
-
-      end if
-
-c     Remind users about the conformal metric
-c     ---------------------------------------
-      if (use_conformal == 1) then
-        write(*,*) "Implements conformal metric"
-        if (use_conformal_derivs == 1) then
-          write(*,*) "  ... and conformal derivatives"
-        else
-          write(*,*) "  ... but no conformal derivatives"
-        end if
-      else 
-        write(*,*) "Implements non-conformal metric" 
-        write(*,*) "  (Not usually a good idea!)"        
-      end if
-
-      return
-      end
-
diff --git a/src/ParamChecker.c b/src/ParamChecker.c
new file mode 100644
index 0000000..e55afc6
--- /dev/null
+++ b/src/ParamChecker.c
@@ -0,0 +1,81 @@
+ /*@@
+   @file      ParamChecker.F
+   @date      March 1999
+   @author    Gabrielle Allen
+   @desc 
+      Check parameters for black hole initial data and give some
+      information
+   @enddesc 
+ @@*/
+
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+ /*@@
+   @routine    ParamChecker
+   @date       March 1999
+   @author     Gabrielle Allen
+   @desc 
+      Check parameters for black hole initial data and give some
+      information
+   @enddesc 
+   @calls     
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void ParamChecker(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  if (CCTK_Equals(initial_data,"schwarzschild") == 1)
+  {
+    printf("One Black Hole:  Throat at : %f\n",mass/2.0);
+  }
+  else if (CCTK_Equals(initial_data,"bl") == 1)
+  {
+    printf("Setting up Brill Lindquist Solution\n");
+    printf("  Number of black holes (bl_nbh) : %d\n",bl_nbh);
+    printf("  Black hole masses (bl_M_?) : %f %f %f %f\n",bl_M_1,bl_M_2,bl_M_3,bl_M_4);
+  }
+  else if (CCTK_Equals(initial_data,"multiple_misner") == 1)
+  {
+    printf("Setting up Misner solution for multiple holes\n");
+  }
+  else if (CCTK_Equals(initial_data,"misner")==1)
+  {
+    printf("Setting up Misner solution for two holes\n");
+    printf("  mu = %f\n");
+  }
+
+  /*     Remind users about the conformal metric
+   *     ---------------------------------------
+   */
+  
+  if (use_conformal == 1)
+  {
+    printf("Implements conformal metric\n");
+    if (use_conformal_derivs == 1)
+    {
+      printf("  ... and conformal derivatives\n");
+    }
+    else
+    {
+      printf("  ... but no conformal derivatives");
+    }
+  }
+  else
+  { 
+    printf("Implements non-conformal metric\n");
+    printf("  (Not usually a good idea!)\n");
+  }        
+}
+
diff --git a/src/Schwarzschild.F b/src/Schwarzschild.F
deleted file mode 100644
index f8088d9..0000000
--- a/src/Schwarzschild.F
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "cctk.h"
-#include "cctk_arguments.h"
-#include "cctk_parameters.h"
-
-c Need include file from Einstein
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
-      subroutine Schwarzschild(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL zero,one,two,three
-
-      CCTK_REAL t1,r0
-
-      integer CCTK_Equals
-      integer i,j,k
-      integer is,js,ks,ie,je,ke
-      CCTK_REAL inval
-
-      zero = 0.d0
-      one = 1.d0
-      two = 2.d0
-      three = 3.d0
-
-c     conformal metric flag
-      if (use_conformal == 1) then
-
-         conformal_state = CONFORMAL_METRIC
-
-c        Compute conformal factor
-         psi = ( one + mass/two/r)
-
-c        derivatives of psi / psi
-         psix = -x*mass/two/r**3 / psi
-         psiy = -y*mass/two/r**3 / psi
-         psiz = -z*mass/two/r**3 / psi
-         
-         psixy = three*x*y*mass/two/r**5 /psi
-         psixz = three*x*z*mass/two/r**5 /psi
-         psiyz = three*y*z*mass/two/r**5 /psi
-         
-         psixx  = (three*x**2 - r**2)*mass/two/r**5/psi
-         psiyy  = (three*y**2 - r**2)*mass/two/r**5/psi
-         psizz  = (three*z**2 - r**2)*mass/two/r**5/psi          
-         
-         gxx = one
-         gyy = one
-         gzz = one
-         gxy = zero 
-         gxz = zero  
-         gyz = zero   
-          
-      else
-
-         conformal_state = NOCONFORMAL_METRIC
-
-         gxx = ( one + mass/two/r)**4
-         gyy = gxx
-         gzz = gxx
-         gxy = zero
-         gxz = zero
-         gyz = zero
-         
-      endif
-
-c     If the initial lapse is not one ...
-      if (CCTK_Equals(initial_lapse,"schwarz")==1) then 
-         print *,"Initialise with schwarzschild lapse"
-         alp = (2.*r - mass)/(2.*r+mass)
-      endif
-      
-c     time symmetric initial slice
-      kxx = zero 
-      kxy = zero
-      kxz = zero
-      kyy = zero
-      kyz = zero
-      kzz = zero
-
-      return
-      end subroutine Schwarzschild
diff --git a/src/Schwarzschild.c b/src/Schwarzschild.c
new file mode 100644
index 0000000..c4515dd
--- /dev/null
+++ b/src/Schwarzschild.c
@@ -0,0 +1,155 @@
+ /*@@
+   @file      Schwarzschild.c
+   @date      Sun Oct 17 10:35:41 1999
+   @author    Tom Goodale
+   @desc 
+   C version of Scwhwarzschild lapse routine
+   @enddesc 
+ @@*/
+
+#include <math.h>
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+
+/* Need include file from Einstein */
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+void Schwarzschild(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL zero,one,two,three;
+  CCTK_REAL tmp;
+  CCTK_REAL r_squared;
+
+  CCTK_REAL t1,r0;
+
+  int i,j,k;
+  int is,js,ks,ie,je,ke;
+  CCTK_REAL inval;
+  int nx,ny,nz;
+  int index;
+
+  zero = 0.0;
+  one = 1.0;
+  two = 2.0;
+  three = 3.0;
+
+  nx = cctk_lsh[0];
+  ny = cctk_lsh[1];
+  nz = cctk_lsh[2];
+  
+  /*     conformal metric flag */
+  if(use_conformal == 1)
+  {
+  
+    *conformal_state = CONFORMAL_METRIC;
+
+
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+
+          /*        Compute conformal factor */
+          psi[index] = ( one + mass/two/r[index]);
+         
+
+          /*        derivatives of psi / psi */
+
+          tmp = mass/two/pow(r[index],3) / psi[index];
+          psix[index] = -x[index]*tmp;
+          psiy[index] = -y[index]*tmp;
+          psiz[index] = -z[index]*tmp;
+         
+          tmp = mass/two/pow(r[index],5)/psi[index];
+          psixy[index] = three*x[index]*y[index]*tmp;
+          psixz[index] = three*x[index]*z[index]*tmp;
+          psiyz[index] = three*y[index]*z[index]*tmp;
+
+          r_squared = r[index]*r[index];
+          psixx[index]  = (three*x[index]*x[index] - r_squared)*tmp;
+          psiyy[index]  = (three*y[index]*y[index] - r_squared)*tmp;
+          psizz[index]  = (three*z[index]*z[index] - r_squared)*tmp;
+         
+          gxx[index] = one;
+          gyy[index] = one;
+          gzz[index] = one;
+          gxy[index] = zero; 
+          gxz[index] = zero;  
+          gyz[index] = zero;  
+        }
+      }
+    }
+  }     
+  else
+  {
+    *conformal_state = NOCONFORMAL_METRIC;
+
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+
+          gxx[index] = pow(( one + mass/two/r[index]), 4);
+          gyy[index] = gxx[index];
+          gzz[index] = gxx[index];
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+
+  /*     If the initial lapse is not one ... */
+  if (CCTK_Equals(initial_lapse,"schwarz"))
+  { 
+    CCTK_INFO("Initialise with schwarzschild lapse");
+
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+
+          alp[index] = (2.*r[index] - mass)/(2.*r[index]+mass);
+        }
+      }
+    }
+  }
+      
+  /*     time symmetric initial slice */
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+        
+        kxx[index] = zero;
+        kxy[index] = zero;
+        kxz[index] = zero;
+        kyy[index] = zero;
+        kyz[index] = zero;
+        kzz[index] = zero;
+      }
+    }
+  }
+    
+  return;
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index 789d02d..9519868 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,11 +2,11 @@
 # $Header$
 
 # Source files in this directory
-SRCS = 	ParamChecker.F\
-	Schwarzschild.F\
-	BrillLindquist.F\
-        Misner_standard.F\
-	Misner_multiple.F\
+SRCS = 	ParamChecker.c\
+	Schwarzschild.c\
+	BrillLindquist.c\
+        Misner_standard.c\
+	Misner_multiple.c\
 	Misner_points.c
 
 # Subdirectories containing source files