add support for computing BH diagnostics (mass, area, centroid)

via surface integrals over the horizon


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

2 files changed, 258 insertions, 76 deletions

diff --git a/param.ccl b/param.ccl
index f899a1c..c7701c0 100644
--- a/param.ccl
+++ b/param.ccl
@@ -3,6 +3,14 @@
 
 ################################################################################
 
+# we may need to look at grid::domain to choose our patch system symmetries
+shares: grid
+USES KEYWORD domain
+
+# all remaining parameters are private to this thorn
+private:
+################################################################################
+
 #
 # overall parameters for the apparent horizon finding algorithm itself
 #
@@ -19,22 +27,27 @@ boolean find_AHs_at_poststep \
 
 keyword method "top-level method used to find the apparent horizon"
 {
-# options mostly for testing/debugging
+# these options are mostly for testing/debugging
 "horizon function"  :: "evaluate the LHS function H(h)"
 "Jacobian test"	    :: \
   "compute/print the J[H(h)] Jacobian matrix by all possible methods"
 "Jacobian test (NP only)" :: \
   "compute/print the J[H(h)] Jacobian matrix by numerical perturbation only"
-# normal AH finding
+
+# this is for normal apparent horizon finding
 "Newton solve"	    :: "find the horizon via Newton's method"
 } "Newton solve"
 
 #
-# at present we support up to 4 horizons
+# At present we support up to 4 horizons; for user-friendliness we
+# number these 1, 2, 3, ... .   Since Cactus arrays are 0-origin, we
+# make the arrays be of size N_horizons+1, and don't use the [0] array
+# elements.
 #
-# to increase this, just raise the upper limit for this variable,
-# and change all occurences of "[4]" in this file to whatever the new
-# upper limit is (no changes are needed in the source code)
+# To change the N_horizons limit, just change the upper limit for
+# N_horizons itself, change all the [N_horizons+1] array sizes in this
+# paramter file, and recompile your configuration.  No changes are
+# needed to the source code.
 #
 int N_horizons "number of apparent horizons to search for"
 {
@@ -51,14 +64,22 @@ int N_horizons "number of apparent horizons to search for"
 keyword verbose_level \
   "controls which (how many) messages to print describing AH finding"
 {
+# 1 line each time step giving number of horizons found and their masses
 "physics highlights"   :: "just a few physics messages"
+
+# 1 line for each horizon giving position/mass/area, + a summary line or two
 "physics details"      :: "more detailed physics messages"
+
+# 1 line giving H(h) norms at each Newton iteration
 "algorithm highlights" :: \
   "physics details + a few messages about the AH-finding algorithm"
+
+# lots of details tracing what the code is doing
 "algorithm details"    :: \
   "physics details + lots of messages about the AH-finding algorithm"
 } "physics details"
 
+# n.b. printing timing stats is independent of  verbose_level
 boolean print_timing_stats \
   "should we print timing stats for the whole apparent-horizon-finding process?"
 {
@@ -76,6 +97,8 @@ keyword Jacobian_method "how do we compute the Jacobian matrix?"
   "*very* slow, but useful for debugging"
 "symbolic differentiation with finite diff d/dr" :: \
   "fast, tricky programming, uses only gij, dx gij, Kij"
+
+# alas, this isn't implemented yet :(
 "symbolic differentiation"  :: \
   "fast, tricky programming, uses gij, dx gij, dxx gij, Kij, dx Kij"
 } "symbolic differentiation with finite diff d/dr"
@@ -191,21 +214,23 @@ string Jacobian_base_file_name "base file name for Jacobian output file(s)"
 #
 private:
 
-real origin_x[4] "global x coordinate of patch system origin"
+real origin_x[5] "global x coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
-real origin_y[4] "global y coordinate of patch system origin"
+real origin_y[5] "global y coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
-real origin_z[4] "global z coordinate of patch system origin"
+real origin_z[5] "global z coordinate of patch system origin"
 {
 *:* :: "any real number"
 } 0.0
 
 keyword patch_system_type "what type of patch system should we use?"
 {
+"match Cactus grid symmetry" :: \
+  "choose automagically based on grid symmetries and the patch system's origin"
 "full sphere"	:: "full sphere, no symmetries"
 "+z hemisphere"	:: "mirror symmetry across z=0 plane"
 "+xy quadrant"	:: "90 degree periodic rotation symmetry about z axis"
@@ -213,7 +238,7 @@ keyword patch_system_type "what type of patch system should we use?"
 		    180 degree periodic rotation symmetry about z axis"
 "+xyz octant"	:: "mirror symmetry across z=0 plane *and* \
 		    90 degree periodic rotation symmetry about z axis"
-} "full sphere"
+} "match Cactus grid symmetry"
 
 int N_ghost_points "number of ghost zones on each side of a patch"
 {
@@ -226,10 +251,38 @@ int N_overlap_points \
 1:*:2 :: "any integer >= 0; current implementation requires that it be odd"
 } 1
 
+#
+# In practice the error in the horizon position is usually dominated
+# by the errors from interpolating the Cactus gij and Kij to the horizon
+# position, not by the angular finite differencing or interpatch interpolation
+# errors.  Thus this parameter can be made quite large (low resolution)
+# for better performance, without seriously affecting the accuracy
+# with which we can locate the horizon.
+#
 real delta_drho_dsigma "angular grid spacing of patches, in degrees"
 {
 (0.0:* :: "any real number > 0.0"
-} 5.0
+} 9.0
+
+################################################################################
+
+#
+# parameters for how we compute surface integrals over the horizon
+#
+
+# ... N is the number of grid zones in a patch, in either direction
+keyword surface_integral_method \
+  "how do we compute surface integrals over the horizon?"
+{
+"trapezoid"		:: ""
+"trapezoid rule"	:: "trapezoid rule (2nd order for smooth functions)"
+"Simpson"		:: ""
+"Simpson's rule"	:: \
+  "Simpson's rule (4th order for smooth fns, requires N to be even)"
+"Simpson (variant)"	:: ""
+"Simpson's rule (variant)":: \
+  "Simpson's rule variant (4th order for smooth fns, requires N >= 7)"
+} "trapezoid"
 
 ################################################################################
 
@@ -242,6 +295,7 @@ keyword geometry_method "how do we compute the slice's geometry?"
 {
 # for normal use
 "interpolate from Cactus grid" :: "interpolate gij and Kij from Cactus grid"
+
 # for testing/debugging
 "Schwarzschild/EF"             :: \
   "hard-wire to Schwarzschild spacetime / Eddington-Finkelstein slice"
@@ -252,6 +306,10 @@ keyword geometry_method "how do we compute the slice's geometry?"
 #
 # parameters for geometry_method = "interpolate from Cactus grid"
 #
+# Note that the interpolated gij and Kij should be at least C1, otherwise
+# the Newton iteration may fail to converge all the way down to tight
+# error tolerances.  In practice a Hermite interpolant works well.
+#
 
 string geometry_interpolator_name \
   "name under which the geometry interpolation operator is registered in Cactus"
@@ -320,7 +378,12 @@ real geometry__Schwarzschild_EF__Delta_xyz \
 ################################################################################
 
 #
-# parameters for the interpatch interpolator
+# parameters for the (1-D angular) interpatch interpolator
+#
+# Note there's no necessary relationship between this interpolator and
+# the geometry interpolator.  In particular, this interpolator could
+# reasonably use a different interpolation operator and/or order from
+# the geometry interpolator.
 #
 string interpatch_interpolator_name \
   "name under which the interpatch interpolation operator is registered in Cactus"
@@ -344,10 +407,11 @@ keyword initial_guess_method \
   "method used to set up initial guess for apparent horizon shape"
 {
 "read from file"   :: "read from input file"
-"sphere"           :: "set up a coordinate sphere"
+"Kerr/Kerr"	   :: "set to the horizon of Kerr spacetime in Kerr coordinates"
+"Kerr/Kerr-Schild" :: \
+  "set to the horizon of Kerr spacetime in Kerr-Schild coordinates"
 "ellipsoid"	   :: "set to a coordinate ellipsoid"
-"Kerr/Kerr"	   :: "set to horizon of Kerr spacetime in Kerr coords"
-"Kerr/Kerr-Schild" :: "set to horizon of Kerr spacetime in Kerr-Schild coords"
+"sphere"           :: "set to a coordinate sphere"
 } "read from file"
 
 boolean output_initial_guess \
@@ -356,69 +420,69 @@ boolean output_initial_guess \
 } "true"
 
 # parameters for initial_guess_method = "sphere"
-real initial_guess__sphere__x_center[4] "x coordinate of sphere center"
+real initial_guess__sphere__x_center[5] "x coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__y_center[4] "y coordinate of sphere center"
+real initial_guess__sphere__y_center[5] "y coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__z_center[4] "z coordinate of sphere center"
+real initial_guess__sphere__z_center[5] "z coordinate of sphere center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__sphere__radius[4] "radius of sphere"
+real initial_guess__sphere__radius[5] "radius of sphere"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
 
 # parameters for initial_guess_method = "ellipsoid"
-real initial_guess__ellipsoid__x_center[4] "x coordinate of ellipsoid center"
+real initial_guess__ellipsoid__x_center[5] "x coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__y_center[4] "y coordinate of ellipsoid center"
+real initial_guess__ellipsoid__y_center[5] "y coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__z_center[4] "z coordinate of ellipsoid center"
+real initial_guess__ellipsoid__z_center[5] "z coordinate of ellipsoid center"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__ellipsoid__x_radius[4] "x radius of ellipsoid"
+real initial_guess__ellipsoid__x_radius[5] "x radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
-real initial_guess__ellipsoid__y_radius[4] "y radius of ellipsoid"
+real initial_guess__ellipsoid__y_radius[5] "y radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
-real initial_guess__ellipsoid__z_radius[4] "z radius of ellipsoid"
+real initial_guess__ellipsoid__z_radius[5] "z radius of ellipsoid"
 {
 (0.0:* :: "any real number > 0.0"
 } 2.0
 
 # parameters for initial_guess_method = "Kerr/Kerr" and "Kerr/Kerr-Schild"
-real initial_guess__Kerr__x_posn[4] "x coordinate of Kerr BH"
+real initial_guess__Kerr__x_posn[5] "x coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__Kerr__y_posn[4] "y coordinate of Kerr BH"
+real initial_guess__Kerr__y_posn[5] "y coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
-real initial_guess__Kerr__z_posn[4] "z coordinate of Kerr BH"
+real initial_guess__Kerr__z_posn[5] "z coordinate of Kerr BH"
 {
 *:* :: "any real number"
 } 0.0
 # n.b. my convention is that a=J/m^2 is dimensionless,
 #      while MTW take a=J/m=m * (my a)
-real initial_guess__Kerr__mass[4] "mass of Kerr BH"
+real initial_guess__Kerr__mass[5] "mass of Kerr BH"
 {
 (0.0:*		:: "BH mass = any real number > 0"
 } 1.0
-real initial_guess__Kerr__spin[4] "dimensionless spin J/m^2 of Kerr BH"
+real initial_guess__Kerr__spin[5] "dimensionless spin J/m^2 of Kerr BH"
 {
 (-1.0:1.0)	:: "BH spin = J/m^2 = any real number with absolute value < 1"
 } 0.6
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index e12ecf1..3f3279e 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -3,7 +3,7 @@
 //
 // <<<access to persistent data>>>
 // <<<prototypes for functions local to this file>>>
-// AHFinderDirect_find_horizons - top-level driver
+// AHFinderDirect_find_horizons - top-level driver to find apparent horizons
 //
 
 #include <stdio.h>
@@ -38,7 +38,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
-#include "AHFinderDirect.hh"
+#include "driver.hh"
 
 //******************************************************************************
 
@@ -54,13 +54,16 @@ extern struct state state;
 //
 namespace {
 bool find_horizon(enum method method,
-		  enum verbose_level verbose_level, int timer_handle,
+		  const struct verbose_info& verbose_info, int timer_handle,
 		  struct IO_info& IO_info,
 		  struct Jacobian_info& Jac_info,
 		  struct solver_info& solver_info,
 		  struct cactus_grid_info& cgi, struct geometry_info& gi,
 		  patch_system& ps, Jacobian* Jac_ptr,
-		  int hn);
+		  int hn, int N_horizons);
+void BH_diagnostics(enum patch::integration_method surface_integral_method,
+		    const struct verbose_info& verbose_info,
+		    struct AH_info& AH_info);
 	  };
 
 //******************************************************************************
@@ -74,25 +77,40 @@ extern "C"
 {
 DECLARE_CCTK_ARGUMENTS
 DECLARE_CCTK_PARAMETERS
+const struct verbose_info& verbose_info = state.verbose_info;
 
 state.IO_info.time_iteration = cctk_iteration;
 
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
 
-	for (int hn = 0 ; hn < state.N_horizons ; ++hn)
+	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
 	{
 	struct AH_info& AH_info = * state.AH_info_ptrs[hn];
+	patch_system& ps = *AH_info.ps_ptr;
 
 	AH_info.AH_found
 	   = find_horizon(state.method,
-			  state.verbose_level, state.timer_handle,
+			  verbose_info, state.timer_handle,
 			  state.IO_info, state.Jac_info, state.solver_info,
 			  state.cgi, state.gi,
-			  * AH_info.ps_ptr, AH_info.Jac_ptr,
-			  hn);
-
-	// FIXME FIXME: compute the centroid, area, mass, ...
+			  ps, AH_info.Jac_ptr,
+			  hn, state.N_horizons);
+
+	if (AH_info.AH_found)
+	   then {
+		BH_diagnostics(state.surface_integral_method,
+			       verbose_info,
+			       AH_info);
+
+		if (verbose_info.print_physics_details)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "AH found: A=%g m=%g at (%g,%g,%g)",
+				   double(AH_info.area), double(AH_info.mass),
+				   double(AH_info.centroid_x),
+				   double(AH_info.centroid_y),
+				   double(AH_info.centroid_z));
+		}
 	}
 
 if (state.timer_handle >= 0)
@@ -114,27 +132,26 @@ if (state.timer_handle >= 0)
 //		  (we only time the computation, not the file I/O)
 // Jac_ptr = may be NULL if no Jacobian is needed (depending on  method)
 // hn = the horizon number (used only in naming output files)
+// N_horizons = the total number of horizon(s) being searched for number
+//		(used only in formatting info messages)
 //
 // Results:
-// This function returns a success flag: this is true if (and only if)
-// it found an h satisfying H(h) = 0 to within the error tolerances,
-// otherwise it's false.
+// This function returns true if the computation succeeds; false if it fails.
+// If  method  specifies finding the (an) apparent horizon, then "success"
+// means finding an h satisfying H(h) = 0 to within the error tolerances.
+// Otherwise, "success" means successfully evaluating the horizon function
+// and/or its Jacobian, as appropriate.
 //
 namespace {
 bool find_horizon(enum method method,
-		  enum verbose_level verbose_level, int timer_handle,
+		  const struct verbose_info& verbose_info, int timer_handle,
 		  struct IO_info& IO_info,
 		  struct Jacobian_info& Jac_info,
 		  struct solver_info& solver_info,
 		  struct cactus_grid_info& cgi, struct geometry_info& gi,
 		  patch_system& ps, Jacobian* Jac_ptr,
-		  int hn)
+		  int hn, int N_horizons)
 {
-bool success = false;
-
-const bool print_io_msg_flag
-	= verbose_level >= verbose_level__algorithm_highlights;
-
 switch	(method)
 	{
 // just evaluate the horizon function
@@ -144,60 +161,70 @@ case method__horizon_function:
 
 	if (timer_handle >= 0)
 	   then CCTK_TimerStartI(timer_handle);
-	horizon_function(ps, cgi, gi, false, true, &H_norms);
+	const bool status
+		= horizon_function(ps, cgi, gi, false, true, &H_norms);
 	if (timer_handle >= 0)
 	   then CCTK_TimerStopI(timer_handle);
+	if (!status)
+	   then return false;				// *** ERROR RETURN ***
 
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   H(h) rms-norm %.2e, infinity-norm %.2e",
-		   H_norms.rms_norm(), H_norms.infinity_norm());
+	if (H_norms.is_nonempty())	// might be empty iv H(h) eval failed
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   H(h) rms-norm %.2e, infinity-norm %.2e",
+			   H_norms.rms_norm(), H_norms.infinity_norm());
 	output_gridfn(ps, gfns::gfn__H,
 		      IO_info, IO_info.H_base_file_name,
 		      hn, true);
-	break;
+	return true;					// *** NORMAL RETURN ***
 	  }
 
 // just compute/print the NP Jacobian
 case method__Jacobian_test_NP_only:
 	  {
 	Jacobian& Jac_NP = *Jac_ptr;
-	horizon_function(ps, cgi, gi, true);
+	if (! horizon_function(ps, cgi, gi, true))
+	   then return false;				// *** ERROR RETURN ***
 	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
-	horizon_Jacobian(ps, Jac_NP,
-			 cgi, gi, Jac_info,
-			 true);
+	if (!  horizon_Jacobian(ps, Jac_NP,
+				cgi, gi, Jac_info,
+				true))
+	   then return false;				// *** ERROR RETURN ***
 
 	print_Jacobians(ps,
 			& Jac_NP, NULL,
 			IO_info, IO_info.Jacobian_base_file_name,
 			hn, true);
-	break;
+	return true;					// *** NORMAL RETURN ***
 	  }
 
 // compute/print the Jacobian by all possible methods
 case method__Jacobian_test:
 	  {
 	Jacobian& Jac_NP = *Jac_ptr;
-	horizon_function(ps, cgi, gi, true);
+	if (!  horizon_function(ps, cgi, gi, true))
+	   then return false;				// *** ERROR RETURN ***
 	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
-	horizon_Jacobian(ps, Jac_NP,
-			 cgi, gi, Jac_info,
-			 true);
+	if (! horizon_Jacobian(ps, Jac_NP,
+			       cgi, gi, Jac_info,
+			       true))
+	   then return false;				// *** ERROR RETURN ***
 
 	// symbolic differentiation with finite diff d/dr
 	Jacobian& Jac_SD_FDdr
 		= new_Jacobian(ps, Jac_info.Jacobian_storage_method);
-	horizon_function(ps, cgi, gi, true);
+	if (! horizon_function(ps, cgi, gi, true))
+	   then return false;				// *** ERROR RETURN ***
 	Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr;
-	horizon_Jacobian(ps, Jac_SD_FDdr,
-			 cgi, gi, Jac_info,
-			 true);
+	if (! horizon_Jacobian(ps, Jac_SD_FDdr,
+			       cgi, gi, Jac_info,
+			       true))
+	   then return false;				// *** ERROR RETURN ***
 
 	print_Jacobians(ps,
 			& Jac_NP, & Jac_SD_FDdr,
 			IO_info, IO_info.Jacobian_base_file_name,
 			hn, true);
-	break;
+	return true;					// *** NORMAL RETURN ***
 	  }
 
 // find the apparent horizon via the Newton solver
@@ -205,22 +232,30 @@ case method__Newton_solve:
 	  {
 	Jacobian& Jac = *Jac_ptr;
 
+	if (verbose_info.print_algorithm_highlights)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "searching for horizon #%d/%d",
+			   hn, N_horizons);
+
 	if (timer_handle >= 0)
 	   then CCTK_TimerStartI(timer_handle);
-	success = Newton_solve(ps, Jac,
+	const bool status
+		= Newton_solve(ps, Jac,
 			       cgi, gi,
 			       Jac_info, solver_info, IO_info,
-			       hn, verbose_level);
+			       hn, verbose_info);
 	if (timer_handle >= 0)
 	   then CCTK_TimerStopI(timer_handle);
+	if (! status)
+	   then return false;				// *** ERROR RETURN ***
 
 	output_gridfn(ps, gfns::gfn__h,
 		      IO_info, IO_info.h_base_file_name,
-		      hn, print_io_msg_flag);
+		      hn, verbose_info.print_algorithm_details);
 	output_gridfn(ps, gfns::gfn__H,
 		      IO_info, IO_info.H_base_file_name,
-		      hn, print_io_msg_flag);
-	break;
+		      hn, verbose_info.print_algorithm_details);
+	return true;					// *** NORMAL RETURN ***
 	  }
 
 default:
@@ -229,9 +264,92 @@ default:
 "   find_horizons(): unknown method=(int)%d!\n"
 "                    (this should never happen!)"
 ,
-		   (int) method);				/*NOTREACHED*/
+		   int(method));				/*NOTREACHED*/
 	}
 
-return success;
+/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+
+//
+// Given that an apparent horizon has been found, this function computes
+// various BH diagnostics.
+//
+// Inputs (gridfns)
+// h		# ghosted
+// one		# nominal
+// global_[xyz]	# nominal
+//
+namespace {
+void BH_diagnostics(enum patch::integration_method surface_integral_method,
+		    const struct verbose_info& verbose_info,
+		    struct AH_info& AH_info)
+{
+const patch_system& ps = * AH_info.ps_ptr;
+
+//
+// compute raw surface integrals
+//
+fp integral_one = ps.integrate_gridfn(gfns::gfn__one,
+				      true,	// area weighting
+				      surface_integral_method);
+fp integral_x = ps.integrate_gridfn(gfns::gfn__global_x,
+				    true,		// area weighting
+				    surface_integral_method);
+fp integral_y = ps.integrate_gridfn(gfns::gfn__global_y,
+				    true,		// area weighting
+				    surface_integral_method);
+fp integral_z = ps.integrate_gridfn(gfns::gfn__global_z,
+				    true,		// area weighting
+				    surface_integral_method);
+
+//
+// adjust integrals to take into account patch system not covering full sphere
+//
+switch	(ps.type())
+	{
+case patch_system::full_sphere_patch_system:
+	break;
+case patch_system::plus_z_hemisphere_patch_system:
+	integral_one *= 2.0;
+	integral_x *= 2.0;
+	integral_y *= 2.0;
+	integral_z = integral_one * ps.origin_z();
+	break;
+case patch_system::plus_xy_quadrant_patch_system:
+	integral_one *= 4.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y = integral_one * ps.origin_y();
+	integral_z *= 4.0;
+	break;
+case patch_system::plus_xz_quadrant_patch_system:
+	integral_one *= 4.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y *= 4.0;
+	integral_z = integral_one * ps.origin_z();
+	break;
+case patch_system::plus_xyz_octant_patch_system:
+	integral_one *= 8.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y = integral_one * ps.origin_y();
+	integral_z = integral_one * ps.origin_z();
+	break;
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   BH_diagnostics(): unknown ps.type()=(int)%d!\n"
+"                     (this should never happen!)"
+,
+		   int(ps.type()));				/*NOTREACHED*/
+	}
+
+AH_info.area = integral_one;
+AH_info.mass = sqrt(AH_info.area / (16.0*PI));
+
+AH_info.centroid_x = integral_x / integral_one;
+AH_info.centroid_y = integral_y / integral_one;
+AH_info.centroid_z = integral_z / integral_one;
 }
 	  }