Adding code to integrate proper distance between horizons. At the moment,

it only finds the distance betwen horizons 2 and 3 if they exist.  If they
don't, I am not sure what it will do yet.  Still testing it.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@246 89daf98e-ef62-4674-b946-b8ff9de2216c

4 files changed, 349 insertions, 2 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index 5880c6b..0033fd0 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -1403,11 +1403,12 @@
          call AHFinder_find3(CCTK_ARGUMENTS,mtype)
       end if
 
+
 !     ****************************
 !     ***   WRITE DATA FILES   ***
 !     ****************************
 
-      call AHFinder_Output (CCTK_ARGUMENTS, report, status, horizon, mtype, intarea_h)
+      call AHFinder_Output(CCTK_ARGUMENTS,report,status,horizon,mtype,intarea_h)
 
 
 !     ****************
@@ -1583,6 +1584,15 @@
       end do
 
 
+!     *******************************************************
+!     ***   CALCULATE PROPER DISTANCES BETWEEN HORIZONS   ***
+!     *******************************************************
+
+      if (find3) then
+         call AHFinder_dis(CCTK_ARGUMENTS)
+      end if
+
+
 !     **************************************************
 !     ***   TRANSFORM PHYSICAL TO CONFORMAL METRIC   ***
 !     **************************************************
diff --git a/src/AHFinder_dis.F b/src/AHFinder_dis.F
new file mode 100644
index 0000000..271d697
--- /dev/null
+++ b/src/AHFinder_dis.F
@@ -0,0 +1,336 @@
+/*@@
+  @file      AHFinder_dis.F
+  @date      September 2001
+  @author    Miguel Alcubierre
+  @desc 
+             Routine to calculate proper distances between horizons.
+  @enddesc 
+@@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine AHFinder_dis(CCTK_ARGUMENTS)
+
+      use AHFinder_dat
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      integer i,l,m
+      integer np,npoints
+      integer interp_handle,coord_system_handle
+      integer num_arrays,ierror
+
+      integer, dimension(6) :: in_array_indices
+
+      CCTK_REAL LEGEN
+      CCTK_REAL zero,half,one,two,pi
+      CCTK_REAL x12,y12,z12,r12,rho12
+      CCTK_REAL theta,phi,cost,sint,cosp,sinp
+      CCTK_REAL x1,y1,z1,r1
+      CCTK_REAL x2,y2,z2,r2
+      CCTK_REAL ddx,ddy,ddz,d12
+      CCTK_REAL aux,cosa,sina
+
+      CCTK_REAL, allocatable, dimension(:) :: xa,ya,za
+      CCTK_REAL, allocatable, dimension(:) :: txx,tyy,tzz,txy,txz,tyz
+
+      character*200 disf
+
+
+!     **************************
+!     ***   DEFINE NUMBERS   ***
+!     **************************
+
+      zero  = 0.0D0
+      half  = 0.5D0
+      one   = 1.0D0
+      two   = 2.0D0
+
+      pi = acos(-one)
+
+
+!     *****************************************
+!     ***   FIND LINE CONNECTING HORIZONS   ***
+!     *****************************************
+
+!     For the moment, horizons 1 and 2 only.
+
+      x12 = ahf_xc_2 - ahf_xc_1
+      y12 = ahf_yc_2 - ahf_yc_1
+      z12 = ahf_zc_2 - ahf_zc_1
+
+      rho12 = sqrt(x12**2 + y12**2)
+      r12   = sqrt(x12**2 + y12**2 + z12**2)
+
+
+!     ***********************************************************
+!     ***   FIND LOCATION OF HORIZONS ALONG CONNECTING LINE   ***
+!     ***********************************************************
+
+!     Find angles for horizon 1.
+
+      if (r12.ne.zero) then
+         cost = z12/r12
+      else
+         return
+      end if
+
+      sint = sqrt(one - cost**2)
+
+      if (y12.ne.zero) then
+         phi = atan2(y12,x12)
+      else
+         phi = zero
+      end if
+
+      sinp = sin(phi)
+      cosp = cos(phi)
+
+!     Find point in horizon 1.
+
+      r1 = c0(0)
+
+      do l=1+stepz,lmax,1+stepz
+         r1 = r1 + c0(l)*LEGEN(l,0,cost)
+      end do
+
+      if (nonaxi) then
+         do m=1,lmax
+            aux = dble(m)*phi
+            sina = sin(aux)
+            cosa = cos(aux)
+            do l=m,lmax
+               aux = LEGEN(l,m,cost)
+               r1 = r1 + aux*cc(l,m)*cosa
+               if (.not.refy) then
+                  r1 = r1 + aux*cs(l,m)*sina
+               end if
+            end do
+         end do
+      end if
+
+      x1 = ahf_xc_1 + r1*sint*cosp
+      y1 = ahf_yc_1 +r1*sint*sinp
+      z1 = ahf_zc_1 +r1*cost
+
+!     Find angles for horizon 2.
+
+      cost = -cost
+      sint = +sint
+
+      cosp = -cosp
+      sinp = -sinp
+
+      phi = phi + pi
+
+!     Find point in horizon 2.
+
+      r2 = c0(0)
+
+      do l=1+stepz,lmax,1+stepz
+         r2 = r2 + c0(l)*LEGEN(l,0,cost)
+      end do
+
+      if (nonaxi) then
+         do m=1,lmax
+            aux = dble(m)*phi
+            sina = sin(aux)
+            cosa = cos(aux)
+            do l=m,lmax
+               aux = LEGEN(l,m,cost)
+               r2 = r2 + aux*cc(l,m)*cosa
+               if (.not.refy) then
+                  r2 = r2 + aux*cs(l,m)*sina
+               end if
+            end do
+         end do
+      end if
+
+      x2 = ahf_xc_2 + r2*sint*cosp
+      y2 = ahf_yc_2 + r2*sint*sinp
+      z2 = ahf_zc_2 + r2*cost
+
+
+!     **************************************
+!     ***   ALLOCATE MEMORY FOR ARRAYS   ***
+!     **************************************
+
+      aux = dsqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)/dx
+
+      np = int(aux)
+
+      if (myproc.eq.0) then
+
+         allocate(xa(0:np))
+         allocate(ya(0:np))
+         allocate(za(0:np))
+
+         allocate(txx(0:np))
+         allocate(tyy(0:np))
+         allocate(tzz(0:np))
+         allocate(txy(0:np))
+         allocate(txz(0:np))
+         allocate(tyz(0:np))
+
+      else
+
+         allocate(xa(1))
+         allocate(ya(1))
+         allocate(za(1))
+
+         allocate(txx(1))
+         allocate(tyy(1))
+         allocate(tzz(1))
+         allocate(txy(1))
+         allocate(txz(1))
+         allocate(tyz(1))
+
+      end if
+
+
+!     *************************************
+!     ***   FIND INTERPOLATING POINTS   ***
+!     *************************************
+
+      if (myproc.eq.0) then
+
+         npoints = np + 1
+
+         do i=0,np
+
+            aux = dble(i)/dble(np)
+
+            xa(i) = x1 + aux*(x2-x1)
+            ya(i) = y1 + aux*(y2-y1)
+            za(i) = z1 + aux*(z2-z1)
+
+         end do
+
+      else
+
+         npoints = 1
+
+         xa = half*(ahf_xc_1 + ahf_xc_2)
+         ya = half*(ahf_xc_1 + ahf_xc_2)
+         za = half*(ahf_xc_1 + ahf_xc_2)
+
+      end if
+
+
+!     ******************************
+!     ***   INTERPOLATE METRIC   ***
+!     ******************************
+
+!     Interpolator and coordinate system handle.
+
+      interp_handle = -1
+      coord_system_handle = -1
+
+      if (interpolation_order.eq.1) then
+        call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
+      else if (interpolation_order .eq. 2) then
+        call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
+      else if (interpolation_order .eq. 3) then
+        call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
+      endif
+
+      call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
+
+      if (interp_handle .lt. 0 .or. coord_system_handle .lt. 0) then
+        call CCTK_WARN (0, "Couldn't get handles for interpolation operator and/or coordinate system")
+      endif
+
+!     Interpolation.
+
+      call CCTK_VarIndex (in_array_indices(1), "einstein::gxx")
+      call CCTK_VarIndex (in_array_indices(2), "einstein::gyy")
+      call CCTK_VarIndex (in_array_indices(3), "einstein::gzz")
+      call CCTK_VarIndex (in_array_indices(4), "einstein::gxy")
+      call CCTK_VarIndex (in_array_indices(5), "einstein::gxz")
+      call CCTK_VarIndex (in_array_indices(6), "einstein::gyz")
+
+      num_arrays = 6
+      call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
+     .                    npoints, num_arrays, num_arrays,
+     .                    xa, ya, za,
+     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+     .                    CCTK_VARIABLE_REAL,
+     .                    in_array_indices(1), in_array_indices(2),
+     .                    in_array_indices(3), in_array_indices(4),
+     .                    in_array_indices(5), in_array_indices(6),
+     .                    txx, tyy, tzz, txy, txz, tyz,
+     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+
+
+!     ******************************
+!     ***   INTEGRATE DISTANCE   ***
+!     ******************************
+
+      ddx = (x2-x1)/dble(np)
+      ddy = (y2-y1)/dble(np)
+      ddz = (z2-z1)/dble(np)
+
+      if (myproc.eq.0) then
+
+!        Integrate using trapezoidal rule.
+
+         d12 = zero
+
+         do i=1,np
+
+            d12 = d12 + half
+     .          *(dsqrt(txx(i-1)*ddx**2  + tyy(i-1)*ddy**2  + tzz(i-1)*ddz**2
+     .          +  two*(txy(i-1)*ddx*ddy + txz(i-1)*ddx*ddz + tyz(i-1)*ddy*ddz))
+     .          + dsqrt(txx(i  )*ddx**2  + tyy(i  )*ddy**2  + tzz(i  )*ddz**2
+     .          +  two*(txy(i  )*ddx*ddy + txz(i  )*ddx*ddz + tyz(i  )*ddy*ddz)))
+
+         end do
+
+      end if
+
+
+!     ******************
+!     ***   OUTPUT   ***
+!     ******************
+
+      disf = filestr(1:nfile) // "/ahf_d12" // ".tl"
+
+      if (myproc.eq.0) then
+
+        if (ahf_ncall.eq.1) then
+           open(1,file=disf,form='formatted',status='replace')
+           write(1,"(A14)") '" Horizon distance'
+        else
+           open(1,file=disf,form='formatted',status='old',position='append')
+        end if
+
+        write(1,"(2ES14.6)") cctk_time,d12
+
+        close(1)
+
+      end if
+
+
+!     *****************************
+!     ***   DEALLOCATE MEMORY   ***
+!     *****************************
+
+      deallocate(xa,ya,za)
+      deallocate(txx,tyy,tzz,txy,txz,tyz)
+
+
+!     ***************
+!     ***   END   ***
+!     ***************
+
+      end subroutine AHFinder_dis
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
index 9f38e65..76d20eb 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # Source files in this directory
 
 SRCS += AHFinder_dat.F AHFinder.F \
-        AHFinder_fun.F AHFinder_exp.F AHFinder_int.F \
+        AHFinder_fun.F AHFinder_exp.F AHFinder_dis.F AHFinder_int.F \
         AHFinder_min.F AHFinder_pow.F AHFinder_leg.F \
         AHFinder_flow.F AHFinder_mask.F AHFinder_output.F \
         AHFinder_find3.F AHFinder_calcsigma.F \
diff --git a/src/make.code.deps b/src/make.code.deps
index 62384ac..0cf58e5 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -4,6 +4,7 @@
 
 AHFinder.F.o:            AHFinder_dat.F.o
 AHFinder_calcsigma.F.o:  AHFinder_dat.F.o
+AHFinder_dis.F.o:        AHFinder_dat.F.o
 AHFinder_exp.F.o:        AHFinder_dat.F.o
 AHFinder_fun.F.o:        AHFinder_dat.F.o
 AHFinder_gau.F.o:        AHFinder_dat.F.o