#include "cctk.h"
#include "cctk_Parameters.h"

      subroutine  boostrotmetric(
     $     x, y, z, t,
     $     gdtt, gdtx, gdty, gdtz, 
     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
     $     gutt, gutx, guty, gutz, 
     $     guxx, guyy, guzz, guxy, guyz, guzx)

      implicit none

      DECLARE_CCTK_PARAMETERS

      CCTK_REAL x, y, z, t,
     $    gdtt, gdtx, gdty, gdtz, 
     $    gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
     $    gutt, gutx, guty, gutz, 
     $    guxx, guyy, guzz, guxy, guyz, guzx
      CCTK_REAL a, b, mu0, mu1, lam1, mu2, lam2,
     $    lam3, mu4, lam4, mu5, lam5, num, div, f,
     $    elam, emu0, delta, gfunc, tmp
      CCTK_REAL h, d, numlim

      external gfunc

      logical firstcall

      data firstcall /.true./
      save firstcall, h, d, numlim

C     Get parameters of the exact solution.

      if (firstcall) then

         h = boostrotscale
         d = boostrotstrength

         numlim = boostrotsafedistance

         firstcall = .false.

      end if

C     Intermediate quantities.

      a = x**2 + y**2
      b = z**2 - t**2

      num = (0.5d0 * (a + b) - h)**2 + 2.d0 * h * a

C     Make sure we are not sitting on one of the two source wordlines,
C     given by x = y = 0, z = +/- sqrt(h^2 + t^2)

      if (num / h**4 .le. numlim) stop 'too close to source wordline'

      div = 1.d0 / sqrt(num**3)
      f = d**2  * ((0.25d0 * (a + b)**2 - h**2)**2 
     $     - 0.5d0 * h**2 * a * b) / num**4

      mu0 = - d * div * (0.5d0 * a**2 + h * a)
      mu1 = - d * div * (0.5d0 * b + a - h)
      lam1 = d * div * (0.5d0 * b - h) - a * f
      mu2 = gfunc(b, mu1)
      lam2 =  gfunc(b, lam1)

      lam3 = d * div * (0.5d0 * b**2 - h * b)
      lam4 = - d * div * (0.5d0 * a + h) - b * f
      mu4 = - d * div * (0.5d0 * a + b + h) 
      mu5 = gfunc(a, - mu4)
      lam5 = gfunc(a, lam4)

      elam = exp(lam3 + a * lam4)
      emu0 = exp(mu0)
      delta = exp(lam3) * (mu5 - lam5)

C     All nonvanishing metric coefficients (downstairs).

      gdxx = elam + y**2 * Delta
      gdyy = elam + x**2 * Delta
      gdxy = - x * y * Delta
      gdzz = emu0 * (1.d0 + lam2 * z**2 - mu2 * t**2)
      gdtz = - emu0 * z * t * (lam2 - mu2)
      gdtt = - emu0 * (1.d0 + mu2 * z**2 - lam2 * t**2)

C     Others.

      gdzx = 0.d0
      gdyz = 0.d0
      gdtx = 0.d0
      gdty = 0.d0

C     It is clear that the 3-metric is always spacelike in the xy plane. So
C     we only need to check that gdzz is positive.

      if (gdzz .le. 0.d0) then
         write(*,*) 'WARNING 3-metric not spacelike in boostrot at'
         write(*,*) 't =', t, 'z =', z
         write(*,*) 'x =', x, 'y =', y
         pause
      end if

C     Calculate inverse metric. That is not too difficult as it is
c     in block-diagonal form.

      tmp = gdtt * gdzz - gdtz**2

      if (tmp .eq. 0.d0) then
         write(*,*) 'boostrot metric inverse failed in tz plane'
         STOP
      end if

      gutt = gdzz / tmp
      guzz = gdtt / tmp
      gutz = - gdtz / tmp
      
      tmp = gdxx * gdyy - gdxy**2

      if (tmp .eq. 0.d0) then
         write(*,*) 'boostrot metric inverse failed in xy plane'
         STOP
      end if

      guxx = gdyy / tmp
      guyy = gdxx / tmp
      guxy = - gdxy / tmp

      guzx = 0.d0
      guyz = 0.d0
      gutx = 0.d0
      guty = 0.d0

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

C     Calculates g = [exp (x f) - 1] / x as a power series for small x, 
C     so that the expression is regular at x = 0.

      function gfunc(x, f)
    
      implicit none

      integer n

      CCTK_REAL x, f, gfunc
      CCTK_REAL sum, tmp

      if (abs(x*f) .ge. 1.d-6) then
         gfunc = (exp(x*f) - 1.d0) / x
      else
         sum = 0.d0
         tmp = f
         do n=1,10
            tmp = tmp / dble(n) 
            sum = sum + tmp
            tmp = tmp * x * f
         end do
         gfunc = sum
      end if

      return
      end