Tidied up and finished (hopefully) black hole initial data

Needs latex documentation now.


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

12 files changed, 648 insertions, 230 deletions

diff --git a/README b/README
index df60551..85dd220 100644
--- a/README
+++ b/README
@@ -1,20 +1,29 @@
 Cactus Code Thorn IDAnalyticBH
-Authors    : ...
-Managed by : ... <...@...........>
-Version    : ...  
+Authors    : Cactus Maintainers (originally written by Paul Walker and Joan Masso)
+Managed by : Cactus Maintainers <cactus-maint@aei-potsdam.mpg.de>
+Version    : 1.0
 CVS info   : $Header$
 --------------------------------------------------------------------------
 
 1. Purpose of the thorn
 
-This thorn does ...
+This thorn calculates analytic initial data for the Einstein grid functions
+(lapse, shift, metric, curv) for various black hole spacetimes:
+
+	Schwarzschild black hole
+ 	2 Misner black holes placed on the z-axis
+	n Misner black holes placed in a circle in the x-y plane
+	Brill-Lindquist data
 
 2. Dependencies of the thorn
 
-This thorn additionally requires thorns ...
+This thorn additionally requires thorns Einstein, which defines
+the Grid Functions.
 
 3. Thorn distribution
 
-This thorn is available to ...
+This thorn is publically available.
 
 4. Additional information
+
+
diff --git a/interface.ccl b/interface.ccl
index fdfbf53..d274e8a 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -4,3 +4,10 @@
 implements: idanalyticbh
 inherits: einstein
 
+real misner_workspace TYPE = GF
+{
+  csch,
+  coth,
+  r1,
+  r2
+} "Workspace arrays for Misner_standard"
diff --git a/param.ccl b/param.ccl
index 6b6c4b6..76f25c5 100644
--- a/param.ccl
+++ b/param.ccl
@@ -3,14 +3,16 @@
 
 private:
 
+# Schwarzschild parameters
+# ------------------------
 REAL mass "Mass of black hole"
 {
 : :: "Not sure if it can be negative or not"
 } 2.0
 
 
-# Misner Parameters
-# -----------------
+# Multiple Misner Parameters
+# --------------------------
 REAL mu "Misner mu value"
 {
 0: :: ""
@@ -100,15 +102,12 @@ REAL bl_M_4 "Mass of 4th BL hole"
 
 friend:einstein
 
-EXTENDS LOGICAL use_conformal ""
-{ 
-}
-
 EXTENDS KEYWORD initial_data ""
 {
-  "bl_bh"         :: "Brill Lindquist black holes"
-  "misner_bh"     :: "Misner black holes"
-  "schwarzschild" :: "One Schwarzshild black hole"
+  "schwarzschild"      :: "One Schwarzshild black hole"
+  "bl_bh"              :: "Brill Lindquist black holes"
+  "misner_bh"          :: "Misner black holes"
+  "multiple_misner_bh" :: "Multiple Misner black holes"
 } 
 
 EXTENDS KEYWORD lapseinit ""
diff --git a/schedule.ccl b/schedule.ccl
index 6812c02..a2b38fe 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,10 +4,6 @@
 
 if (CCTK_Equals(initial_data,"schwarzschild")) 
 { 
-   if (use_conformal)
-   {
-     STORAGE: confac, confac_1derivs, confac_2derivs
-   }
 	
    schedule Schwarzschild at CACTUS_INITIAL
    {
@@ -18,8 +14,7 @@ if (CCTK_Equals(initial_data,"schwarzschild"))
 
 if (CCTK_Equals(initial_data,"bl_bh")) 
 { 
-   STORAGE: confac, confac_1derivs, confac_2derivs
-
+   STORAGE: confac
    schedule BrillLindquist at CACTUS_INITIAL
    {
      LANG: Fortran
@@ -29,12 +24,20 @@ if (CCTK_Equals(initial_data,"bl_bh"))
 
 if (CCTK_Equals(initial_data,"misner_bh")) 
 { 
-   STORAGE: confac, confac_1derivs, confac_2derivs
-
-   schedule Misner at CACTUS_INITIAL
+   STORAGE: confac
+   schedule Misner_standard at CACTUS_INITIAL
    {
      LANG: Fortran
-   } "Construct initial data for Misner black holes"
+   } "Construct initial data for multiple Misner black holes"
 
 }
  
+if (CCTK_Equals(initial_data,"multiple_misner_bh")) 
+{ 
+   STORAGE: confac
+   schedule Misner_multiple at CACTUS_INITIAL
+   {
+     LANG: Fortran
+   } "Construct initial data for two Misner black holes"
+
+}
diff --git a/src/BrillLindquist.F b/src/BrillLindquist.F
index 38db1d8..287e531 100644
--- a/src/BrillLindquist.F
+++ b/src/BrillLindquist.F
@@ -1,8 +1,32 @@
+ /*@@
+   @file      BrillLindquist.F
+   @date      
+   @author    
+   @desc 
+      Set up initial data for Brill Lindquist black holes
+   @enddesc 
+ @@*/
+
 #include "cctk.h"
 #include "declare_arguments.h"
 #include "declare_parameters.h"
 
-#include "../../Einstein/src/Einstein.h"
+#include "../../../packages/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)
 
@@ -12,19 +36,11 @@
       DECLARE_PARAMETERS
 
       integer i
-      REAL hole_mass(4)
-      REAL hole_x0(4)
-      REAL hole_y0(4)
-      REAL hole_z0(4)
-
-c     Must use conformal metric
-c     -------------------------
-      if (conformal_state == NOCONFORMAL_METRIC) then
-         call CCTK_Warn(1,"Changing to conformal metric in BrillLindquist")
-         conformal_state = CONFORMAL_METRIC
-      end if
+      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
@@ -46,28 +62,39 @@ c     Put parameters into arrays for following calculations
       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
-      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
-
-      do i=1,bl_nbh
+
+      if (use_conformal_derivs) 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) 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
@@ -107,36 +134,57 @@ c     Maple Output
      &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
-c     End Maple Output
+
+      end if
 
       end do
 
-C     Cactus conventions:
-      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
-
-C     The conformal metric is flat.
-      gxx = 1d0
-      gyy = 1d0
-      gzz = 1d0
-      gxy = 1d0 
-      gxz = 1d0  
-      gyz = 1d0   
-
-C     Time-symmetric data.
-      hxx = 1d0
-      hyy = 1d0
-      hzz = 1d0
-      hxy = 1d0 
-      hxz = 1d0  
-      hyz = 1d0   
+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     -------------------
+      hxx = 0d0
+      hyy = 0d0
+      hzz = 0d0
+      hxy = 0d0 
+      hxz = 0d0  
+      hyz = 0d0   
 
       return
       end
diff --git a/src/Misner.F b/src/Misner.F
deleted file mode 100644
index 532c1ce..0000000
--- a/src/Misner.F
+++ /dev/null
@@ -1,109 +0,0 @@
-C     This wrapper written by Carsten Gundlach. 
-
-#include "cctk.h"
-#include "declare_arguments.h"
-#include "declare_parameters.h"
-
-#include "../../Einstein/src/Einstein.h"
-
-      subroutine Misner(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_PARAMETERS
-
-      integer i, j, k
-      REAL one, zero, tmp0, tmp1, tmp2, tmp3, tmp4, nm_eps
-
-C     Finite differencing step.
-      parameter(nm_eps = 1.d-6)
-
-      one = 1.0D0
-      zero = 0.0D0
-
-c     Must use conformal metric
-c     -------------------------
-      if (conformal_state == NOCONFORMAL_METRIC) then
-         call CCTK_Warn(1,"Changing to conformal metric in Misner")
-         conformal_state = CONFORMAL_METRIC
-      end if
-
-C     Initialize some C global variables before calls to 
-C     misner_psi_eval.
-      call nmisner_init(misner_nbh,mu,nmax)
-
-C     Get value of psi pointwise from a C function.
-      do k=1,sh(1)
-         do j=1,sh(2)
-            do i=1,sh(3)
- 
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k),z(i,j,k),tmp0)
-               psi(i,j,k) = tmp0
-               
-               call nmisner_eval_psi(x(i,j,k)+nm_eps,y(i,j,k),z(i,j,k),tmp1)
-               call nmisner_eval_psi(x(i,j,k)-nm_eps,y(i,j,k),z(i,j,k),tmp2)
-               psix(i,j,k) = 0.5d0 * (tmp1 - tmp2) / nm_eps
-               psixx(i,j,k) = (tmp1 + tmp2 - 2.d0 * tmp0) / nm_eps**2
-
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)+nm_eps,z(i,j,k),tmp1)
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)-nm_eps,z(i,j,k),tmp2)
-               psiy(i,j,k) = 0.5d0 * (tmp1 - tmp2) / nm_eps
-               psiyy(i,j,k) = (tmp1 + tmp2 - 2.d0 * tmp0) / nm_eps**2
-
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k),z(i,j,k)+nm_eps,tmp1)
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k),z(i,j,k)-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 nmisner_eval_psi(x(i,j,k)+nm_eps,y(i,j,k)+nm_eps,z(i,j,k),tmp1)
-               call nmisner_eval_psi(x(i,j,k)+nm_eps,y(i,j,k)-nm_eps,z(i,j,k),tmp2)
-               call nmisner_eval_psi(x(i,j,k)-nm_eps,y(i,j,k)+nm_eps,z(i,j,k),tmp3)
-               call nmisner_eval_psi(x(i,j,k)-nm_eps,y(i,j,k)-nm_eps,z(i,j,k),tmp4)
-               psixy(i,j,k) = 0.25d0 * (tmp1 - tmp2 - tmp3 + tmp4) / nm_eps**2
-
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)+nm_eps,z(i,j,k)+nm_eps,tmp1)
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)-nm_eps,z(i,j,k)+nm_eps,tmp2)
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)+nm_eps,z(i,j,k)-nm_eps,tmp3)
-               call nmisner_eval_psi(x(i,j,k),y(i,j,k)-nm_eps,z(i,j,k)-nm_eps,tmp4)
-               psiyz(i,j,k) = 0.25d0 * (tmp1 - tmp2 - tmp3 + tmp4) / nm_eps**2
-
-               call nmisner_eval_psi(x(i,j,k)+nm_eps,y(i,j,k),z(i,j,k)+nm_eps,tmp1)
-               call nmisner_eval_psi(x(i,j,k)+nm_eps,y(i,j,k),z(i,j,k)-nm_eps,tmp2)
-               call nmisner_eval_psi(x(i,j,k)-nm_eps,y(i,j,k),z(i,j,k)+nm_eps,tmp3)
-               call nmisner_eval_psi(x(i,j,k)-nm_eps,y(i,j,k),z(i,j,k)-nm_eps,tmp4)
-               psixz(i,j,k) = 0.25d0 * (tmp1 - tmp2 - tmp3 + tmp4) / nm_eps**2
-            end do
-         end do
-      end do
-
-C     Cactus conventions:
-      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
-
-C     The conformal metric is flat.
-      gxx = one
-      gyy = one
-      gzz = one
-      gxy = zero 
-      gxz = zero  
-      gyz = zero   
-
-C     Time-symmetric data.
-      hxx = zero
-      hyy = zero
-      hzz = zero
-      hxy = zero 
-      hxz = zero  
-      hyz = zero   
-
-      return
-      end
-
diff --git a/src/Misner_multiple.F b/src/Misner_multiple.F
new file mode 100644
index 0000000..04e7348
--- /dev/null
+++ b/src/Misner_multiple.F
@@ -0,0 +1,153 @@
+ /*@@
+   @file      Misner_multiple.F
+   @date      
+   @author    Carsten Gundlach
+   @desc 
+      Set up initial data for multiple Misner black holes
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "declare_arguments.h"
+#include "declare_parameters.h"
+
+#include "../../../packages/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_PARAMETERS
+
+      integer :: i, j, k
+      REAL :: xval, yval, zval
+      REAL :: one, zero, tmp0, tmp1, tmp2, tmp3, tmp4
+      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,sh(1)
+         do j=1,sh(2)
+            do i=1,sh(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) 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     -------------------
+      hxx = zero
+      hyy = zero
+      hzz = zero
+      hxy = zero 
+      hxz = zero  
+      hyz = zero   
+
+      return
+      end
+
diff --git a/src/Misner_points.c b/src/Misner_points.c
index d0da40d..36964c1 100644
--- a/src/Misner_points.c
+++ b/src/Misner_points.c
@@ -1,13 +1,15 @@
-/* This code written by Steve Brandt. */
-
-/* This calculates the conformal factor for nbh black holes, 
-with naked mass m0 = 2 csch(mu) each, and placed on a circle in the
-xy plane around the origin, of radius coth(mu). 
-One of them sits on the positive x axis, the others are evenly spaced.
-Naked mass here corresponds to the term m0 / (2 |r - r0|) in the expansion.
-For nbh =2, these are the same masses as in the routine misner.F, 
-with the difference that here two BHs sit on the positive and negative
-x axis, there on the z-axis. */
+ /*@@
+   @file      Misner_points.c
+   @date      
+   @author    Steve Brandt
+   @desc 
+      	      This calculates the conformal factor for nbh black holes, 
+              with naked mass m0 = 2 csch(mu) each, and placed on a circle in the
+              xy plane around the origin, of radius coth(mu). 
+              One of them sits on the positive x axis, the others are evenly spaced.
+              Naked mass here corresponds to the term m0 / (2 |r - r0|) in the expansion.
+   @enddesc 
+ @@*/
 
 #include "cctk.h"
 
@@ -16,12 +18,12 @@ x axis, there on the z-axis. */
 #include <math.h>
 
 int nbholes;
-double mu;
+Double mu;
 
 /* Basic data about a brill-lindquist black hole term. */
 struct bhole {
-  double x,y;
-  double mass;
+  Double x,y;
+  Double mass;
 
   /* i gives either the number of the seed
      black hole we are starting with, or
@@ -35,16 +37,30 @@ struct bhole {
 /* The seed black holes. */
 struct bhole bholes[MAXBHOLES];
 
-double csch(double mu) {
+Double csch(Double mu) {
   return 1.0/sinh(mu); 
 }
-double coth(double mu) {
+Double coth(Double mu) {
   return cosh(mu)/sinh(mu); 
 }
 
-/* Isometrize black hole a1 through hole a2 */
-void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3) {
-  double rad,radtwo;
+ /*@@
+   @routine    fill_iso
+   @date       
+   @author     Steve Brandt
+   @desc 
+     	       Isometrize black hole a1 through hole a2
+   @enddesc 
+   @calls      fill_iso   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3) 
+{
+  Double rad,radtwo;
   radtwo=(
       (a1->x - a2->x)*(a1->x - a2->x)+
       (a1->y - a2->y)*(a1->y - a2->y)
@@ -55,17 +71,33 @@ void iso(struct bhole *a1, struct bhole *a2, struct bhole *a3) {
   a3->y = a2->y+(a2->mass*a2->mass)*(a1->y - a2->y)/radtwo;
 }
 
-/* Fills in the iso structure of a given black hole.  Applies
-   recursively to the number of terms desired. */
-void fill_iso(struct bhole *b,int n) {
+ /*@@
+   @routine    fill_iso
+   @date       
+   @author     Steve Brandt
+   @desc 
+     	       Fills in the iso structure of a given black hole.  
+               Applies recursively to the number of terms desired.
+   @enddesc 
+   @calls      fill_iso   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void fill_iso(struct bhole *b, int n) 
+{
   int i,j;
-  if(n==0) {
+  if(n==0) 
+  {
     b->isos = 0;
     return;
   }
   b->isos = (struct bhole *)malloc(sizeof(struct bhole)*(nbholes-1));
   assert(b->isos != 0);
-  for(j=0, i=0;i<nbholes;i++) {
+  for(j=0, i=0;i<nbholes;i++) 
+  {
     if(i != b->i) {
       iso(b,&bholes[i],&b->isos[j]);
       b->isos[j].i = i;
@@ -75,32 +107,70 @@ void fill_iso(struct bhole *b,int n) {
   }
 }
 
-/* Initializes the black holes, then makes the
-   isometry black holes. */
-void FORTRAN_NAME(nmisner_init)(int *n, double *mu,int *terms) 
+ /*@@
+   @routine    Misner_init
+   @date       
+   @author     Steve Brandt
+   @desc 
+      Initialises the black holes then makes the isometry black holes
+   @enddesc 
+   @calls     
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void FORTRAN_NAME(Misner_init)(int *n, Double *mu, int *terms) 
 {
+
   int i;
-  double pi,ang;
+  Double pi,ang;
+
   assert((nbholes=*n) < MAXBHOLES);
+
   pi = 4.0*atan(1.);
+
   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].i = i;
     bholes[i].isos = 0;
   }
+
   for(i=0;i<*n;i++)
     fill_iso(&bholes[i],*terms);
+
 }
 
-double eval_bh_psi(struct bhole *b,double x,double y,double z) {
+
+ /*@@
+   @routine    eval_bh_psi
+   @date       
+   @author     Steve Brandt
+   @desc 
+      
+   @enddesc 
+   @calls      eval_bh_psi   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+Double eval_bh_psi(struct bhole *b, Double x, Double y, Double z) 
+{
   int i;
-  double res;
+  Double res;
   res = 0.0;
-  if(b->isos != 0) {
-    for(i=0;i<nbholes-1;i++) {
+  if(b->isos != 0) 
+  {
+    for(i=0;i<nbholes-1;i++) 
+    {
       res += eval_bh_psi(&b->isos[i],x,y,z);
     }
   }
@@ -112,8 +182,23 @@ double eval_bh_psi(struct bhole *b,double x,double y,double z) {
   return res;
 }
 
-/* use this function to evaluate psi at a point. */
-void FORTRAN_NAME(nmisner_eval_psi)(double *x,double *y,double *z,double *res) {
+
+ /*@@
+   @routine    MisnerEvalPsi
+   @date       
+   @author     Steve Brandt
+   @desc 
+      Evaluate psi at a point
+   @enddesc 
+   @calls      eval_bh_psi
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void FORTRAN_NAME(MisnerEvalPsi)(Double *x, Double *y, Double *z, Double *res) 
+{
   int i;
   *res = 1;
   for(i=0;i<nbholes;i++)
diff --git a/src/Misner_standard.F b/src/Misner_standard.F
new file mode 100644
index 0000000..08d9c67
--- /dev/null
+++ b/src/Misner_standard.F
@@ -0,0 +1,177 @@
+ /*@@
+   @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 "declare_arguments.h"
+#include "declare_parameters.h"
+
+#include "../../../packages/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_PARAMETERS
+
+      integer :: n
+      INTEGER :: CCTK_Equals
+      REAL :: zero, one, three
+
+      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) 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
+      write(*,*) "   >>ADM mass",mass
+
+c     Should initialize lapse to Cadez value if possible
+c     --------------------------------------------------
+      if (CCTK_Equals(lapseinit,"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     -------------------
+      hxx = zero
+      hyy = zero
+      hzz = zero
+      hxy = zero 
+      hxz = zero  
+      hyz = zero   
+
+      return
+      end
diff --git a/src/ParamChecker.F b/src/ParamChecker.F
index 10c9caa..b4e7e5e 100644
--- a/src/ParamChecker.F
+++ b/src/ParamChecker.F
@@ -1,7 +1,34 @@
+ /*@@
+   @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 "declare_arguments.h"
 #include "declare_parameters.h"
 
+#include "../../../packages/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
@@ -12,21 +39,39 @@
       integer CCTK_Equals
 
       if (CCTK_Equals(initial_data,"schwarzschild")==1) then
+
          write(*,*) 'One Black Hole:  Throat at :',mass/2d0
-         if (use_conformal == 1) then
-            write(*,*) "Uses conformal metric"
-         else 
-            write(*,*) "Uses non-conformal metric"         
-         end if
+
       else if (CCTK_Equals(initial_data,"bl")==1) then
+
          write(*,*) "Setting up Brill Lindquist Solution"
-         write(*,*) "   Uses conformal metric"
-         write(*,*) "   Number of black holes (bl_nbh) :",bl_nbh
-         write(*,*) "   Black hole masses (bl_M_?) :",bl_M_1,bl_M_2,
+         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"
-         write(*,*) "   Uses conformal metric"
+
+         write(*,*) "Setting up Misner solution for two holes"
+         write(*,*) "  mu = ",mu
+
+      end if
+
+c     Remind users about the conformal metric
+c     ---------------------------------------
+      if (conformal_state == CONFORMAL_METRIC) 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
diff --git a/src/Schwarzschild.F b/src/Schwarzschild.F
index d597e45..ec10c59 100644
--- a/src/Schwarzschild.F
+++ b/src/Schwarzschild.F
@@ -25,7 +25,7 @@ c Need include file from Einstein
       one = 1.d0
       two = 2.d0
       three = 3.d0
-            
+
 c     conformal metric flag
       if (use_conformal == 1) then
 
diff --git a/src/make.code.defn b/src/make.code.defn
index 2bdae0b..789d02d 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -5,7 +5,8 @@
 SRCS = 	ParamChecker.F\
 	Schwarzschild.F\
 	BrillLindquist.F\
-	Misner.F\
+        Misner_standard.F\
+	Misner_multiple.F\
 	Misner_points.c
 
 # Subdirectories containing source files