Added code for calculating circumferences. Changed storing of results, since

I switched to 1D grid array with default distribution.


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

3 files changed, 276 insertions, 78 deletions

diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index ece1b1c..655109b 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -205,13 +205,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       s_symz = .true.; s_symy = .true.
     else if ( crange_min_glob(3) .lt. zero ) then
       ltheta = zero; utheta = half * pi
-      lphi = 0; uphi = two * pi
+      lphi = zero; uphi = two * pi
       sym_factor = two
       center(3) = zero
       s_symz = .true.
     else if ( crange_min_glob(2) .lt. zero ) then
       ltheta = zero; utheta = pi
-      lphi = 0; uphi = pi
+      lphi = zero; uphi = pi
       sym_factor = two
       center(2) = zero
       s_symy = .true.
@@ -229,7 +229,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       s_symz = .true.; s_symx = .true.
     else if ( crange_min_glob(3) .lt. zero ) then
       ltheta = zero; utheta = half * pi
-      lphi = 0; uphi = two * pi
+      lphi = zero; uphi = two * pi
       sym_factor = two
       center(3) = zero
       s_symz = .true.
@@ -253,7 +253,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       s_symy = .true.; s_symx = .true.
     else if ( crange_min_glob(2) .lt. zero ) then
       ltheta = zero; utheta = pi
-      lphi = 0; uphi = pi
+      lphi = zero; uphi = pi
       sym_factor = two
       center(2) = zero
       s_symy = .true.
@@ -318,6 +318,21 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     end if
   end if
 
+  if ( .not. s_symz ) then
+    githeta = ( ntheta - 1 ) / 2 + 1
+  else
+    githeta = ntheta
+  end if
+  
+  if ( ( .not. s_symy ) .and. ( s_symx ) ) then
+    gjphi = ( nphi - 1 ) / 2 + 1
+  else
+    gjphi = 1
+  end if
+
+  theta_sym_factor = two * pi / ( utheta - ltheta )
+  phi_sym_factor = two * pi / ( uphi - lphi )
+
 !  print*
 !  print*,ltheta/pi,utheta/pi,lphi/pi,uphi/pi
 !  print*,sym_factor
@@ -633,6 +648,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 !  print*,'Calculated derivatives'
 !  print*,rsurf
   call Util_TableDestroy ( status, table_handle )
+  print*,'gjphi = ',gjphi
   
 end subroutine EHFinder_FindSurface
 
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index 53848f9..0b06ffb 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -25,8 +25,8 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   CCTK_REAL :: dxdth, dxdph, dydth, dydph, dzdth, dzdph
 !  CCTK_REAL :: gtt, gtp, gpp
   CCTK_POINTER, dimension(3) :: interp_coords
-  CCTK_POINTER, dimension(7) :: out_array
-  CCTK_INT, dimension(7) :: in_array
+  CCTK_POINTER, dimension(7) :: in_array, out_array
+!  CCTK_INT, dimension(7) :: in_array
   CCTK_INT, dimension(7), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL, &
                                                       CCTK_VARIABLE_REAL, &
@@ -176,6 +176,8 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
                     dxdph * dzdph * gxzi(i,j) + &
                     dydph * dzdph * gyzi(i,j) )
       da(i,j) = dtheta * dphi * sqrt ( gtt(i,j) * gpp(i,j) - gtp(i,j)**2 )
+      dltheta(i,j) = dtheta * sqrt ( gtt(i,j) )
+      dlphi(i,j) = dphi * sqrt ( gpp(i,j) )
     end do
   end do
 !  print*,da
@@ -203,7 +205,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   CCTK_INT :: int_var, sn
   CCTK_REAL :: eh_area_tmp
   character(len=30) :: info_message
-!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
   call CCTK_INFO ( 'Integrating area' )
 
@@ -212,35 +214,37 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
 
   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' )
+    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(sn), 1, int_var )
-  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn)
-
 !  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-!                     eh_area_tmp, 1, int_var )
-!  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp
+!                     eh_area2(sn), 1, int_var )
+!  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn)
+
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     eh_area_tmp, 1, int_var )
+  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp
 
   call CCTK_INFO(info_message)
 
-!  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
-!  end if
-!  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
-!  end if
-!  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get local size for area array" )
-!  end if
-!
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for area array" )
+  end if
+
 !  if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
-!    eh_area2(sn-lbnd(1)) = eh_area_tmp
-!  end if
-!  print*,lbnd,ubnd,lsh
+  if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then
+    eh_area2(sn-lbnd(1)) = eh_area_tmp
+  end if
+  print*,'Debug information : ',sn,lbnd,ubnd,lsh
+!  print*,eh_area2
 
 end subroutine EHFinder_IntegrateArea
 
@@ -257,7 +261,7 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
 
   CCTK_INT :: int_var, sn, k
   character(len=50) :: info_message
-!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
   CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
   call CCTK_INFO ( 'Integrating centroid' )
@@ -270,74 +274,221 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
 
   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' )
+    call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
   end if
 
+!  int_tmp = sym_factor * interp_x * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     eh_centroid2_x(sn), 1, int_var )
+!
+!  int_tmp = sym_factor * interp_y * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     eh_centroid2_y(sn), 1, int_var )
+!
+!  int_tmp = sym_factor * interp_z * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     eh_centroid2_z(sn), 1, int_var )
+
+!  if ( s_symx ) eh_centroid2_x(sn) = zero
+!  if ( s_symy ) eh_centroid2_y(sn) = zero
+!  if ( s_symz ) eh_centroid2_z(sn) = zero
+
   int_tmp = sym_factor * interp_x * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid2_x(sn), 1, int_var )
+                     centroid_x, 1, int_var )
 
   int_tmp = sym_factor * interp_y * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid2_y(sn), 1, int_var )
+                     centroid_y, 1, int_var )
 
   int_tmp = sym_factor * interp_z * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid2_z(sn), 1, int_var )
+                     centroid_z, 1, int_var )
 
-  if ( s_symx ) eh_centroid2_x(sn) = zero
-  if ( s_symy ) eh_centroid2_y(sn) = zero
-  if ( s_symz ) eh_centroid2_z(sn) = zero
+  if ( s_symx ) centroid_x = zero
+  if ( s_symy ) centroid_y = zero
+  if ( s_symz ) centroid_z = zero
 
-!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-!                     centroid_x, 1, int_var )
-!
-!  int_tmp = interp_y * weights * da
-!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-!                     centroid_y, 1, int_var )
-!
-!  int_tmp = interp_z * weights * da
-!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-!                     centroid_z, 1, int_var )
-!
-!  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
-!  end if
-!  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
-!  end if
-!  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_area2" )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN ( 0, "cannot get local size for area array" )
-!  end if
-
-    eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn)
-    eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn)
-    eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn)
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_centroid2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_centroid2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_centroid2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for area array" )
+  end if
 
-!  if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
-!    k = sn - lbnd(1)
-!    eh_centroid2_x(k) = eh_centroid2_x(k) / eh_area2(k)
-!    eh_centroid2_y(k) = eh_centroid2_y(k) / eh_area2(k)
-!    eh_centroid2_z(k) = eh_centroid2_z(k) / eh_area2(k)
-!  end if
+!    eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn)
+!    eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn)
+!    eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn)
 
-  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
-                                       eh_centroid2_x(sn), &
-                                       eh_centroid2_y(sn), &
-                                       eh_centroid2_z(sn)
+  if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then
+    k = sn - lbnd(1)
+    eh_centroid2_x(k) = centroid_x / eh_area2(k)
+    eh_centroid2_y(k) = centroid_y / eh_area2(k)
+    eh_centroid2_z(k) = centroid_z / eh_area2(k)
+  end if
 
 !  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
-!                                       centroid_x, &
-!                                       centroid_y, &
-!                                       centroid_z
+!                                       eh_centroid2_x(sn), &
+!                                       eh_centroid2_y(sn), &
+!                                       eh_centroid2_z(sn)
+
+  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+                                       centroid_x, &
+                                       centroid_y, &
+                                       centroid_z
   call CCTK_INFO(info_message)
 
 end subroutine EHFinder_IntegrateCentroid
 
 
+subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: int_var, sn, k
+  character(len=50) :: info_message
+  CCTK_INT, dimension(1) :: lbnd1, ubnd1, lsh1
+  CCTK_INT, dimension(2) :: lbnd, lsh, nghost
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha
+  CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ
+
+  sn = surface_counter - integrate_counter
+
+  ! Find out the lower bounds of the distributed integration grid arrays.
+  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
+
+  ! Find out the local size of the distributed integration grid arrays
+  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
+
+  ! Find out the bounding box of the distributed integration grid arrays
+  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
+
+  ! Find out the ghost size of the distributed integration grid arrays
+  call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, "ehfinder::surface_arrays" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get ghost zones for surface arrays" )
+  end if
+ 
+!  print*,lbnd
+!  print*,lsh
+!  print*,bbox
+!  print*,nghost
+
+  ithl = 1; ithr = lsh(1)
+  jphl = 1; jphr = lsh(2)
+  
+  if ( bbox(1) .eq. 0 ) ithl = ithl + nghost(1)
+  if ( bbox(2) .eq. 0 ) ithr = ithr - nghost(1)
+  if ( bbox(3) .eq. 0 ) jphl = jphl + nghost(2)
+  if ( bbox(4) .eq. 0 ) jphr = jphr - nghost(2)
+
+!  print*,ithl,ithr
+!  print*,jphl,jphr
+
+  call CCTK_INFO ( 'Integrating equatorial circumference' )
+
+  int_tmp = phi_sym_factor * phiweights * dlphi
+
+  if ( ( ithl+lbnd(1) .le. githeta ) .and. ( githeta .le. ithr+lbnd(1) ) ) then
+    itha = githeta - lbnd(1)
+    eq_circ_loc = sum ( int_tmp(itha,jphl:jphr) )
+  else
+    eq_circ_loc = zero
+  end if
+
+!  print*,itha,githeta,eq_circ_loc
+!  print*
+!  print*,int_tmp(itha,jphl:jphr)
+!  print*
+!  print*,phi_sym_factor
+!  print*
+!  print*,phiweights(itha,jphl:jphr)
+!  print*
+!  print*,dlphi(itha,jphl:jphr)
+
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                              eq_circ_loc, eq_circ, CCTK_VARIABLE_REAL )
+
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd1, "ehfinder::eh_circumference2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get lower bounds for area array" )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd1, "ehfinder::eh_circumference2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get upper bounds for area array" )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh1, "ehfinder::eh_circumference2" )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, "cannot get local size for area array" )
+  end if
+
+  if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then
+    k = sn - lbnd1(1)
+    eh_circ_eq2(k) = eq_circ
+  end if
+!  eh_circ_eq2(sn) = eq_circ
+
+  write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ',eq_circ
+  call CCTK_INFO(info_message)
+
+  call CCTK_INFO ( 'Integrating polar circumference' )
+
+  int_tmp = theta_sym_factor * thetaweights * dltheta
+
+!  print*,jphl+lbnd(2),jphr+lbnd(2),gjphi
+  if ( ( jphl+lbnd(2) .le. gjphi ) .and. ( gjphi .le. jphr+lbnd(2) ) ) then
+    jpha = gjphi - lbnd(2)
+    pol_circ_loc = sum ( int_tmp(ithl:ithr,jpha) )
+  else
+    pol_circ_loc = zero
+  end if
+
+!  print*,jpha,gjphi,pol_circ_loc
+!  print*
+!  print*,int_tmp(ithl:ithr,jpha)
+!  print*
+!  print*,theta_sym_factor
+!  print*,thetaweights(ithl:ithr,jpha)
+!  print*
+!  print*,dltheta(ithl:ithr,jpha)
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                              pol_circ_loc, pol_circ, CCTK_VARIABLE_REAL )
+
+
+  if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then
+    k = sn - lbnd1(1)
+    eh_circ_pol2(k) = pol_circ
+  end if
+!  eh_circ_pol2(sn) = pol_circ
+
+  write(info_message,'(a30,f10.5)') 'Horizon polar circumference = ',pol_circ
+  call CCTK_INFO(info_message)
+  
+end subroutine EHFinder_IntegrateEquatorial
+
+
 ! This routine sets up the weights for the Simpsons rule integration
 ! over the surface.
 subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
@@ -381,10 +532,13 @@ subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
     ! the same theta coordinate are set to the same weight.
     if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then
       weights(:,j) = onethird
+      phiweights(:,j) = onethird
     else if ( mod(lbnd(2)+j,2) .eq. 0 ) then
       weights(:,j) = fourthirds
+      phiweights(:,j) = fourthirds
     else
       weights(:,j) = twothirds
+      phiweights(:,j) = twothirds
     end if
 
     ! Then it is done in the theta direction with the one-dimensional
@@ -392,10 +546,13 @@ subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
     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)
+        thetaweights(i,j) = onethird
       else if ( mod(lbnd(1)+i,2) .eq. 0 ) then
         weights(i,j) = fourthirds * weights(i,j)
+        thetaweights(i,j) = fourthirds
       else
         weights(i,j) = twothirds * weights(i,j)
+        thetaweights(i,j) = twothirds
       end if
     end do
   end do
@@ -424,6 +581,9 @@ subroutine EHFinder_CopyArea(CCTK_ARGUMENTS)
 
   call CCTK_INFO('Copying area data')
 
+!  print*,Xeh_area0
+!  print*,Xeh_area20
+
   eh_area = eh_area2
 
 end subroutine EHFinder_CopyArea
@@ -441,8 +601,29 @@ subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS)
 
   call CCTK_INFO('Copying centroid data')
 
+!  print*,Xeh_centroid0
+!  print*,Xeh_centroid20
+
   eh_centroid_x = eh_centroid2_x
   eh_centroid_y = eh_centroid2_y
   eh_centroid_z = eh_centroid2_z
 
 end subroutine EHFinder_CopyCentroid
+
+
+subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  call CCTK_INFO('Copying circumference data')
+
+  eh_circ_eq = eh_circ_eq2
+  eh_circ_pol = eh_circ_pol2
+
+end subroutine EHFinder_CopyCircumference
diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90
index 0588cbe..6d6b319 100644
--- a/src/EHFinder_mod.F90
+++ b/src/EHFinder_mod.F90
@@ -29,6 +29,7 @@ module EHFinder_mod
   CCTK_INT, dimension(3) :: imin_glob, imax_glob
   CCTK_REAL, dimension(3) :: fimin_glob, fimax_glob
   CCTK_INT :: current_iteration
+  CCTK_INT :: githeta, gjphi
   logical :: mask_first = .true.
   logical :: symx, symy, symz
   logical :: s_symx, s_symy, s_symz