Add the capability to rotate initial data, similar to the existing

boost transformation.

Add spherical Kerr-Schild initial data, i.e., Kerr-Schild data in
which the horizon is a coordinate sphere.

Add the capability to smooth Kerr-Schild data with a parabolic term.

Re-indent param.ccl consistently.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@246 e296648e-0e4f-0410-bd07-d597d9acff87

9 files changed, 597 insertions, 234 deletions

diff --git a/src/boost.F77 b/src/boost.F77
index d81b100..dc26340 100644
--- a/src/boost.F77
+++ b/src/boost.F77
@@ -1,5 +1,7 @@
 c This subroutine calculates the 4-metric and its inverse at an event,
-c taking into account an optional Lorentz boost.
+c taking into account an optional Lorentz boost and an optional rotation.
+c The model is first rotated and then boosted, such that the boost is
+c applied to the rotated model.
 c $Header$
 c
 c The coordinates are
@@ -9,6 +11,12 @@ c The 4-metrics are
 c   Cgdd(a,b) = Cactus $g_{ab}$            Cguu(a,b) = Cactus $g^{ab}$
 c   Mgdd(a,b) = Model  $g_{ab}$            Mguu(a,b) = Model  $g^{ab}$
 c
+c For a definition of the Euler angles in the conventions used below, see
+c    http://mathworld.wolfram.com/EulerAngles.html
+c Another useful resource may be
+c    http://en.wikipedia.org/wiki/Euler_angles
+c although this uses (on 2006-11-29) different conventions.
+c
 c This file is copyright (c) 2003 by Jonathan Thornburg <jthorn@aei.mpg.de>.
 c This file is covered by the GNU GPL license; see the files ../README
 c and ../COPYING for details. 
@@ -21,7 +29,7 @@ c
 
 #include "param_defs.inc"
 
-	subroutine Exact__metric(
+      subroutine Exact__metric(
      $     decoded_exact_model,
      $     x, y, z, t,
      $     gdtt, gdtx, gdty, gdtz,
@@ -30,83 +38,94 @@ c
      $     guxx, guyy, guzz, guxy, guyz, guxz,
      $     psi)
 
-	implicit none
-	DECLARE_CCTK_FUNCTIONS
-	DECLARE_CCTK_PARAMETERS
+      implicit none
+      DECLARE_CCTK_FUNCTIONS
+      DECLARE_CCTK_PARAMETERS
 
 c input arguments
-	CCTK_INT decoded_exact_model
-	CCTK_REAL x, y, z, t
+      CCTK_INT decoded_exact_model
+      CCTK_REAL x, y, z, t
 
 c output arguments
-	CCTK_REAL gdtt, gdtx, gdty, gdtz,
-     $		  gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
-     $		  gutt, gutx, guty, gutz,
-     $		  guxx, guyy, guzz, guxy, guyz, guxz,
-     $            psi
+      CCTK_REAL gdtt, gdtx, gdty, gdtz,
+     $          gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $          gutt, gutx, guty, gutz,
+     $          guxx, guyy, guzz, guxy, guyz, guxz,
+     $          psi
 
 c intrinsic functions called
-	CCTK_REAL sqrt
+      CCTK_REAL sqrt
 
 c static local variables describing Lorentz transformation
-	logical   firstcall
-        data      firstcall /.true./
-	CCTK_REAL gamma
-	CCTK_REAL vv(3), nn(3)
-	CCTK_REAL parallel(3,3), perp(3,3)
-	CCTK_REAL Cx_par(3), Cx_perp(3)
-	CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3)
-	CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3)
-        save      firstcall
-        save      gamma
-        save      vv, nn
-        save      parallel, perp
-	save      Cx_par, Cx_perp
-        save      partial_Mx_wrt_Cx
-        save      partial_Cx_wrt_Mx
+      logical   firstcall
+      data      firstcall /.true./
+      CCTK_REAL gamma
+      CCTK_REAL vv(3), nn(3)
+      CCTK_REAL parallel(3,3), perp(3,3)
+      CCTK_REAL Cx_par(3), Cx_perp(3)
+      CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3)
+      CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3)
+      CCTK_REAL R(0:3,0:3)
+      save      firstcall
+      save      gamma
+      save      vv, nn
+      save      parallel, perp
+      save      Cx_par, Cx_perp
+      save      partial_Mx_wrt_Cx
+      save      partial_Cx_wrt_Mx
+      save      R
 
 c coordinates and 4-metric
-	CCTK_REAL Ct, Cx(3)
-	CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3)
-	CCTK_REAL Mt, Mx(3)
-	CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3)
+      CCTK_REAL Cx(0:3)
+      CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3)
+      CCTK_REAL Mx(0:3)
+      CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3)
+      CCTK_REAL Nx(0:3)
+      CCTK_REAL Ngdd(0:3,0:3), Nguu(0:3,0:3)
 
 c misc temps
-	CCTK_REAL vnorm, vnormsq
-	CCTK_REAL delta_ij
-	CCTK_REAL Cx_par_i, Cx_perp_i
-	CCTK_REAL vdotCx
-	CCTK_REAL Cgdd_ab, Cguu_ab
-	character*100 warn_buffer
+      CCTK_REAL vnorm, vnormsq
+      CCTK_REAL delta_ij
+      CCTK_REAL Cx_par_i, Cx_perp_i
+      CCTK_REAL vdotCx
+      CCTK_REAL Cgdd_ab, Cguu_ab
+      CCTK_REAL cos_phi, sin_phi
+      CCTK_REAL cos_theta, sin_theta
+      CCTK_REAL cos_psi, sin_psi
+      CCTK_REAL R_phi(0:3,0:3), R_theta(0:3,0:3), R_psi(0:3,0:3)
+      character*1000 warn_buffer
 
 c flags, array indices, etc
-	logical Tmunu_flag
-	integer i, j
-	integer Ca, Cb, MA, MB
+      logical Tmunu_flag
+      integer i, j, k, l
+      integer Ca, Cb, MA, MB
 
 c constants
-	integer n
-	parameter (n = 3)
+      integer n
+      parameter (n = 3)
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
 c
-c optimized fast-path if no Lorentz boost
+c optimized fast-path if no Lorentz boost and no rotation
 c
-	if (       (boost_vx .eq. 0.0)
-     $	     .and. (boost_vy .eq. 0.0)
-     $	     .and. (boost_vz .eq. 0.0)       ) then
-		call Exact__metric_for_model(
-     $			decoded_exact_model,
-     $			x, y, z, t,
-     $			gdtt, gdtx, gdty, gdtz,
-     $			gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
-     $			gutt, gutx, guty, gutz,
-     $			guxx, guyy, guzz, guxy, guyz, guxz,
-     $			psi,
-     $			Tmunu_flag)
-		return
-	endif
+      if (       (boost_vx .eq. 0.0)
+     $     .and. (boost_vy .eq. 0.0)
+     $     .and. (boost_vz .eq. 0.0)
+     $     .and. (rotation_euler_phi   .eq. 0.0)
+     $     .and. (rotation_euler_theta .eq. 0.0)
+     $     .and. (rotation_euler_psi   .eq. 0.0)) then
+         call Exact__metric_for_model(
+     $        decoded_exact_model,
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz,
+     $        psi,
+     $        Tmunu_flag)
+         return
+      end if
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
@@ -125,203 +144,321 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c
 c compute Lorentz transformation information on first call
 c
-	if (firstcall) then
-		firstcall = .false.
+      if (firstcall) then
+         firstcall = .false.
 
 c boost velocity
-		vv(1) = boost_vx
-		vv(2) = boost_vy
-		vv(3) = boost_vz
+         vv(1) = boost_vx
+         vv(2) = boost_vy
+         vv(3) = boost_vz
 
 c Lorentz gamma factor, unit vector in direction of boost velocity
-		vnormsq = 0.0d0
-			do 100 i = 1,n
-			vnormsq = vnormsq + vv(i)*vv(i)
-100			continue
-		gamma = 1.0 / sqrt(1.0 - vnormsq)
-		vnorm = sqrt(vnormsq)
-			do 110 i = 1,n
-			nn(i) = vv(i) / vnorm
-110			continue
+         vnormsq = 0.0
+         do i = 1,n
+            vnormsq = vnormsq + vv(i)*vv(i)
+         end do
+         gamma = 1.0 / sqrt(1.0 - vnormsq)
+         vnorm = sqrt(vnormsq)
+         if (       (boost_vx .eq. 0.0)
+     $        .and. (boost_vy .eq. 0.0)
+     $        .and. (boost_vz .eq. 0.0)) then
+            nn(1) = 1
+            nn(2) = 0
+            nn(3) = 0
+         else
+            do i = 1,n
+               nn(i) = vv(i) / vnorm
+            end do
+         end if
 
 c projection operators parallel(*,*) and perp(*,*)
-			do 210 j = 1,n
-				do 200 i = 1,n
-				parallel(i,j) = nn(i) * nn(j)
-				if (i .eq. j) then
-					delta_ij = 1.0d0
-				   else
-					delta_ij = 0.0d0
-				endif
-				perp(i,j) = delta_ij - parallel(i,j)
-200				continue
-210			continue
+         do j = 1,n
+            do i = 1,n
+               parallel(i,j) = nn(i) * nn(j)
+               if (i .eq. j) then
+                  delta_ij = 1.0
+               else
+                  delta_ij = 0.0
+               end if
+               perp(i,j) = delta_ij - parallel(i,j)
+            end do
+         end do
 
 c partial derivatives of Model coordinates with respect to Cactus coordinates
-		partial_Mx_wrt_Cx(0,0) = gamma
-			do 300 i = 1,n
-			partial_Mx_wrt_Cx(0,i) = -gamma*vv(i)
-300			continue
-			do 320 i = 1,n
-			partial_Mx_wrt_Cx(i,0) = -gamma*vv(i)
-				do 310 j=1,n
-				partial_Mx_wrt_Cx(i,j)
-     $					= gamma*parallel(i,j)
-     $					  + perp(i,j)
-310				continue
-320			continue
+         partial_Mx_wrt_Cx(0,0) = gamma
+         do i = 1,n
+            partial_Mx_wrt_Cx(0,i) = -gamma*vv(i)
+         end do
+         do i = 1,n
+            partial_Mx_wrt_Cx(i,0) = -gamma*vv(i)
+            do j=1,n
+               partial_Mx_wrt_Cx(i,j) = gamma*parallel(i,j) + perp(i,j)
+            end do
+         end do
 
 c partial derivatives of Cactus coordinates with respect to Model coordinates
-		partial_Cx_wrt_Mx(0,0) = gamma
-			do 400 i = 1,n
-			partial_Cx_wrt_Mx(0,i) = + gamma*vv(i)
-400			continue
-			do 420 i = 1,n
-			partial_Cx_wrt_Mx(i,0) = + gamma*vv(i)
-				do 410 j=1,n
-				partial_Cx_wrt_Mx(i,j)
-     $					= gamma*parallel(i,j)
-     $					  + perp(i,j)
-410				continue
-420			continue
-
-	endif
+         partial_Cx_wrt_Mx(0,0) = gamma
+         do i = 1,n
+            partial_Cx_wrt_Mx(0,i) = + gamma*vv(i)
+         end do
+         do i = 1,n
+            partial_Cx_wrt_Mx(i,0) = + gamma*vv(i)
+            do j=1,n
+               partial_Cx_wrt_Mx(i,j) = gamma*parallel(i,j) + perp(i,j)
+            end do
+         end do
+
+c Sines and cosines of rotation angles
+         cos_phi   = cos (rotation_euler_phi)
+         sin_phi   = sin (rotation_euler_phi)
+         cos_theta = cos (rotation_euler_theta)
+         sin_theta = sin (rotation_euler_theta)
+         cos_psi   = cos (rotation_euler_psi)
+         sin_psi   = sin (rotation_euler_psi)
+
+c Set up individual rotation matrices
+         R_phi(0,0) = 1
+         R_phi(0,1) = 0
+         R_phi(0,2) = 0
+         R_phi(0,3) = 0
+         R_phi(1,0) = 0
+         R_phi(1,1) = + cos_phi
+         R_phi(1,2) = + sin_phi
+         R_phi(1,3) = 0
+         R_phi(2,0) = 0
+         R_phi(2,1) = - sin_phi
+         R_phi(2,2) = + cos_phi
+         R_phi(2,3) = 0
+         R_phi(3,0) = 0
+         R_phi(3,1) = 0
+         R_phi(3,2) = 0
+         R_phi(3,3) = 1
+
+         R_theta(0,0) = 1
+         R_theta(0,1) = 0
+         R_theta(0,2) = 0
+         R_theta(0,3) = 0
+         R_theta(1,0) = 0
+         R_theta(1,1) = 1
+         R_theta(1,2) = 0
+         R_theta(1,3) = 0
+         R_theta(2,0) = 0
+         R_theta(2,1) = 0
+         R_theta(2,2) = + cos_theta
+         R_theta(2,3) = + sin_theta
+         R_theta(3,0) = 0
+         R_theta(3,1) = 0
+         R_theta(3,2) = - sin_theta
+         R_theta(3,3) = + cos_theta
+
+         R_psi(0,0) = 1
+         R_psi(0,1) = 0
+         R_psi(0,2) = 0
+         R_psi(0,3) = 0
+         R_psi(1,0) = 0
+         R_psi(1,1) = + cos_psi
+         R_psi(1,2) = + sin_psi
+         R_psi(1,3) = 0
+         R_psi(2,0) = 0
+         R_psi(2,1) = - sin_psi
+         R_psi(2,2) = + cos_psi
+         R_psi(2,3) = 0
+         R_psi(3,0) = 0
+         R_psi(3,1) = 0
+         R_psi(3,2) = 0
+         R_psi(3,3) = 1
+
+c Combine individual rotation matrices
+         do i = 0,n
+            do j = 0,n
+               R(i,j) = 0
+               do k = 0,n
+                  do l = 0,n
+                     R(i,j) = R(i,j) + R_psi(i,k) * R_theta(k,l) * R_phi(l,j)
+                  end do
+               end do
+            end do
+         end do
+
+c     Notes that help me (Erik Schnetter) think:
+c     This considers a rotation with phi=0, theta=pi/2, psi=0.
+c     Nx(i) = Nx(i) + R(j,i) * Mx(j)
+c     Nx(1) = - Mx(3)
+c     Nx(3) =   Mx(1)
+c     Mgxx(1,1) =   Ngxx(3,3)
+c     Mgxx(1,3) = - Ngxx(1,3)
+c     Mgxx(3,3) =   Ngxx(1,1)
+c     Mgxx(i,j) = R(i,k) R(j,l) Ngxx(k,l)  [correct]
+c     Mgxx(i,j) = R(k,i) R(l,j) Ngxx(k,l)  [correct]
+c     Mbetax(1) = - Nbetax(3)
+c     Mbetax(3) =   Nbetax(1)
+c     Mbetax(i) + R(i,j) Nbetax(j)   [wrong]
+c     Mbetax(i) + R(j,i) Nbetax(j)   [correct]
+
+      end if
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
 c
 c compute flat-space components of Cx(*) parallel and perpendicular to vv(*)
 c
-	Ct    = t
-	Cx(1) = x
-	Cx(2) = y
-	Cx(3) = z
-
-			do 530 i=1,n
-			Cx_par_i  = 0.0d0
-			Cx_perp_i = 0.0d0
-				do 520 j=1,n
-				Cx_par_i  = Cx_par_i
-     $					    + parallel(i,j)*Cx(j)
-				Cx_perp_i = Cx_perp_i
-     $					    + perp    (i,j)*Cx(j)
-520				continue
-			Cx_par (i) = Cx_par_i
-			Cx_perp(i) = Cx_perp_i
-530			continue
+      Cx(0) = t
+      Cx(1) = x
+      Cx(2) = y
+      Cx(3) = z
+
+      do i=1,n
+         Cx_par_i  = 0.0
+         Cx_perp_i = 0.0
+         do j=1,n
+            Cx_par_i  = Cx_par_i  + parallel(i,j)*Cx(j)
+            Cx_perp_i = Cx_perp_i + perp    (i,j)*Cx(j)
+         end do
+         Cx_par (i) = Cx_par_i
+         Cx_perp(i) = Cx_perp_i
+      end do
 
 c
-c Lorentz-transform the Cactus coordinate to get the Model coordinates
+c Lorentz-transform and rotate the Cactus coordinate
+c to get the Model coordinates
 c
-	vdotCx = 0.0
-		do 600 i = 1,n
-		vdotCx = vdotCx + vv(i)*Cx(i)
-600		continue
 
-	Mt = gamma * (Ct - vdotCx)
-		do 610 i=1,n
-		Mx(i) = gamma * (Cx_par(i) - vv(i)*Ct)
-     $			+ Cx_perp(i)
-610		continue
+c Boost
+      vdotCx = 0.0
+      do i = 1,n
+         vdotCx = vdotCx + vv(i)*Cx(i)
+      end do
+
+      Mx(0) = gamma * (Cx(0) - vdotCx)
+      do i=1,n
+         Mx(i) = gamma * (Cx_par(i) - vv(i)*Cx(0)) + Cx_perp(i)
+      end do
+
+c Rotation
+      do i=0,n
+         Nx(i) = 0
+         do j = 0,n
+            Nx(i) = Nx(i) + R(i,j) * Mx(j)
+         end do
+      end do
+
 
 c
 c compute the Model 4-metric and inverse 4-metric at the Model coordinates
 c
-	call Exact__metric_for_model(
-     $			decoded_exact_model,
-     $			Mx(1), Mx(2), Mx(3), Mt,
-     $			Mgdd(0,0), Mgdd(0,1), Mgdd(0,2), Mgdd(0,3),
-     $			Mgdd(1,1), Mgdd(2,2), Mgdd(3,3),
-     $				Mgdd(1,2), Mgdd(2,3), Mgdd(1,3),
-     $			Mguu(0,0), Mguu(0,1), Mguu(0,2), Mguu(0,3),
-     $			Mguu(1,1), Mguu(2,2), Mguu(3,3),
-     $				Mguu(1,2), Mguu(2,3), Mguu(1,3),
-     $			psi, Tmunu_flag)
-
-	if (Tmunu_flag) then
-		write (warn_buffer, '(a,i8,a,a)')
-     $		      'exact_model = ', decoded_exact_model,
-     $		      'sets the stress-energy tensor',
-     $		      ' ==> we cannot Lorentz-boost it! :('
-		call CCTK_WARN(0, warn_buffer)
-	endif
+      call Exact__metric_for_model(
+     $     decoded_exact_model,
+     $     Nx(1), Nx(2), Nx(3), Nx(0),
+     $     Ngdd(0,0), Ngdd(0,1), Ngdd(0,2), Ngdd(0,3),
+     $     Ngdd(1,1), Ngdd(2,2), Ngdd(3,3),
+     $     Ngdd(1,2), Ngdd(2,3), Ngdd(1,3),
+     $     Nguu(0,0), Nguu(0,1), Nguu(0,2), Nguu(0,3),
+     $     Nguu(1,1), Nguu(2,2), Nguu(3,3),
+     $     Nguu(1,2), Nguu(2,3), Nguu(1,3),
+     $     psi, Tmunu_flag)
+
+      if (Tmunu_flag) then
+         write (warn_buffer, '(a,i8,a,a)')
+     $        'exact_model = ', decoded_exact_model,
+     $        'sets the stress-energy tensor',
+     $        ' ==> we cannot Lorentz-boost or rotate it!'
+         call CCTK_WARN(0, warn_buffer)
+      end if
 
 c
 c symmetrize the Model 4-metric and inverse 4-metric arrays
 c (the Exact__metric_for_model() call only set the upper triangles)
 c
-	Mgdd(1,0) = Mgdd(0,1)
-	Mgdd(2,0) = Mgdd(0,2)
-	Mgdd(2,1) = Mgdd(1,2)
-	Mgdd(3,0) = Mgdd(0,3)
-	Mgdd(3,1) = Mgdd(1,3)
-	Mgdd(3,2) = Mgdd(2,3)
-
-	Mguu(1,0) = Mguu(0,1)
-	Mguu(2,0) = Mguu(0,2)
-	Mguu(2,1) = Mguu(1,2)
-	Mguu(3,0) = Mguu(0,3)
-	Mguu(3,1) = Mguu(1,3)
-	Mguu(3,2) = Mguu(2,3)
+      Ngdd(1,0) = Ngdd(0,1)
+      Ngdd(2,0) = Ngdd(0,2)
+      Ngdd(2,1) = Ngdd(1,2)
+      Ngdd(3,0) = Ngdd(0,3)
+      Ngdd(3,1) = Ngdd(1,3)
+      Ngdd(3,2) = Ngdd(2,3)
+
+      Nguu(1,0) = Nguu(0,1)
+      Nguu(2,0) = Nguu(0,2)
+      Nguu(2,1) = Nguu(1,2)
+      Nguu(3,0) = Nguu(0,3)
+      Nguu(3,1) = Nguu(1,3)
+      Nguu(3,2) = Nguu(2,3)
 
 c
 c tensor-transorm (the upper triangle of) the 4-metric and inverse 4-metric
 c
-		do 730 Ca = 0,n
-		do 720 Cb = Ca,n
-		Cgdd_ab = 0.0d0
-			do 710 Ma = 0,n
-			do 700 Mb = 0,n
-			Cgdd_ab = Cgdd_ab
-     $				  +  Mgdd(Ma,Mb)
-     $				     * partial_Mx_wrt_Cx(Ma,Ca)
-     $				     * partial_Mx_wrt_Cx(Mb,Cb)
-700			continue
-710			continue
-		Cgdd(Ca,Cb) = Cgdd_ab
-720		continue
-730		continue
-
-		do 830 Ca = 0,n
-		do 820 Cb = Ca,n
-		Cguu_ab = 0.0d0
-			do 810 Ma = 0,n
-			do 800 Mb = 0,n
-			Cguu_ab = Cguu_ab
-     $				  +  Mguu(Ma,Mb)
-     $				     * partial_Cx_wrt_Mx(Ca,Ma)
-     $				     * partial_Cx_wrt_Mx(Cb,Mb)
-800			continue
-810			continue
-		Cguu(Ca,Cb) = Cguu_ab
-820		continue
-830		continue
+
+c Rotations
+      do i = 0,n
+         do j = 0,n
+            Mgdd(i,j) = 0
+            Mguu(i,j) = 0
+            do k = 0,n
+               do l = 0,n
+                  Mgdd(i,j) = Mgdd(i,j) + R(k,i) * R(l,j) * Ngdd(k,l)
+c     The inverse of R is also its transpose.  That means that the
+c     transpose of the inverse, which you would use for g^ij, is just R
+c     again.
+                  Mguu(i,j) = Mguu(i,j) + R(k,i) * R(l,j) * Nguu(k,l)
+               end do
+            end do
+         end do
+      end do
+
+c Boost
+      do Ca = 0,n
+         do Cb = Ca,n
+            Cgdd_ab = 0.0
+            do Ma = 0,n
+               do Mb = 0,n
+                  Cgdd_ab = Cgdd_ab
+     $                 +  Mgdd(Ma,Mb)
+     $                    * partial_Mx_wrt_Cx(Ma,Ca)
+     $                    * partial_Mx_wrt_Cx(Mb,Cb)
+               end do
+            end do
+            Cgdd(Ca,Cb) = Cgdd_ab
+         end do
+      end do
+
+      do Ca = 0,n
+         do Cb = Ca,n
+            Cguu_ab = 0.0
+            do Ma = 0,n
+               do Mb = 0,n
+                  Cguu_ab = Cguu_ab
+     $                 +  Mguu(Ma,Mb)
+     $                    * partial_Cx_wrt_Mx(Ca,Ma)
+     $                    * partial_Cx_wrt_Mx(Cb,Mb)
+               end do
+            end do
+            Cguu(Ca,Cb) = Cguu_ab
+         end do
+      end do
 
 c
 c unpack the Cactus-coordinates 4-metric and inverse 4-metric
 c into the corresponding output arguments
 c
-	gdtt = Cgdd(0,0)
-	gdtx = Cgdd(0,1)
-	gdty = Cgdd(0,2)
-	gdtz = Cgdd(0,3)
-	gdxx = Cgdd(1,1)
-	gdxy = Cgdd(1,2)
-	gdxz = Cgdd(1,3)
-	gdyy = Cgdd(2,2)
-	gdyz = Cgdd(2,3)
-	gdzz = Cgdd(3,3)
-
-	gutt = Cguu(0,0)
-	gutx = Cguu(0,1)
-	guty = Cguu(0,2)
-	gutz = Cguu(0,3)
-	guxx = Cguu(1,1)
-	guxy = Cguu(1,2)
-	guxz = Cguu(1,3)
-	guyy = Cguu(2,2)
-	guyz = Cguu(2,3)
-	guzz = Cguu(3,3)
-
-	return
-	end
+      gdtt = Cgdd(0,0)
+      gdtx = Cgdd(0,1)
+      gdty = Cgdd(0,2)
+      gdtz = Cgdd(0,3)
+      gdxx = Cgdd(1,1)
+      gdxy = Cgdd(1,2)
+      gdxz = Cgdd(1,3)
+      gdyy = Cgdd(2,2)
+      gdyz = Cgdd(2,3)
+      gdzz = Cgdd(3,3)
+
+      gutt = Cguu(0,0)
+      gutx = Cguu(0,1)
+      guty = Cguu(0,2)
+      gutz = Cguu(0,3)
+      guxx = Cguu(1,1)
+      guxy = Cguu(1,2)
+      guxz = Cguu(1,3)
+      guyy = Cguu(2,2)
+      guyz = Cguu(2,3)
+      guzz = Cguu(3,3)
+
+      end
diff --git a/src/decode_pars.F77 b/src/decode_pars.F77
index 1b7da46..d51bb4a 100644
--- a/src/decode_pars.F77
+++ b/src/decode_pars.F77
@@ -76,6 +76,8 @@ c black hole spacetimes
 	     decoded_exact_model = EXACT__Kerr_BoyerLindquist
 	elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild") .ne. 0) then
 	     decoded_exact_model = EXACT__Kerr_KerrSchild
+	elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild/spherical") .ne. 0) then
+	     decoded_exact_model = EXACT__Kerr_KerrSchild_spherical
 	elseif (CCTK_Equals(exact_model, "Schwarzschild-Lemaitre") .ne. 0) then
 	     decoded_exact_model = EXACT__Schwarzschild_Lemaitre
 	elseif (CCTK_Equals(exact_model, "multi-BH") .ne. 0) then
diff --git a/src/gauge.F77 b/src/gauge.F77
index 56088fb..b10bfa1 100644
--- a/src/gauge.F77
+++ b/src/gauge.F77
@@ -19,6 +19,7 @@ C $Header$
       integer i,j,k
       integer nx,ny,nz
 
+      CCTK_REAL tt, xx, yy, zz
       CCTK_REAL gxxjunk, gyyjunk, gzzjunk, 
      $     gxyjunk, gyzjunk, gxzjunk,
      $     hxxjunk, hyyjunk, hzzjunk, 
@@ -84,6 +85,11 @@ C
             do j=1,ny
                do i=1,nx
 
+                  tt = cctk_time
+                  xx = x(i,j,k) - cctk_time * shift_add_x
+                  yy = y(i,j,k) - cctk_time * shift_add_y
+                  zz = z(i,j,k) - cctk_time * shift_add_z
+
 C     Initialize the psi of exact
 C     (also to tell the models about the conformal_state)
                   if (conformal_state .ne. 0) then
@@ -103,7 +109,7 @@ C     (also to tell the models about the conformal_state)
 
                   call Exact__Bona_Masso_data(
      $                 decoded_exact_model,
-     $                 x(i,j,k), y(i,j,k), z(i,j,k), cctk_time,
+     $                 xx, yy, zz, tt,
      $                 gxxjunk, gyyjunk, gzzjunk, 
      $                 gxyjunk, gyzjunk, gxzjunk,
      $                 hxxjunk, hyyjunk, hzzjunk, 
@@ -147,6 +153,11 @@ C
             do j=1,ny
                do i=1,nx
 
+                  tt = cctk_time
+                  xx = x(i,j,k) - cctk_time * shift_add_x
+                  yy = y(i,j,k) - cctk_time * shift_add_y
+                  zz = z(i,j,k) - cctk_time * shift_add_z
+
 C     Initialize the psi of exact
 C     (also to tell the models about the conformal_state)
                   if (conformal_state .ne. 0) then
@@ -166,7 +177,7 @@ C     (also to tell the models about the conformal_state)
 
                   call Exact__Bona_Masso_data(
      $                 decoded_exact_model,
-     $                 x(i,j,k), y(i,j,k), z(i,j,k), cctk_time,
+     $                 xx, yy, zz, tt,
      $                 gxxjunk, gyyjunk, gzzjunk, 
      $                 gxyjunk, gyzjunk, gxzjunk,
      $                 hxxjunk, hyyjunk, hzzjunk, 
@@ -208,6 +219,11 @@ C
             do j=1,ny
                do i=1,nx
 
+                  tt = cctk_time
+                  xx = x(i,j,k) - cctk_time * shift_add_x
+                  yy = y(i,j,k) - cctk_time * shift_add_y
+                  zz = z(i,j,k) - cctk_time * shift_add_z
+
 C     Initialize the psi of exact
 C     (also to tell the models about the conformal_state)
                   if (conformal_state .ne. 0) then
@@ -227,7 +243,7 @@ C     (also to tell the models about the conformal_state)
 
                   call Exact__Bona_Masso_data(
      $                 decoded_exact_model,
-     $                 x(i,j,k), y(i,j,k), z(i,j,k), cctk_time,
+     $                 xx, yy, zz, tt,
      $                 gxxjunk, gyyjunk, gzzjunk, 
      $                 gxyjunk, gyzjunk, gxzjunk,
      $                 hxxjunk, hyyjunk, hzzjunk, 
diff --git a/src/include/Scalar_CalcTmunu.inc b/src/include/Scalar_CalcTmunu.inc
index 346fe09..e4ae8da 100644
--- a/src/include/Scalar_CalcTmunu.inc
+++ b/src/include/Scalar_CalcTmunu.inc
@@ -61,6 +61,8 @@ c        no stress-energy tensor in this model
 c        no stress-energy tensor in this model
        elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild) then
 c        no stress-energy tensor in this model
+       elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild_spherical) then
+c        no stress-energy tensor in this model
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
diff --git a/src/include/param_defs.inc b/src/include/param_defs.inc
index 10bbefa..d2374cb 100644
--- a/src/include/param_defs.inc
+++ b/src/include/param_defs.inc
@@ -38,16 +38,17 @@ c Minkowski spacetime
 #define EXACT__Minkowski_conf_wave		6
 
 c black hole spacetimes
-#define EXACT__Schwarzschild_EF		10
-#define EXACT__Schwarzschild_PG		11
-#define EXACT__Schwarzschild_BL		12
-#define EXACT__Schwarzschild_Novikov	13
-#define EXACT__Kerr_BoyerLindquist	14
-#define EXACT__Kerr_KerrSchild		15
-#define EXACT__Schwarzschild_Lemaitre	16
-#define EXACT__multi_BH			17
-#define EXACT__Alvi			18
-#define EXACT__Thorne_fakebinary	19
+#define EXACT__Schwarzschild_EF			10
+#define EXACT__Schwarzschild_PG			11
+#define EXACT__Schwarzschild_BL			12
+#define EXACT__Schwarzschild_Novikov		13
+#define EXACT__Kerr_BoyerLindquist		14
+#define EXACT__Kerr_KerrSchild			15
+#define EXACT__Kerr_KerrSchild_spherical	16
+#define EXACT__Schwarzschild_Lemaitre		17
+#define EXACT__multi_BH				18
+#define EXACT__Alvi				19
+#define EXACT__Thorne_fakebinary		20
 
 c cosmological spacetimes
 #define EXACT__Lemaitre			50
diff --git a/src/metric.F77 b/src/metric.F77
index a93ea12..f28c948 100644
--- a/src/metric.F77
+++ b/src/metric.F77
@@ -154,6 +154,15 @@ c
      $        guxx, guyy, guzz, guxy, guyz, guxz,
      $        psi, Tmunu_flag)
 
+      elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild_spherical) then
+         call Exact__Kerr_KerrSchild_spherical(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz,
+     $        psi, Tmunu_flag)
+
       elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Lemaitre) then
          call Exact__Schwarzschild_Lemaitre(
      $        x, y, z, t,
diff --git a/src/metrics/Kerr_KerrSchild.F77 b/src/metrics/Kerr_KerrSchild.F77
index f2fb45a..11b6ce1 100644
--- a/src/metrics/Kerr_KerrSchild.F77
+++ b/src/metrics/Kerr_KerrSchild.F77
@@ -38,6 +38,7 @@ c output arguments
 
 c local variables
       CCTK_REAL boostv, eps, m, a
+      integer   power
 
 c local variables
       CCTK_REAL gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
@@ -54,6 +55,7 @@ C     to the J/m definition used in the code here.
 
       boostv = Kerr_KerrSchild__boost_v
       eps    = Kerr_KerrSchild__epsilon
+      power  = Kerr_KerrSchild__power
       m      = Kerr_KerrSchild__mass
       a      = m*Kerr_KerrSchild__spin
 
@@ -79,7 +81,18 @@ C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
 
       r02 = 0.5d0 * (rho02 - a**2) 
      $     + sqrt(0.25d0 * (rho02 - a**2)**2 + a**2 * z0**2)
-      r0 = (r02**2 + eps**4)**0.25d0
+      r0 = sqrt(r02)
+      if (Kerr_KerrSchild__parabolic .eq. 0) then
+C        Use a power law to avoid the singularity
+         r0 = (r0**power + eps**power)**(1.0d0/power)
+      else
+         if (r0 .lt. eps) then
+            r0 = (eps + r0**2 / eps) / 2
+         end if
+      end if
+C     Another idea:
+C        r0 = r0 + eps * exp(-x/eps)
+
       costheta0 = z0 / r0
       
 C     Coefficient H. Note this transforms as a scalar, so does not carry
diff --git a/src/metrics/Kerr_KerrSchild_spherical.F77 b/src/metrics/Kerr_KerrSchild_spherical.F77
new file mode 100644
index 0000000..417b224
--- /dev/null
+++ b/src/metrics/Kerr_KerrSchild_spherical.F77
@@ -0,0 +1,182 @@
+C     Kerr-Schild form of boosted rotating black hole.
+C     Program g_ab = eta_ab + H l_a l_b, g^ab = eta^ab - H l^a l^b.
+C     Here eta_ab is Minkowski in Cartesian coordinates, H is a scalar,
+C     and l is a null vector.
+C     
+C     The coordinates are distorted, such that the event horizon is
+C     a coordinate sphere.
+C
+C     Author: Erik Schnetter <schnetter@cct.lsu.edu>
+C     This formulation was invented by Nils Dorband <dorband@cct.lsu.edu>
+C
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+      subroutine Exact__Kerr_KerrSchild_spherical (
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz, 
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $     gutt, gutx, guty, gutz, 
+     $     guxx, guyy, guzz, guxy, guyz, guzx,
+     $     psi, Tmunu_flag)
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+c input arguments
+      CCTK_REAL x, y, z, t
+
+c output arguments
+      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
+     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $       gutt, gutx, guty, gutz, 
+     $       guxx, guyy, guzz, guxy, guyz, guzx
+      CCTK_REAL psi
+      LOGICAL   Tmunu_flag
+
+c local variables
+      CCTK_REAL boostv, eps, m, a
+
+c local variables
+      CCTK_REAL t1, x1, y1, z1, rho1,
+     $     gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
+     $     lt0, lx0, ly0, lz0, hh, lt, lx, ly, lz
+
+      CCTK_REAL gd(3,3), gdt(3,3), det, jac(3,3)
+
+      integer i, j, k, l
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+C This is a vacuum spacetime with no cosmological constant
+      Tmunu_flag = .false.
+
+C     Get parameters of the exact solution,
+C     and convert from parameter file spin parameter J/m^2
+C     to the J/m definition used in the code here.
+
+      boostv = Kerr_KerrSchild__boost_v
+      eps    = Kerr_KerrSchild__epsilon
+      m      = Kerr_KerrSchild__mass
+      a      = m*Kerr_KerrSchild__spin
+
+C     Distort the coordinates such that the event horizon is a
+C     coordinate sphere
+      rho1 = sqrt (x**2 + y**2 + z**2)
+      rho1 = (rho1**4 + eps**4) ** 0.25d0
+      t1 = t
+      x1 = x - a * y / rho1
+      y1 = y + a * x / rho1
+      z1 = z
+
+C     Boost factor.
+
+      gamma = 1 / sqrt(1 - boostv**2)
+
+C     Lorentz transform t,x,y,z -> t0,x0,y0,z0. 
+C     t0 is never used, but is here for illustration, and we introduce
+C     x0 and y0 also only for clarity.
+C     Note that z0 = 0 means z = vt for the BH.
+
+      t0 = gamma * ((t1 - Kerr_KerrSchild__t) - boostv * (z1 - Kerr_KerrSchild__z))
+      z0 = gamma * ((z1 - Kerr_KerrSchild__z) - boostv * (t1 - Kerr_KerrSchild__t))
+      x0 = x1 - Kerr_KerrSchild__x
+      y0 = y1 - Kerr_KerrSchild__y
+
+C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
+
+      r0 = rho1
+      costheta0 = z0 / r0
+      
+C     Coefficient H.  Note this transforms as a scalar, so does not carry
+C     the suffix 0.
+      hh = m * r0 / (r0**2 + a**2 * costheta0**2)
+      
+C     Components of l_a in rest frame. Note indices down.
+      lt0 = 1.d0
+      lx0 = (r0 * x0 + a * y0) / (r0**2 + a**2)
+      ly0 = (r0 * y0 - a * x0) / (r0**2 + a**2)
+      lz0 = z0 / r0
+
+C     Now boost it to coordinates x, y, z, t.
+C     This is the reverse Lorentz transformation, but applied
+C     to a one-form, so the sign of boostv is the same as the forward
+C     Lorentz transformation applied to the coordinates.
+
+      lt = gamma * (lt0 - boostv * lz0) 
+      lz = gamma * (lz0 - boostv * lt0) 
+      lx = lx0
+      ly = ly0
+
+C     Down metric.  g_ab = flat_ab + H l_a l_b
+
+      gdtt = - 1.d0 + 2.d0 * hh * lt * lt
+      gdtx = 2.d0 * hh * lt * lx
+      gdty = 2.d0 * hh * lt * ly
+      gdtz = 2.d0 * hh * lt * lz
+
+      gd(1,1) = 1.d0 + 2.d0 * hh * lx * lx
+      gd(2,2) = 1.d0 + 2.d0 * hh * ly * ly
+      gd(3,3) = 1.d0 + 2.d0 * hh * lz * lz
+      gd(1,2) = 2.d0 * hh * lx * ly
+      gd(2,3) = 2.d0 * hh * ly * lz
+      gd(3,1) = 2.d0 * hh * lz * lx
+      gd(2,1) = gd(1,2)
+      gd(3,2) = gd(2,3)
+      gd(1,3) = gd(3,1)
+
+C     Transform the tensor basis back
+      jac(1,1) = 1 + (a*x*y) / (rho1**3)
+      jac(1,2) = - a * (x**2+z**2) / (rho1**3)
+      jac(1,3) = a*y*z / (rho1**3)
+      
+      jac(2,1) = a * (y**2+z**2) / (rho1**3)
+      jac(2,2) = 1 - a*x*y / (rho1**3)
+      jac(2,3) = - a*x*z / (rho1**3)
+      
+      jac(3,1) = 0
+      jac(3,2) = 0
+      jac(3,3) = 1
+
+      do i = 1, 3
+         do j = 1, 3
+            gdt(i,j) = 0
+            do k = 1, 3
+               do l = 1, 3
+                  gdt(i,j) = gdt(i,j) + gd(k,l) * jac(k,i) * jac(l,j)
+               end do
+            end do
+         end do
+      end do
+
+      gdxx = gdt(1,1)
+      gdyy = gdt(2,2)
+      gdzz = gdt(3,3)
+      gdxy = gdt(1,2)
+      gdyz = gdt(2,3)
+      gdzx = gdt(3,1)
+
+C     Up metric.  g^ab = flat^ab - H l^a l^b.
+C     Notice that g^ab = g_ab and l^i = l_i and l^0 = - l_0 in flat spacetime.
+      gutt = - 1.d0 - 2.d0 * hh * lt * lt 
+      gutx = 2.d0 * hh * lt * lx
+      guty = 2.d0 * hh * lt * ly
+      gutz = 2.d0 * hh * lt * lz
+
+      det = gdxx*gdyy*gdzz + 2*gdxy*gdzx*gdyz
+     .    - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2
+
+      guxx = (gdyy*gdzz - gdyz**2)/det
+      guyy = (gdxx*gdzz - gdzx**2)/det
+      guzz = (gdxx*gdyy - gdxy**2)/det
+
+      guxy = (gdzx*gdyz - gdzz*gdxy)/det
+      guyz = (gdxy*gdzx - gdxx*gdyz)/det
+      guzx = (gdxy*gdyz - gdyy*gdzx)/det
+
+      end
diff --git a/src/metrics/make.code.defn b/src/metrics/make.code.defn
index 12226e2..1e90de0 100644
--- a/src/metrics/make.code.defn
+++ b/src/metrics/make.code.defn
@@ -19,6 +19,7 @@ SRCS = Minkowski.F77			\
        Schwarzschild_Novikov.F77	\
        Kerr_BoyerLindquist.F77		\
        Kerr_KerrSchild.F77		\
+       Kerr_KerrSchild_spherical.F77	\
        Schwarzschild_Lemaitre.F77	\
        multi_BH.F77			\
        Thorne_fakebinary.F77		\