Minor change in mask initialization, plus added a new output file for the

horizon radius (the c0(0) coefficient).  I have found this very useful to
get an idea of the size of the horizon.


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

1 files changed, 37 insertions, 23 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index 4aa722d..f18a890 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -48,7 +48,7 @@
       CCTK_REAL, parameter :: half = 0.5D0
       CCTK_REAL, parameter :: one  = 1.0D0
 
-      character*200 almf,logf,areaf,massf,circeqf,merip1f,merip2f
+      character*200 almf,logf,areaf,massf,radf,circeqf,merip1f,merip2f
 
       save ahfafter,ahfso,atime
 
@@ -89,6 +89,21 @@
 !                  4.  Integral of H^2 is positive out and in.
 
 
+!     ***************************
+!     ***   INITIALIZE MASK   ***
+!     ***************************
+
+!     Initialize mask to one, but only on first iteration.
+!     Notice that I do this even if we do not want to find
+!     a horizon from the very beginning.  This is to make
+!     sure that the mask is 1 and not 0 in case some other
+!     thorn wants to use it.
+
+      if (cctk_iteration.eq.0) then
+         ahmask = one
+      end if
+
+
 !     ***********************************************************
 !     ***   CHECK IF WE WANT TO FIND A HORIZON AT THIS TIME   ***
 !     ***********************************************************
@@ -112,17 +127,6 @@
 
 
 !     ***************************
-!     ***   INITIALIZE MASK   ***
-!     ***************************
-
-!     Initialize mask to one, but only on first call.
-
-      if (ahf_ncall.eq.0) then
-         ahmask = one
-      end if
-
-
-!     ***************************
 !     ***   GRID PARAMETERS   ***
 !     ***************************
 
@@ -241,6 +245,7 @@
 
             areaf = filestr(1:nfile)//"/ahf_area_0.tl"
             massf = filestr(1:nfile)//"/ahf_mass_0.tl"
+            radf  = filestr(1:nfile)//"/ahf_rad_0.tl"
 
             circeqf = filestr(1:nfile)//"/ahf_circ_eq_0.tl"
             merip1f = filestr(1:nfile)//"/ahf_meri_p1_0.tl"
@@ -253,6 +258,7 @@
 
             areaf = filestr(1:nfile)//"/ahf_area_1.tl"
             massf = filestr(1:nfile)//"/ahf_mass_1.tl"
+            radf  = filestr(1:nfile)//"/ahf_rad_1.tl"
 
             circeqf = filestr(1:nfile)//"/ahf_circ_eq_1.tl"
             merip1f = filestr(1:nfile)//"/ahf_meri_p1_1.tl"
@@ -265,6 +271,7 @@
 
             areaf = filestr(1:nfile)//"/ahf_area_2.tl"
             massf = filestr(1:nfile)//"/ahf_mass_2.tl"
+            radf  = filestr(1:nfile)//"/ahf_rad_2.tl"
 
             circeqf = filestr(1:nfile)//"/ahf_circ_eq_2.tl"
             merip1f = filestr(1:nfile)//"/ahf_meri_p1_2.tl"
@@ -279,7 +286,8 @@
             
          areaf = filestr(1:nfile)//"/ahf_area.tl"
          massf = filestr(1:nfile)//"/ahf_mass.tl"
-         
+         radf  = filestr(1:nfile)//"/ahf_rad.tl"
+
          circeqf = filestr(1:nfile)//"/ahf_circ_eq.tl"
          merip1f = filestr(1:nfile)//"/ahf_meri_p1.tl"
          merip2f = filestr(1:nfile)//"/ahf_meri_p2.tl"
@@ -1352,16 +1360,16 @@
      .   0.141047396D0*dsqrt(intarea_h)
          close(1)
 
-c         if (ahf_ncall.eq.1) then
-c            open(1,file='ahf_radius.tl',form='formatted',
-c     .      status='replace')
-c            write(1,"(A16)") '" Horizon radius'
-c         else
-c            open(1,file='ahf_radius.tl',form='formatted',
-c     .      status='old',position='append')
-c         end if
-c         if (status) write(1,"(2ES14.6)") cctk_time,c0(0)
-c         close(1)
+         if (ahf_ncall.eq.1) then
+            open(1,file=radf,form='formatted',
+     .      status='replace')
+            write(1,"(A16)") '" Horizon radius'
+         else
+            open(1,file=radf,form='formatted',
+     .      status='old',position='append')
+         end if
+         if (status) write(1,"(2ES14.6)") cctk_time,c0(0)
+         close(1)
 
 !        Circumferences.
 
@@ -1459,6 +1467,12 @@ c         close(1)
 
          call AHFinder_mask(CCTK_FARGUMENTS,rhor)
 
+!        If desired, copy ahf_mask to einstein mask (emask).
+
+c         if (use_mask.eq.1) then
+c            emask = ahf_mask
+c         end if
+
 #ifdef DEVELOPMENT_SIMPLEEXCISION
 
 !        Here I alter the parameters for SimpleExcision.  Notice that