c#######################################################################
c#
c# Set initial data for two black hole collision in axisymmetry.
c# Sets everything except Psi,Lapse,ShiftE,ShiftX
c#
c#######################################################################
      subroutine cadez_id(pzp,pzm,szp,szm,isom,
     .                    KK,RHOzc,Zzc,CJACzc,ext,fil)
      implicit none
c-----------------------------------------------------------------------
      integer ext(4),fil(4)
      real*8 pzp,pzm,szp,szm,isom
      real*8 KK(ext(1):ext(3),ext(2):ext(4))
      real*8 RHOzc(ext(1):ext(3),ext(2):ext(4))
      real*8 Zzc(ext(1):ext(3),ext(2):ext(4))
      real*8 CJACzc(ext(1):ext(3),ext(2):ext(4))
c-----------------------------------------------------------------------
      integer i,j
      real*8 crho,cz
      real*8 axx,ayy,azz,axy,axz,ayz
      real*8 ZERO,TWO
      parameter(ZERO = 0.D0,TWO = 2.0D0)
c-----------------------------------------------------------------------
c
      do j=fil(2),fil(4)
        do i=fil(1),fil(3)
          crho = RHOzc(i,j)
          cz = Zzc(i,j)
c
          call excrvimg(crho,cz,
     .                  axx,ayy,azz,axy,axz,ayz,
     .                  pzp,pzm,szp,szm,isom,10)
          KK(i,j) = axx*axx + ayy*ayy + azz*azz +
     .                  TWO*(axy*axy + axz*axz + ayz*ayz)
c
        end do
      end do
c
      return
      end
c#######################################################################
c#
c# Compute the INVERSION SYMMETRIC cartesian components of the
c# extrinsic curvature for two black holes with axisymmetric linear
c# and angular momenta.
c# [compute on phi=0 plane (y=0 in cartesian coordinates)]
c#
c#######################################################################
      subroutine excrvimg(rho,z,axx,ayy,azz,axy,axz,ayz,
     .                    pzp,pzm,szp,szm,isom,nterms)
      implicit none
c-----------------------------------------------------------------------
      real*8 rho,z,axx,ayy,azz,axy,axz,ayz,pzp,pzm,szp,szm,isom
      integer nterms
c-----------------------------------------------------------------------
      integer n
      real*8 rhop,rhom,zp,zm,fp,fm
      real*8 mp(3,3),mm(3,3),acurv(6),atmp(6),atmm(6)
      real*8 ZERO,HALF,ONE
      integer Ixx,Iyy,Izz,Ixy,Ixz,Iyz
      parameter(ZERO = 0.D0,HALF = 0.5D0,ONE = 1.D0,
     .          Ixx = 1,Iyy = 2,Izz = 3,Ixy = 4,Ixz = 5,Iyz = 6)
c-----------------------------------------------------------------------
c- Common block from "chicord.F" subroutines for Cadez Coordinates.
c- Values for alpha,zetap,zetam are needed and are set by calling
c- "cadef".
c-
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c- Initialize quantities for Inversion Iteration.
c-----------------------------------------------------------------------
        rhop = rho
        rhom = rho
        zp = z
        zm = z
        fp = ONE
        fm = ONE
        mp(1,1) = ONE
        mp(1,2) = ZERO
        mp(1,3) = ZERO
        mp(2,1) = ZERO
        mp(2,2) = ONE
        mp(2,3) = ZERO
        mp(3,1) = ZERO
        mp(3,2) = ZERO
        mp(3,3) = ONE
        mm(1,1) = ONE
        mm(1,2) = ZERO
        mm(1,3) = ZERO
        mm(2,1) = ZERO
        mm(2,2) = ONE
        mm(2,3) = ZERO
        mm(3,1) = ZERO
        mm(3,2) = ZERO
        mm(3,3) = ONE
c
c-----------------------------------------------------------------------
c- Compute the "bare" extrinsic curvature.
c-----------------------------------------------------------------------
      call excrv(rho,z,acurv,pzp,pzm,szp,szm,zetap,zetam)
c
c-----------------------------------------------------------------------
c- Add NTERMS inversion terms.
c-----------------------------------------------------------------------
      do n = 1,nterms
        if(mod(n,2).eq.1) then
c n is odd
          call update(fp,mp,rhop,zp,ONE,zetap)
          call update(fm,mm,rhom,zm,alpha,zetam)
          call excrv(rhop,zp,atmp,pzp,pzm,szp,szm,zetap,zetam)
          call excrv(rhom,zm,atmm,pzp,pzm,szp,szm,zetap,zetam)
          call imagcomb(fp,mp,atmp)
          call imagcomb(fm,mm,atmm)
          acurv(Ixx) = acurv(Ixx) + isom*(atmp(Ixx) + atmm(Ixx))
          acurv(Ixy) = acurv(Ixy) + isom*(atmp(Ixy) + atmm(Ixy))
          acurv(Ixz) = acurv(Ixz) + isom*(atmp(Ixz) + atmm(Ixz))
          acurv(Iyy) = acurv(Iyy) + isom*(atmp(Iyy) + atmm(Iyy))
          acurv(Iyz) = acurv(Iyz) + isom*(atmp(Iyz) + atmm(Iyz))
          acurv(Izz) = acurv(Izz) + isom*(atmp(Izz) + atmm(Izz))
        else
c n is even
          call update(fp,mp,rhop,zp,alpha,zetam)
          call update(fm,mm,rhom,zm,ONE,zetap)
          call excrv(rhop,zp,atmp,pzp,pzm,szp,szm,zetap,zetam)
          call excrv(rhom,zm,atmm,pzp,pzm,szp,szm,zetap,zetam)
          call imagcomb(fp,mp,atmp)
          call imagcomb(fm,mm,atmm)
          acurv(Ixx) = acurv(Ixx) + atmp(Ixx) + atmm(Ixx)
          acurv(Ixy) = acurv(Ixy) + atmp(Ixy) + atmm(Ixy)
          acurv(Ixz) = acurv(Ixz) + atmp(Ixz) + atmm(Ixz)
          acurv(Iyy) = acurv(Iyy) + atmp(Iyy) + atmm(Iyy)
          acurv(Iyz) = acurv(Iyz) + atmp(Iyz) + atmm(Iyz)
          acurv(Izz) = acurv(Izz) + atmp(Izz) + atmm(Izz)
        end if
      end do
      axx = acurv(Ixx)
      ayy = acurv(Iyy)
      azz = acurv(Izz)
      axy = acurv(Ixy)
      axz = acurv(Ixz)
      ayz = acurv(Iyz)
c
      return
      end
c#######################################################################
c#
c# Compute the cartesian components of the base extrinsic curvature
c# for two black holes with axisymmetric linear and angular momenta.
c# [compute on phi=0 plane (y=0 in cartesian coordinates)]
c#
c#######################################################################
      subroutine excrv(rho,z,acurv,pzp,pzm,szp,szm,zetap,zetam)
      implicit none
c-----------------------------------------------------------------------
      real*8 rho,z,acurv(6),pzp,pzm,szp,szm,zetap,zetam
c-----------------------------------------------------------------------
      real*8 rpm1,rmm1,rpm2,rmm2
      real*8 nxp,nyp,nzp,nxm,nym,nzm
      real*8 hpdotnp,hpdotnm,pdotn
      real*8 fxp,fyp,fzp,fxm,fym,fzm
      real*8 ZERO,HALF,ONE,THREEHALF,TWO
      integer Ixx,Iyy,Izz,Ixy,Ixz,Iyz
      parameter(ZERO = 0.D0,HALF = 0.5D0,ONE = 1.D0,
     .          THREEHALF = 1.5D0,TWO = 2.D0,
     .          Ixx = 1,Iyy = 2,Izz = 3,Ixy = 4,Ixz = 5,Iyz = 6)
c-----------------------------------------------------------------------
c
      rpm1 = ONE/sqrt(rho*rho + (z - zetap)**2)
      rmm1 = ONE/sqrt(rho*rho + (z - zetam)**2)
      rpm2 = THREEHALF*(rpm1**2)
      rmm2 = THREEHALF*(rmm1**2)
      nxm = rho*rmm1
      nxp = rho*rpm1
      nzp = (z - zetap)*rpm1
      nzm = (z - zetam)*rmm1
      hpdotnp = pzp*nzp
      hpdotnm = pzm*nzm
      pdotn = -rpm2*hpdotnp - rmm2*hpdotnm
      hpdotnp = HALF*hpdotnp
      hpdotnm = HALF*hpdotnm
      fxp = nxp*hpdotnp*rpm2
      fxm = nxm*hpdotnm*rmm2
      fyp = nxp*szp*rpm1*rpm2
      fym = nxm*szm*rmm1*rmm2
      fzp = (pzp + nzp*hpdotnp)*rpm2
      fzm = (pzm + nzm*hpdotnm)*rmm2
      acurv(Ixx) = TWO*(nxp*fxp + nxm*fxm) + pdotn
      acurv(Iyy) = pdotn
      acurv(Izz) = TWO*(nzp*fzp + nzm*fzm) + pdotn
      acurv(Ixy) = nxp*fyp + nxm*fym
      acurv(Ixz) = nzp*fxp + nzm*fxm + nxp*fzp + nxm*fzm
      acurv(Iyz) = nzp*fyp + nzm*fym
c
      return
      end
c#######################################################################
c#
c# Update the Inversion Iteration quantities assuming axisymmetry.
c# [compute on phi=0 plane (y=0 in cartesian coordinates)]
c#
c#######################################################################
      subroutine update(f,m,rho,z,a,zeta)
      implicit none
c-----------------------------------------------------------------------
      real*8 f,m(3,3),rho,z,a,zeta
c-----------------------------------------------------------------------
      integer k
      real*8 rm2,fxx,fxz,fzz,tmpx
      real*8 ONE,TWO
      parameter(ONE = 1.D0,TWO = 2.D0)
c-----------------------------------------------------------------------
c
      rm2 = ONE/(rho*rho + (z - zeta)**2)
      fxx = ONE - TWO*rm2*(rho**2)
      fxz = z - zeta
      fzz = ONE - TWO*rm2*(fxz**2)
      fxz = -TWO*rho*rm2*fxz
      do 10 k=1,3
        tmpx = m(k,1)
        m(k,1) = tmpx*fxx + m(k,3)*fxz
        m(k,3) = tmpx*fxz + m(k,3)*fzz
10    continue
      rm2 = rm2/(a**2)
      f = f*(rm2**3)
      rho = rm2*rho
      z = rm2*(z - zeta) + zeta
c
      return
      end
c#######################################################################
c#
c# Combine the inversion matrix and scale factor to give an image term.
c#
c#######################################################################
      subroutine imagcomb(f,m,acurv)
      implicit none
c-----------------------------------------------------------------------
      real*8 f,m(3,3),acurv(6)
c-----------------------------------------------------------------------
      real*8 axx,ayy,azz,axy,axz,ayz
      real*8 TWO
      integer Ixx,Iyy,Izz,Ixy,Ixz,Iyz
      parameter(TWO=2.D0,
     .          Ixx = 1,Iyy = 2,Izz = 3,Ixy = 4,Ixz = 5,Iyz = 6)
c-----------------------------------------------------------------------
c
      axx = m(1,1)*m(1,1)*acurv(Ixx) 
     .       + m(1,3)*(m(1,3)*acurv(Izz) + TWO*m(1,1)*acurv(Ixz))
      ayy = m(2,1)*(m(2,1)*acurv(Ixx) + TWO*acurv(Ixy))
     .       + acurv(Iyy) + TWO*m(2,3)*acurv(Iyz)
     .       + m(2,3)*(m(2,3)*acurv(Izz) + TWO*m(2,1)*acurv(Ixz))
      azz = m(3,1)*m(3,1)*acurv(Ixx)
     .       + m(3,3)*(m(3,3)*acurv(Izz) + TWO*m(3,1)*acurv(Ixz))
      axy = m(1,1)*(m(2,1)*acurv(Ixx) +
     .                acurv(Ixy) + m(2,3)*acurv(Ixz)) +
     .      m(1,3)*(m(2,1)*acurv(Ixz) + acurv(Iyz) + m(2,3)*acurv(Izz))
      axz = m(1,1)*(m(3,1)*acurv(Ixx) + m(3,3)*acurv(Ixz)) +
     .         m(1,3)*(m(3,1)*acurv(Ixz) + m(3,3)*acurv(Izz))
      ayz = m(2,1)*(m(3,1)*acurv(Ixx) + m(3,3)*acurv(Ixz)) +
     .         m(3,1)*acurv(Ixy) + m(3,3)*acurv(Iyz) +
     .         m(2,3)*(m(3,1)*acurv(Ixz) + m(3,3)*acurv(Izz))
      acurv(Ixx) = f*axx
      acurv(Ixy) = f*axy
      acurv(Ixz) = f*axz
      acurv(Iyy) = f*ayy
      acurv(Iyz) = f*ayz
      acurv(Izz) = f*azz
c
      return
      end
