************************************************************************
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

*  Computes the lapse and shift vector
************************************************************************
      subroutine setgauge(nx,ny,nz,mask,x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      implicit none

      integer  nx,ny,nz
      real*8  mask(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz)
      real*8   bha(2), bhr(2,3)
      real*8   alpha(nx,ny,nz),
     &         beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

c      call setgauge_ef(nx,ny,nz,mask,
c     &           x,y,z,bha,bhr,
c     &           alpha,beta1,beta2,beta3)

c      call setgauge_iso(nx,ny,nz,mask,
c     &           x,y,z,bha,bhr,
c     &           alpha,beta1,beta2,beta3)

c      call setgauge_alg(nx,ny,nz,mask,
c     &           x,y,z,bha,bhr,
c     &           alpha,beta1,beta2,beta3)

      call setgauge_geo(nx,ny,nz,mask,
     &           x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      return
      end
************************************************************************
      subroutine setgauge_iso(nx,ny,nz,mask,
     &           x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      implicit none

      integer  nx,ny,nz
      real*8  mask(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz)
      real*8   bha(2), bhr(2,3)
      real*8   alpha(nx,ny,nz),
     &         beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

c  local variables

      integer  i,j,k
      real*8   r1,r2,alpha1,alpha2

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx

         r1 = sqrt( (x(i)-bhr(1,1))**2
     &             +(y(j)-bhr(1,2))**2
     &             +(z(k)-bhr(1,3))**2 )
         r2 = sqrt( (x(i)-bhr(2,1))**2
     &             +(y(j)-bhr(2,2))**2
     &             +(z(k)-bhr(2,3))**2 )

            alpha1 = 1.0e0
         if ( bha(1) .ne. 0.0e0 ) then
            alpha1 = (r1-bha(1))/(r1+bha(1))
         end if
    
            alpha2 = 1.0e0
         if ( bha(2) .ne. 0.0e0 ) then
            alpha2 = (r2-bha(2))/(r2+bha(2))
         end if

         alpha(i,j,k) = mask(i,j,k)*alpha1*alpha2
         beta1(i,j,k) = 0.0e0
         beta2(i,j,k) = 0.0e0
         beta3(i,j,k) = 0.0e0

      end do
      end do
      end do

      return
      end
************************************************************************
      subroutine setgauge_ef(nx,ny,nz,mask,x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      implicit none

      integer  nx,ny,nz
      real*8  mask(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz)
      real*8   bha(2), bhr(2,3)
      real*8   alpha(nx,ny,nz),
     &         beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

c  local variables

      integer  i,j,k
      real*8   rr

      do k = 1, nz
        do j = 1, ny
          do i = 1, nx

         rr = sqrt( (x(i)-bhr(1,1))**2
     &             +(y(j)-bhr(1,2))**2
     &             +(z(k)-bhr(1,3))**2 )

			if(mask(i,j,k) .ne. 0.0)then
				alpha(i,j,k) = 1.0e0  / sqrt(1.0e0 + 2.0e0 / rr)

				beta1(i,j,k) = -2.0e0 * x(i) / (rr * (2.0e0 + rr))
            beta2(i,j,k) = -2.0e0 * y(j) / (rr * (2.0e0 + rr))
            beta3(i,j,k) = -2.0e0 * z(k) / (rr * (2.0e0 + rr))
			else
				alpha(i,j,k) = 0.0e0
				beta1(i,j,k) = 0.0e0
				beta2(i,j,k) = 0.0e0
				beta3(i,j,k) = 0.0e0
			end if


          end do
        end do
      end do

      return
      end
************************************************************************
*  algebraic slicings
      subroutine setgauge_alg(nx,ny,nz,mask,
     &           x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      implicit none

      integer  nx,ny,nz
      real*8  mask(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz)
      real*8   bha(2), bhr(2,3)
      real*8   alpha(nx,ny,nz),
     &         beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

c  local variables

      integer  i,j,k
      real*8   rr

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx

         rr = sqrt( (x(i)-bhr(1,1))**2
     &             +(y(j)-bhr(1,2))**2
     &             +(z(k)-bhr(1,3))**2 )


			if(mask(i,j,k) .ne. 0.0)then
         	alpha(i,j,k) =1.0 + 6.0* log(1.0 + bha(1)/rr)
			else
         	alpha(i,j,k) = 0.0
			end if
         beta1(i,j,k) = 0.0e0
         beta2(i,j,k) = 0.0e0
         beta3(i,j,k) = 0.0e0

      end do
      end do
      end do

      return
      end
************************************************************************
*  geodesic slicings
      subroutine setgauge_geo(nx,ny,nz,mask,
     &           x,y,z,bha,bhr,
     &           alpha,beta1,beta2,beta3)

      implicit none

      integer  nx,ny,nz
      real*8  mask(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz)
      real*8   bha(2), bhr(2,3)
      real*8   alpha(nx,ny,nz),
     &         beta1(nx,ny,nz), beta2(nx,ny,nz), beta3(nx,ny,nz)

c  local variables

      integer  i,j,k
      real*8   rr

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx

         rr = sqrt( (x(i)-bhr(1,1))**2
     &             +(y(j)-bhr(1,2))**2
     &             +(z(k)-bhr(1,3))**2 )


			if(mask(i,j,k) .ne. 0.0)then
         	alpha(i,j,k) =1.0 
			else
         	alpha(i,j,k) = 0.0
			end if
         beta1(i,j,k) = 0.0e0
         beta2(i,j,k) = 0.0e0
         beta3(i,j,k) = 0.0e0

      end do
      end do
      end do

      return
      end
************************************************************************