add compuatation and printing of mean radius for each AH found

output format for BH info in stdout log is now

INFO (AHFinderDirect): AH 1/1: r=1.86525 at (-0.000000,0.000000,-0.000000)
INFO (AHFinderDirect): AH 1/1: area=45.23876324 m_irreducible=0.9486815054

and in BH_diagnostics file is now

# apparent horizon 1 of 1
#
# column 1 = cctk_iteration
# column 2 = cctk_time
# column 3 = centroid_x
# column 4 = centroid_y
# column 5 = centroid_z
# column 6 = mean radius
# column 7 = area
# column 8 = m_irreducible


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@879 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 23 insertions, 13 deletions

diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index bfeedef..6ff104b 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -219,18 +219,21 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
 				= AH_info.BH_diagnostics;
 
 			// print BH diagnostics?
-			if (verbose_info.print_physics_details)
-			   then {
-				CCTK_VInfo(CCTK_THORNSTRING,
-					   "AH found: A=%.10g at (%f,%f,%f)",
-					   double(BH_diagnostics.area),
-					   double(BH_diagnostics.centroid_x),
-					   double(BH_diagnostics.centroid_y),
-					   double(BH_diagnostics.centroid_z));
-				CCTK_VInfo(CCTK_THORNSTRING,
-					   "m_irreducible=%.10g",
-					  double(BH_diagnostics.m_irreducible));
-				}
+if (verbose_info.print_physics_details)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: r=%g at (%f,%f,%f)",
+		   hn, state.N_horizons,
+		   double(BH_diagnostics.mean_radius),
+		   double(BH_diagnostics.centroid_x),
+		   double(BH_diagnostics.centroid_y),
+		   double(BH_diagnostics.centroid_z));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
+		   hn, state.N_horizons,
+		   double(BH_diagnostics.area),
+		   double(BH_diagnostics.m_irreducible));
+	}
 
 			// store BH diagnostics in Cactus array variables
 			centroid_x[hn]    = BH_diagnostics.centroid_x;
@@ -519,7 +522,7 @@ Jacobian& Jac = *Jac_ptr;
 
 if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "searching for horizon %d of %d",
+		   "searching for horizon %d/%d",
 		   hn, N_horizons);
 
 if (timer_handle >= 0)
@@ -572,6 +575,8 @@ void compute_BH_diagnostics
 //
 fp integral_one = surface_integral(ps, gfns::gfn__one,
 				   BH_diagnostics_info.surface_integral_method);
+fp integral_h = surface_integral(ps, gfns::gfn__h,
+				 BH_diagnostics_info.surface_integral_method);
 fp integral_x = surface_integral(ps, gfns::gfn__global_x,
 				 BH_diagnostics_info.surface_integral_method);
 fp integral_y = surface_integral(ps, gfns::gfn__global_y,
@@ -588,6 +593,7 @@ case patch_system::patch_system__full_sphere:
 	break;
 case patch_system::patch_system__plus_z_hemisphere:
 	integral_one *= 2.0;
+	integral_h *= 2.0;
 	integral_x *= 2.0;
 	integral_y *= 2.0;
 	integral_z = integral_one * ps.origin_z();
@@ -595,12 +601,14 @@ case patch_system::patch_system__plus_z_hemisphere:
 case patch_system::patch_system__plus_xy_quadrant_mirrored:
 case patch_system::patch_system__plus_xy_quadrant_rotating:
 	integral_one *= 4.0;
+	integral_h *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
 	integral_z *= 4.0;
 	break;
 case patch_system::patch_system__plus_xz_quadrant_rotating:
 	integral_one *= 4.0;
+	integral_h *= 4.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y *= 4.0;
 	integral_z = integral_one * ps.origin_z();
@@ -608,6 +616,7 @@ case patch_system::patch_system__plus_xz_quadrant_rotating:
 case patch_system::patch_system__plus_xyz_octant_mirrored:
 case patch_system::patch_system__plus_xyz_octant_rotating:
 	integral_one *= 8.0;
+	integral_h *= 8.0;
 	integral_x = integral_one * ps.origin_x();
 	integral_y = integral_one * ps.origin_y();
 	integral_z = integral_one * ps.origin_z();
@@ -626,6 +635,7 @@ BH_diagnostics.centroid_y = integral_y / integral_one;
 BH_diagnostics.centroid_z = integral_z / integral_one;
 
 BH_diagnostics.area = integral_one;
+BH_diagnostics.mean_radius = integral_h / integral_one;
 BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
 }
 	  }