Added support for finding surface area. However only for one surface and

in full mode. It is parallelized though. Need to put in a routine to
detect how many surface are present, their symmetries and handle various
symmetries.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@79 2a26948c-0e4f-0410-aee8-f1d3e353619c

7 files changed, 493 insertions, 62 deletions

diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90
index 0a72481..3c60ec8 100644
--- a/src/EHFinder_Constants.F90
+++ b/src/EHFinder_Constants.F90
@@ -12,6 +12,9 @@ module EHFinder_Constants
   CCTK_REAL, parameter :: eight = 8.0d0
   CCTK_REAL, parameter :: ten = 10.0d0
   CCTK_REAL, parameter :: half = 0.5d0
+  CCTK_REAL, parameter :: onethird = one / three
+  CCTK_REAL, parameter :: twothirds = two / three
+  CCTK_REAL, parameter :: fourthirds = four / three
   CCTK_REAL, parameter :: quarter = 0.25d0
   CCTK_REAL, parameter :: tenth = 0.1d0
   CCTK_REAL, parameter :: huge = 1d23
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index e51542d..4a62ae0 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -1,4 +1,4 @@
-! Finding the surface(s). Right now only on one processor and in full mode.
+! Finding the surface(s). Right now only in full mode.
 ! $Header$
 
 #include "cctk.h"
@@ -15,18 +15,19 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k
+  CCTK_INT :: i, j, k, im, jm
   CCTK_REAL, parameter :: eps = 1.0d-10
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
+  CCTK_INT :: interp_handle2
 
   CCTK_INT, dimension(4) :: bbox
-  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl
 
-  CCTK_REAL, dimension(3) :: center_loc, center
+  CCTK_REAL, dimension(3) :: center_loc
   CCTK_INT :: nc_loc, nc
   CCTK_REAL :: initial_radius, ncinv, maxdr_loc, maxdr, maxf_loc, maxf
 
-  CCTK_REAL :: dtheta, dphi, delta, rloc
+  CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv, delta, rloc
   CCTK_POINTER, dimension(3) :: interp_coords
   CCTK_POINTER, dimension(4) :: out_array
   CCTK_POINTER :: in_array
@@ -69,41 +70,52 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
   end if
   
-  nc_loc = 0; center_loc = zero
-  delta = minval ( cctk_delta_space )
-  do k = kzl, kzr
-    do j = jyl, jyr
-      do i = ixl, ixr
-        if ( eh_mask(i,j,k) .eq. 0 ) then
-          if ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
-               ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
-               ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
-               ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
-               ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
-               ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) then
-            nc_loc = nc_loc + 1
-            center_loc(1) = center_loc(1) + x(i,j,k)
-            center_loc(2) = center_loc(2) + y(i,j,k)
-            center_loc(3) = center_loc(3) + z(i,j,k)
+  if ( use_user_center ) then
+    center(1) = center_x
+    center(2) = center_y
+    center(3) = center_z
+  else
+    nc_loc = 0; center_loc = zero
+    delta = minval ( cctk_delta_space )
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+          if ( eh_mask(i,j,k) .eq. 0 ) then
+            if ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
+                 ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
+                 ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
+                 ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
+                 ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
+                 ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) then
+              nc_loc = nc_loc + 1
+              center_loc(1) = center_loc(1) + x(i,j,k)
+              center_loc(2) = center_loc(2) + y(i,j,k)
+              center_loc(3) = center_loc(3) + z(i,j,k)
+            end if
           end if
-        end if
+        end do
       end do
     end do
-  end do
 
-  call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, &
+    call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, &
                                     nc_loc, nc, CCTK_VARIABLE_INT )
-  call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, &
-                                      center_loc, center, 3, &
-                                      CCTK_VARIABLE_REAL )
-  ncinv = one / nc
-  center = ncinv * center
-  if ( symx ) center(1) = zero
-  if ( symy ) center(2) = zero
-  if ( symz ) center(3) = zero
+    call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, &
+                                        center_loc, center, 3, &
+                                        CCTK_VARIABLE_REAL )
+    ncinv = one / nc
+    center = ncinv * center
+!    print*,center
+    if ( symx ) center(1) = zero
+    if ( symy ) center(2) = zero
+    if ( symz ) center(3) = zero
+  end if
+!  print*,'Center : ',center
+!  print*,'use_user_center : ',use_user_center
 
   dtheta = pi / ( ntheta - 1 )
   dphi = two * pi / ( nphi - 1 )
+  dthetainv = one / dtheta
+  dphiinv = one / dphi
   do i = 1, lsh(1)
     ctheta(i,:) = dtheta * ( i + lbnd(1) - 1 )
   end do
@@ -111,25 +123,31 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     cphi(:,j) = dphi * ( j +lbnd(2) - 1 )
   end do
 
-  do j = 1, lsh(2)
-    do i = 1, lsh(1)
-      sinthetacosphi(i,j) = sin(ctheta(i,j)) * cos(cphi(i,j))
-      sinthetasinphi(i,j) = sin(ctheta(i,j)) * sin(cphi(i,j))
-      costheta(i,j) = cos(ctheta(i,j))
-    end do
-  end do
+  sintheta = sin(ctheta)
+  costheta = cos(ctheta)
+  sinphi = sin(cphi)
+  cosphi = cos(cphi)
+!  do j = 1, lsh(2)
+!    do i = 1, lsh(1)
+!      sinthetacosphi(i,j) = sin(ctheta(i,j)) * cos(cphi(i,j))
+!      sinthetasinphi(i,j) = sin(ctheta(i,j)) * sin(cphi(i,j))
+!      costheta(i,j) = cos(ctheta(i,j))
+!    end do
+!  end do
   
   drsurf = two * delta
   rsurf = zero
   n_since_last_reduction = -1
 
-  call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
+!  call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
+  call CCTK_InterpHandle ( interp_handle, "Hermite polynomial interpolation" )
   if ( interp_handle .lt. 0 ) then
     call CCTK_WARN( 0, "Cannot get handle for interpolation. &
  &  Forgot to activate an implementation providing interpolation operators??" )
   end if
 
-  call Util_TableCreateFromString ( table_handle, "order=3" )
+!  call Util_TableCreateFromString ( table_handle, "order=3" )
+  call Util_TableCreateFromString ( table_handle, "order=2" )
   if ( table_handle .lt. 0 ) then
     call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
   end if
@@ -155,8 +173,8 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     do j = 1, lsh(2)
       do i = 1, lsh(1)
         rloc = rsurf(i,j) + drsurf(i,j)
-        interp_x(i,j) = center(1) + rloc * sinthetacosphi(i,j)
-        interp_y(i,j) = center(2) + rloc * sinthetasinphi(i,j)
+        interp_x(i,j) = center(1) + rloc * sintheta(i,j) * cosphi(i,j)
+        interp_y(i,j) = center(2) + rloc * sintheta(i,j) * sinphi(i,j)
         interp_z(i,j) = center(3) + rloc * costheta(i,j)
         if ( n_since_last_reduction(i,j) .ge. 0 ) then 
           n_since_last_reduction(i,j) = n_since_last_reduction(i,j) + 1
@@ -173,7 +191,15 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     maxdr_loc = maxval ( drsurf )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
                                     maxdr_loc, maxdr, CCTK_VARIABLE_REAL )
+
+    maxf_loc = maxval ( abs ( f_interp ) )
+    call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
+                                    maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
+    maxl = maxloc ( abs ( f_interp ) )
+!    print*,maxl,rsurf(maxl(1),maxl(2)),&
+!           drsurf(maxl(1),maxl(2)),&
+!           f_interp(maxl(1),maxl(2))
     if ( maxdr .lt. delta ) exit find_starting_points
 
     do j = 1, lsh(2)
@@ -191,7 +217,9 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     end do
     
   end do find_starting_points
- 
+!  print*,'Max values: ',maxval(rsurf),maxval(f_interp)
+!  print*,'Min values: ',minval(rsurf),minval(f_interp)
+!  print*,'found starting points' 
   out_array(2) = CCTK_PointerTo(dfdx_interp)
   out_array(3) = CCTK_PointerTo(dfdy_interp)
   out_array(4) = CCTK_PointerTo(dfdz_interp)
@@ -209,12 +237,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
 ! find the surface points in all directions simultaneously.
 
-  find_surface_points: do
+  find_surface_points: do k = 1, 15
 
     do j = 1, lsh(2)
       do i = 1, lsh(1)
-        interp_x(i,j) = center(1) + rsurf(i,j) * sinthetacosphi(i,j)
-        interp_y(i,j) = center(2) + rsurf(i,j) * sinthetasinphi(i,j)
+        interp_x(i,j) = center(1) + rsurf(i,j) * sintheta(i,j) * cosphi(i,j)
+        interp_y(i,j) = center(2) + rsurf(i,j) * sintheta(i,j) * sinphi(i,j)
         interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j)
       end do
     end do
@@ -229,12 +257,21 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
                                     maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
+    maxl = maxloc ( abs ( f_interp ) )
+!    print*,k,maxf
+!    print*,maxf,maxl,f_interp(maxl(1),maxl(2)),rsurf(maxl(1),maxl(2)), &
+!                     dfdx_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* &
+!                                                  cosphi(maxl(1),maxl(2))+ &
+!                     dfdy_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* &
+!                                                  sinphi(maxl(1),maxl(2))+ &
+!                     dfdz_interp(maxl(1),maxl(2))*costheta(maxl(1),maxl(2))
+!    print*
     if ( maxf .gt. eps ) then
       do j = 1, lsh(2)
         do i = 1, lsh(1)
           rsurf(i,j) = rsurf(i,j) - f_interp(i,j) / &
-                       ( dfdx_interp(i,j) * sinthetacosphi(i,j) + &
-                         dfdy_interp(i,j) * sinthetasinphi(i,j) + &
+                       ( dfdx_interp(i,j) * sintheta(i,j) * cosphi(i,j) + &
+                         dfdy_interp(i,j) * sintheta(i,j) * sinphi(i,j) + &
                          dfdz_interp(i,j) * costheta(i,j) )
         end do
       end do
@@ -244,7 +281,110 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 !    call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
 
   end do find_surface_points
+  if ( k .gt. 15 ) call CCTK_WARN(1,'Newton iteration did not converge! Area may be inaccurate')
+!  print*,'found_surface_points'
+
+  drdtheta = zero
+  drdphi = zero
+! Calculate the derivatives of r.
+
+! First for interior points
+  do j = 2, lsh(2) - 1
+    do i = 2, lsh(1) - 1
+      drdtheta(i,j) = half * dthetainv * ( rsurf(i+1,j) - rsurf(i-1,j) )
+      drdphi(i,j) = half * dphiinv * ( rsurf(i,j+1) - rsurf(i,j-1) )
+    end do
+  end do
+
+! If boundary 1 is physical update those points.
+  if ( bbox(1) .eq. 1 ) then
+    do j = 2, lsh(2) - 1
+      drdtheta(1,j) = half * dthetainv * ( - three * rsurf(1,j) + &
+                                             four * rsurf(2,j) - &
+                                             rsurf(3,j) )
+      drdphi(1,j) = half * dphiinv * ( rsurf(1,j+1) - rsurf(1,j-1) )
+    end do
+  end if
+
+! If boundary 2 is physical update those points.
+  if ( bbox(2) .eq. 1 ) then
+    im = lsh(1)
+    do j = 2, lsh(2) - 1
+      drdtheta(im,j) = half * dthetainv * ( three * rsurf(im,j) - &
+                                            four * rsurf(im-1,j) + &
+                                            rsurf(im-2,j) )
+      drdphi(im,j) = half * dphiinv * ( rsurf(im,j+1) - rsurf(im,j-1) )
+    end do
+  end if
+
+! If boundary 3 is physical update those points.
+  if ( bbox(3) .eq. 1 ) then
+    do i = 2, lsh(1) - 1
+      drdtheta(i,1) = half * dthetainv * ( rsurf(i+1,1) - rsurf(i-1,1) )
+      drdphi(i,1) = half * dphiinv * ( - three * rsurf(i,1) + &
+                                         four * rsurf(i,2) - &
+                                         rsurf(i,3) )
+    end do
+  end if
+
+! If boundary 4 is physical update those points.
+  if ( bbox(4) .eq. 1 ) then
+    jm = lsh(2)
+    do i = 2, lsh(1) - 1
+      drdtheta(i,jm) = half * dthetainv * ( rsurf(i+1,jm) - rsurf(i-1,jm) )
+      drdphi(i,jm) = half * dphiinv * ( three * rsurf(i,jm) - &
+                                        four * rsurf(i,jm-1) + &
+                                        rsurf(i,jm-2) )
+    end do
+  end if
+
+! If corner 1 is physical update that point.
+  if ( bbox(1) .eq. 1 .and. bbox(3) .eq. 1 ) then
+    drdtheta(1,1) = half * dthetainv * ( - three * rsurf(1,1) + &
+                                             four * rsurf(2,1) - &
+                                             rsurf(3,1) )
+    drdphi(1,1) = half * dphiinv * ( - three * rsurf(1,1) + &
+                                       four * rsurf(1,2) - &
+                                       rsurf(1,3) )
+  end if
+
+! If corner 2 is physical update that point.
+  if ( bbox(1) .eq. 1 .and. bbox(4) .eq. 1 ) then
+    jm = lsh(2)
+    drdtheta(1,jm) = half * dthetainv * ( - three * rsurf(1,jm) + &
+                                             four * rsurf(2,jm) - &
+                                             rsurf(3,jm) )
+    drdphi(1,jm) = half * dphiinv * ( three * rsurf(1,jm) - &
+                                      four * rsurf(1,jm-1) + &
+                                      rsurf(1,jm-2) )
+  end if
+ 
+! If corner 3 is physical update that point.
+  if ( bbox(2) .eq. 1 .and. bbox(3) .eq. 1 ) then
+    im = lsh(1)
+    drdtheta(im,1) = half * dthetainv * ( three * rsurf(im,1) - &
+                                          four * rsurf(im-1,1) + &
+                                          rsurf(im-2,1) )
+    drdphi(im,1) = half * dphiinv * ( - three * rsurf(im,1) + &
+                                        four * rsurf(im,2) - &
+                                        rsurf(im,3) )
+  end if
+ 
+! If corner 4 is physical update that point.
+  if ( bbox(2) .eq. 1 .and. bbox(4) .eq. 1 ) then
+    im = lsh(1)
+    jm = lsh(2)
+    drdtheta(im,jm) = half * dthetainv * ( three * rsurf(im,jm) - &
+                                           four * rsurf(im-1,jm) + &
+                                           rsurf(im-2,jm) )
+    drdphi(im,jm) = half * dphiinv * ( three * rsurf(im,jm) - &
+                                       four * rsurf(im,jm-1) + &
+                                       rsurf(im,jm-2) )
+  end if
+
+  call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
 
+!  print*,'Calculated derivatives'
 !  print*,rsurf
   call Util_TableDestroy ( status, table_handle )
   
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 95d0e14..4aa43ff 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -95,5 +95,25 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
     call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
   end if
 
+
+  call CCTK_CoordRegisterSystem ( ierr, 2, 'cart2d' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN(0,'Could not register a 2D coordinate system as "cart2d"')
+  end if
+  call CCTK_CoordRegisterData ( ierr, 1, 'ehfinder::ctheta', &
+                                                 'x', 'cart2d' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN(0,'Could not register ctheta as the 1st coordinate for "cart2d"')
+  end if
+  call CCTK_CoordRegisterData ( ierr, 2, 'ehfinder::cphi', &
+                                                 'y', 'cart2d' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN(0,'Could not register cphi as the 2nd coordinate for "cart2d"')
+  end if
+
+  call CCTK_INFO('2d coordinate system registered')
+
+! set up weights for Simpsons rule integration
+
   return
 end subroutine EHFinder_Init
diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90
index c46a43a..04a2a46 100644
--- a/src/EHFinder_Integrate.F90
+++ b/src/EHFinder_Integrate.F90
@@ -28,7 +28,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
   CCTK_REAL :: dfsquare, dfsquare_i
   CCTK_REAL, dimension(3) :: vec1, vec2
   CCTK_REAL :: tmp1, tmp2, tmp3
-  CCTK_REAL :: idetg, sqrtdg2
+  CCTK_REAL :: detg, idetg, sqrtdg2
   CCTK_REAL :: a, b, c, r1, r2, r3, q
 
   interface
@@ -63,8 +63,8 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
             tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
   
-            idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + &
-                                              gxz(i,j,k)*tmp3 )
+            detg = gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3
+            idetg = one / detg
 
 ! Invert the spatial 3-metric
 
@@ -84,8 +84,10 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
                               ginv(2,2) * dfy(i,j,k)**2 + &
                               two * ginv(2,3) * dfy(i,j,k)*dfz(i,j,k) + &
                               ginv(3,3) *  dfz(i,j,k)**2 )
+!            print*,dfsquare
 
             dfsquare_i = one / dfsquare
+!            print*,dfsquare_i
 !            dfsquare_i = one / sqrt ( dfx(i,j,k)**2 + &
 !                                      dfy(i,j,k)**2 + &
 !                                      dfz(i,j,k)**2 )
@@ -93,6 +95,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             norm(1) = dfx(i,j,k) * dfsquare_i
             norm(2) = dfy(i,j,k) * dfsquare_i
             norm(3) = dfz(i,j,k) * dfsquare_i
+!            print*,norm
 
 !            norm(1) = dfx(i,j,k)
 !            norm(2) = dfy(i,j,k)
@@ -131,6 +134,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             g2yy = gyy(i,j,k) - norm(2)**2
             g2yz = gyz(i,j,k) - norm(2) * norm(3)
             g2zz = gzz(i,j,k) - norm(3)**2
+!            print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
 
 !            a = - ( g2xx + g2yy + g2zz )
             b = g2xx * g2yy + g2xx * g2zz + g2yy * g2zz - &
@@ -139,11 +143,14 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
 !                  g2xx * g2yy * g2zz - two * g2xy * g2xz * g2yz
 
 !            q = - half * ( a + sign ( one, a ) * sqrt ( a**2 - 4 * b ) )
-!            print*,q,b/q,sqrt(b)
-            sqrtdg2 = sqrt ( sqrt (b) )
-!            pause
+!            print*,'factors : ',a,b,c,q,b/q,sqrt(b)
+!            print*
+!            sqrtdg2 = sqrt ( sqrt (b) )
+            sqrtdg2 = sqrt (b)
 !            call cubic ( a, b, c, r1, r2, r3 )
-!           print*
+!            print*,'roots : ',r1,r2,r3
+!            print*
+!            pause
 !           print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
 !           print*
 
@@ -189,6 +196,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
           tmp1 = g2yy*g2zz - g2yz**2
           tmp2 = g2xz*g2yz - g2xy*g2zz
           tmp3 = g2xy*g2yz - g2xz*g2yy        
+!          print*,tmp1,tmp2,tmp3
 
 !          print*,tmp1,tmp2,tmp3
 !
@@ -198,23 +206,25 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
 !            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
 !            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
 
-!           print*,sqrtdg2 !,g22xx,g22xy,g22yy
+!           print*,sqrtdg2,g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3
 !           pause
 
             local_area = local_area + sqrtdg2 * dirac(f(i,j,k),eps,eps_inv) * &
-                         sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
+                                      dfsquare * sqrt(detg)
+!                         sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
           end if
         end do
       end do
     end do
   end if
   local_area = dv * local_area
+!  print*,'local_area = ',local_area,sqrtdg2
 
   call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
                                     local_area, area, CCTK_VARIABLE_REAL )
 
   eh_area = area
-!  print*,'Area = ',eh_area
+  print*,'Area = ',eh_area
 
   return
 end subroutine EHFinder_Integrate
@@ -312,12 +322,19 @@ subroutine cubic ( a, b, c, x1, x2, x3 )
 
   CCTK_REAL, intent(in) :: a, b, c
   CCTK_REAL, intent(out) :: x1, x2, x3
-  CCTK_REAL :: q, r
+  CCTK_REAL :: q, r, theta
+  CCTK_REAL, parameter :: third = 1.0d0 / 3.0d0, &
+                          pi2 = 6.2831853071795864769d0
 
   q = ( a**2 - 3.0d0 * b ) / 9.0d0
   r = ( 2.0d0 * a**3 - 9.0d0 * a * b + 27.0d0 * c ) / 54.0d0
+  theta = acos ( r / sqrt( q**3 ) )
+
+  x1 = - 2.0d0 * sqrt( q ) * cos( third * theta ) - third * a
+  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta + pi2 ) ) - third * a 
+  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta - pi2 ) ) - third * a 
 
   print*,a,b,c
-  print*,q,r
-  pause
+  print*,q,r,q**3-r**2
+!  pause
 end subroutine cubic
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
new file mode 100644
index 0000000..24ad3e9
--- /dev/null
+++ b/src/EHFinder_Integrate2.F90
@@ -0,0 +1,249 @@
+! Integrating over the surface(s). Right now only in full mode.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Integrate2(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k, im, jm
+  CCTK_INT :: interp_handle, table_handle, coord_system_handle, int_var
+
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+
+  CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv
+  CCTK_REAL :: dxdth, dxdph, dydth, dydph, dzdth, dzdph
+!  CCTK_REAL :: gtt, gtp, gpp
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(7) :: in_array, out_array
+  character(len=30) :: info_message
+  CCTK_INT, dimension(7), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL /)
+  
+  call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get bounding box for surface arrays" )
+  end if
+!  print*,bbox
+  call CCTK_GroupgshGN ( ierr, cctkGH, 2, gsh, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get global size for surface arrays" )
+  end if
+!  print*,gsh
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" )
+  end if
+!  print*,lbnd
+  call CCTK_GroupubndGN ( ierr, cctkGH, 2, ubnd, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get upper bounds for surface arrays" )
+  end if
+!  print*,ubnd
+  call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+!  print*,lsh
+  call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+  
+  dtheta = pi / ( ntheta - 1 )
+  dphi = two * pi / ( nphi - 1 )
+  dthetainv = one / dtheta
+  dphiinv = one / dphi
+  
+  call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
+  if ( interp_handle .lt. 0 ) then
+    call CCTK_WARN( 0, "Cannot get handle for interpolation. &
+ &  Forgot to activate an implementation providing interpolation operators??" )
+  end if
+
+  call Util_TableCreateFromString ( table_handle, "order=3" )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+  end if
+  
+  call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+  if ( coord_system_handle .lt. 0) then
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. &
+ &  Forgot to activate an implementation providing coordinates ??" )
+  endif
+
+  interp_coords(1) = CCTK_PointerTo(interp_x)
+  interp_coords(2) = CCTK_PointerTo(interp_y)
+  interp_coords(3) = CCTK_PointerTo(interp_z)
+
+  out_array(1) = CCTK_PointerTo(gxxi)
+  out_array(2) = CCTK_PointerTo(gxyi)
+  out_array(3) = CCTK_PointerTo(gxzi)
+  out_array(4) = CCTK_PointerTo(gyyi)
+  out_array(5) = CCTK_PointerTo(gyzi)
+  out_array(6) = CCTK_PointerTo(gzzi)
+  out_array(7) = CCTK_PointerTo(psii)
+
+! find the cartesian coordinates for the interpolation points
+
+!  print*,center
+  do j = 1, lsh(2)
+    do i = 1, lsh(1)
+      interp_x(i,j) = center(1) + rsurf(i,j) * sintheta(i,j) * cosphi(i,j)
+      interp_y(i,j) = center(2) + rsurf(i,j) * sintheta(i,j) * sinphi(i,j)
+      interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j)
+    end do
+  end do
+
+  call CCTK_VarIndex ( in_array(1), "admbase::gxx" )
+  call CCTK_VarIndex ( in_array(2), "admbase::gxy" )
+  call CCTK_VarIndex ( in_array(3), "admbase::gxz" )
+  call CCTK_VarIndex ( in_array(4), "admbase::gyy" )
+  call CCTK_VarIndex ( in_array(5), "admbase::gyz" )
+  call CCTK_VarIndex ( in_array(6), "admbase::gzz" )
+
+  if ( CCTK_EQUALS ( metric_type, "static conformal" ) ) then
+    call CCTK_VarIndex ( in_array(7), "staticconformal::psi" )
+
+    call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 7, in_array, &
+                                 7, out_types, out_array )
+    gxxi = psii**4 * gxxi
+    gxyi = psii**4 * gxyi
+    gxzi = psii**4 * gxzi
+    gyyi = psii**4 * gyyi
+    gyzi = psii**4 * gyzi
+    gzzi = psii**4 * gzzi
+  else
+    call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 6, in_array(1:6), &
+                                 6, out_types(1:6), out_array(1:6) )
+  end if
+
+  do j = 1, lsh(2)
+    do i = 1, lsh(1)
+!      print*,i,j,ctheta(i,j),cphi(i,j)
+      dxdth = ( drdtheta(i,j) * sintheta(i,j) + &
+                rsurf(i,j) * costheta(i,j) ) * cosphi(i,j)
+!      print*,dxdth,rsurf(i,j)*costheta(i,j)*cosphi(i,j)
+      dydth = ( drdtheta(i,j) * sintheta(i,j) + &
+                rsurf(i,j) * costheta(i,j) ) * sinphi(i,j)
+!      print*,dydth,rsurf(i,j)*costheta(i,j)*sinphi(i,j)
+      dzdth = drdtheta(i,j) * costheta(i,j) - rsurf(i,j) * sintheta(i,j)
+!      print*,dzdth,-rsurf(i,j)*sintheta(i,j)
+      dxdph = ( drdphi(i,j) * cosphi(i,j) - &
+                rsurf(i,j) * sinphi(i,j) ) * sintheta(i,j)
+!      print*,dxdph,-rsurf(i,j)*sintheta(i,j)*sinphi(i,j)
+      dydph = ( drdphi(i,j) * sinphi(i,j) + &
+                rsurf(i,j) * cosphi(i,j) ) * sintheta(i,j)
+!      print*,dydph,rsurf(i,j)*sintheta(i,j)*cosphi(i,j)
+      dzdph = drdphi(i,j) * costheta(i,j)
+!      print*,dzdph,0
+!      print*
+
+      gtt(i,j) = dxdth**2 * gxxi(i,j) + dydth**2 * gyyi(i,j) + &
+            dzdth**2 * gzzi(i,j) + &
+            two * ( dxdth * dydth * gxyi(i,j) + &
+                    dxdth * dzdth * gxzi(i,j) + &
+                    dydth * dzdth * gyzi(i,j) )
+      gtp(i,j) = dxdth * dxdph * gxxi(i,j) + &
+            ( dxdth * dydph + dydth * dxdph ) * gxyi(i,j) + &
+            ( dxdth * dzdph + dzdth * dxdph ) * gxzi(i,j) + &
+            dydth * dydph * gyyi(i,j) + &
+            ( dydth * dzdph + dzdth * dydph ) * gyzi(i,j) + &
+            dzdth * dzdph * gzzi(i,j)
+      gpp(i,j) = dxdph**2 * gxxi(i,j) + dydph**2 * gyyi(i,j) + &
+            dzdph**2 * gzzi(i,j) + &
+            two * ( dxdph * dydph * gxyi(i,j) + &
+                    dxdph * dzdph * gxzi(i,j) + &
+                    dydph * dzdph * gyzi(i,j) )
+      da(i,j) = sqrt ( gtt(i,j) * gpp(i,j) - gtp(i,j)**2 )
+    end do
+  end do
+!  print*,da
+!  do j = 1, lsh(2)
+!    print*,'phi = ',cphi(1,j)
+!    print*,da(:,j)
+!    print*,gtt(:,j)
+!    print*,gtp(:,j)
+!    print*,gpp(:,j)
+!  end do
+ 
+  int_tmp = weights * da
+
+  call CCTK_VarIndex ( int_var, "ehfinder::int_tmp" )
+  if ( int_var .lt. 0 ) then
+    Call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
+  end if
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     eh_area2, 1, int_var )
+  eh_area2 = eh_area2 * dtheta * dphi
+  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2
+  call CCTK_INFO(info_message)
+!  print*,eh_area2
+  
+end subroutine EHFinder_Integrate2
+
+subroutine EHFinder_Init_Weights(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k, im, jm
+
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+  end if
+
+  weights = one
+  do j = 1, lsh(2)
+    if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then
+      weights(:,j) = onethird
+    else if ( mod(lbnd(2)+j,2) .eq. 0 ) then
+      weights(:,j) = fourthirds
+    else
+      weights(:,j) = twothirds
+    end if
+    do i = 1, lsh(1)
+      if ( ( lbnd(1)+i .eq. 1 ) .or. ( lbnd(1)+i .eq. ntheta ) ) then
+        weights(i,j) = onethird * weights(i,j)
+      else if ( mod(lbnd(1)+i,2) .eq. 0 ) then
+        weights(i,j) = fourthirds * weights(i,j)
+      else
+        weights(i,j) = twothirds * weights(i,j)
+      end if
+    end do
+  end do
+end subroutine EHFinder_Init_Weights
diff --git a/src/make.code.defn b/src/make.code.defn
index 2d63268..c54b696 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -15,6 +15,7 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_SetSym.F90 \
        EHFinder_Check.F90 \
        EHFinder_Integrate.F90 \
+       EHFinder_Integrate2.F90 \
        EHFinder_FindSurface.F90 \
        EHFinder_ReadData.F90
 
diff --git a/src/make.code.deps b/src/make.code.deps
index c487d7d..1c09018 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -13,4 +13,5 @@ EHFinder_SetSym.F90.o:		EHFinder_mod.F90.o
 EHFinder_ReadData.F90.o:	EHFinder_mod.F90.o
 EHFinder_Check.F90.o:		EHFinder_mod.F90.o
 EHFinder_Integrate.F90.o:	EHFinder_mod.F90.o
+EHFinder_Integrate2.F90.o:	EHFinder_mod.F90.o
 EHFinder_FindSurface.F90.o:	EHFinder_mod.F90.o