* add interface to Erik's SphericalSurface thorn

  ==> With this commit, AHFinderDirect now inherits from
      AEIThorns/SphericalSurface, so you must have that thorn in
      your configuration to be able to compile.
* add computation of surface quadrupole moments and areal radius
* expand BH_diagnostics file format to accomodate quadrupole moments
  and areal radius, and also to include not-implemented-yet columns
  for 9 more diagnostics which Erik has implemented in his branch
* some other small cleanups


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

10 files changed, 521 insertions, 56 deletions

diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc
index a641dde..163ef49 100644
--- a/src/driver/BH_diagnostics.cc
+++ b/src/driver/BH_diagnostics.cc
@@ -7,7 +7,7 @@
 // BH_diagnostics::copy_from_buffer - copy buffer to diagnostics
 //
 // BH_diagnostics::compute - compute BH diagnostics after an AH has been found
-/// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere
+// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere
 //
 // print - print a line or two summarizing the diagnostics
 // setup_output_file - create/open output file, write header describing fields
@@ -18,6 +18,8 @@
 #include <assert.h>
 #include <math.h>
 
+#include <vector>
+
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -57,17 +59,27 @@ namespace AHFinderDirect
 //******************************************************************************
 
 //
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
 // This function initializes a  struct BH_diagnostics  to all zeros.
 //
 BH_diagnostics::BH_diagnostics()
 	: centroid_x(0.0), centroid_y(0.0), centroid_z(0.0),
-	  min_radius(0.0), max_radius(0.0), mean_radius(0.0),
+	  quadrupole_xx(0.0), quadrupole_xy(0.0), quadrupole_xz(0.0),
+			      quadrupole_yy(0.0), quadrupole_yz(0.0),
+						  quadrupole_zz(0.0),
+	  min_radius(0.0), max_radius(0.0),
+	  mean_radius(0.0),
 	  min_x(0.0), max_x(0.0),
 	  min_y(0.0), max_y(0.0),
 	  min_z(0.0), max_z(0.0),
 	  circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0),
-	  area(0.0),
-	  m_irreducible(0.0)
+	  area(0.0), irreducible_mass(0.0), areal_radius(0.0) // no comma
 { }
 
 //******************************************************************************
@@ -84,6 +96,13 @@ buffer[posn__centroid_x] = centroid_x;
 buffer[posn__centroid_y] = centroid_y;
 buffer[posn__centroid_z] = centroid_z;
 
+buffer[posn__quadrupole_xx] = quadrupole_xx;
+buffer[posn__quadrupole_xy] = quadrupole_xy;
+buffer[posn__quadrupole_xz] = quadrupole_xz;
+buffer[posn__quadrupole_yy] = quadrupole_yy;
+buffer[posn__quadrupole_xz] = quadrupole_yz;
+buffer[posn__quadrupole_zz] = quadrupole_zz;
+
 buffer[posn__min_radius]  = min_radius;
 buffer[posn__max_radius]  = max_radius;
 buffer[posn__mean_radius] = mean_radius;
@@ -98,8 +117,10 @@ buffer[posn__max_z] = max_z;
 buffer[posn__circumference_xy] = circumference_xy;
 buffer[posn__circumference_xz] = circumference_xz;
 buffer[posn__circumference_yz] = circumference_yz;
-buffer[posn__area]             = area;
-buffer[posn__m_irreducible]    = m_irreducible;
+
+buffer[posn__area]          = area;
+buffer[posn__irreducible_mass] = irreducible_mass;
+buffer[posn__areal_radius]  = areal_radius;
 }
 
 //******************************************************************************
@@ -113,6 +134,13 @@ centroid_x = buffer[posn__centroid_x];
 centroid_y = buffer[posn__centroid_y];
 centroid_z = buffer[posn__centroid_z];
 
+quadrupole_xx = buffer[posn__quadrupole_xx];
+quadrupole_xy = buffer[posn__quadrupole_xy];
+quadrupole_xz = buffer[posn__quadrupole_xz];
+quadrupole_yy = buffer[posn__quadrupole_yy];
+quadrupole_yz = buffer[posn__quadrupole_yz];
+quadrupole_zz = buffer[posn__quadrupole_zz];
+
  min_radius = buffer[posn__min_radius];
  max_radius = buffer[posn__max_radius];
 mean_radius = buffer[posn__mean_radius];
@@ -127,8 +155,10 @@ max_z = buffer[posn__max_z];
 circumference_xy = buffer[posn__circumference_xy];
 circumference_xz = buffer[posn__circumference_xz];
 circumference_yz = buffer[posn__circumference_yz];
+
             area = buffer[posn__area];
-   m_irreducible = buffer[posn__m_irreducible];
+   irreducible_mass = buffer[posn__irreducible_mass];
+    areal_radius = buffer[posn__areal_radius];
 }
 
 //******************************************************************************
@@ -137,12 +167,13 @@ circumference_yz = buffer[posn__circumference_yz];
 
 //
 // Given that an apparent horizon has been found, this function computes
-// various black hole diagnostics.
+// the black hole diagnostics.
 //
 // Inputs (gridfns)
-// h		# ghosted
-// one		# nominal
-// global_[xyz]	# nominal
+// h					# ghosted
+// one					# nominal
+// global_[xyz]				# nominal
+// global_{xx,xy,xz,yy,yz,zz}		# nominal
 //
 // Bugs:
 // The computation is rather inefficient -- we make many passes over the
@@ -165,20 +196,19 @@ max_radius = h_norms.max_abs_value();
 // xyz bounding box of horizon
 //
 
-// compute bounding box of nominal grid
-// ... this is only the stored part of the horizon if there are symmetries
+// compute bounding box of nominal grid (for stored part of the horizon only)
 jtutil::norm<fp> x_norms;
+jtutil::norm<fp> y_norms;
+jtutil::norm<fp> z_norms;
+
 ps.gridfn_norms(gfns::gfn__global_x, x_norms);
+ps.gridfn_norms(gfns::gfn__global_y, y_norms);
+ps.gridfn_norms(gfns::gfn__global_z, z_norms);
+
 min_x = x_norms.min_value();
 max_x = x_norms.max_value();
-
-jtutil::norm<fp> y_norms;
-ps.gridfn_norms(gfns::gfn__global_y, y_norms);
 min_y = y_norms.min_value();
 max_y = y_norms.max_value();
-
-jtutil::norm<fp> z_norms;
-ps.gridfn_norms(gfns::gfn__global_z, z_norms);
 min_z = z_norms.min_value();
 max_z = z_norms.max_value();
 
@@ -233,6 +263,24 @@ const fp integral_y = surface_integral(ps,
 const fp integral_z = surface_integral(ps,
 				       gfns::gfn__global_z, false, true, true,
 				       BH_diagnostics_info.integral_method);
+const fp integral_xx = surface_integral(ps,
+				       gfns::gfn__global_xx, true, true, true,
+				       BH_diagnostics_info.integral_method);
+const fp integral_xy = surface_integral(ps,
+				       gfns::gfn__global_xy, true, false, false,
+				       BH_diagnostics_info.integral_method);
+const fp integral_xz = surface_integral(ps,
+				       gfns::gfn__global_xz, false, true, false,
+				       BH_diagnostics_info.integral_method);
+const fp integral_yy = surface_integral(ps,
+				       gfns::gfn__global_yy, true, true, true,
+				       BH_diagnostics_info.integral_method);
+const fp integral_yz = surface_integral(ps,
+				       gfns::gfn__global_yz, false, false, true,
+				       BH_diagnostics_info.integral_method);
+const fp integral_zz = surface_integral(ps,
+				       gfns::gfn__global_zz, true, true, true,
+				       BH_diagnostics_info.integral_method);
 
 
 //
@@ -244,15 +292,32 @@ centroid_z = integral_z / integral_one;
 
 
 //
-// area, mean radius, and mass
+// quadrupoles (taken about centroid position)
+//
+quadrupole_xx = integral_xx / integral_one  -  centroid_x * centroid_x;
+quadrupole_xy = integral_xy / integral_one  -  centroid_x * centroid_y;
+quadrupole_xz = integral_xz / integral_one  -  centroid_x * centroid_z;
+quadrupole_yy = integral_yy / integral_one  -  centroid_y * centroid_y;
+quadrupole_yz = integral_yz / integral_one  -  centroid_y * centroid_z;
+quadrupole_zz = integral_zz / integral_one  -  centroid_z * centroid_z;
+
+
+//
+// mean radius of horizon
 //
-area = integral_one;
 mean_radius = integral_h / integral_one;
-m_irreducible = sqrt(area / (16.0*PI));
 
 
 //
-// circumferences
+// surface area and quantities derived from it
+//
+area             = integral_one;
+irreducible_mass = sqrt(area / (16.0*PI));
+areal_radius     = sqrt(area / ( 4.0*PI));
+
+
+//
+// proper circumferences
 //
 circumference_xy
   = ps.circumference("xy", gfns::gfn__h,
@@ -310,15 +375,16 @@ return ps.integrate_gridfn
 void BH_diagnostics::print(int N_horizons, int hn)
 	const
 {
+const fp irreducible_mass = sqrt(area / (16*PI));
 CCTK_VInfo(CCTK_THORNSTRING,
 	   "AH %d/%d: r=%g at (%f,%f,%f)",
 	   hn, N_horizons,
 	   double(mean_radius),
 	   double(centroid_x), double(centroid_y), double(centroid_z));
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "AH %d/%d: area=%.10g m_irreducible=%.10g",
+	   "AH %d/%d: area=%.10g irreducible_mass=%.10g",
 	   hn, N_horizons,
-	   double(area), double(m_irreducible));
+	   double(area), double(irreducible_mass));
 }
 
 //******************************************************************************
@@ -375,19 +441,35 @@ fprintf(fileptr, "# column  5 = centroid_z\n");
 fprintf(fileptr, "# column  6 = min radius\n");
 fprintf(fileptr, "# column  7 = max radius\n");
 fprintf(fileptr, "# column  8 = mean radius\n");
-fprintf(fileptr, "# column  9 = min x\n");
-fprintf(fileptr, "# column 10 = max x\n");
-fprintf(fileptr, "# column 11 = min y\n");
-fprintf(fileptr, "# column 12 = max y\n");
-fprintf(fileptr, "# column 13 = min z\n");
-fprintf(fileptr, "# column 14 = max z\n");
-fprintf(fileptr, "# column 15 = xy-plane circumference\n");
-fprintf(fileptr, "# column 16 = xz-plane circumference\n");
-fprintf(fileptr, "# column 17 = yz-plane circumference\n");
-fprintf(fileptr, "# column 18 = ratio of xz/xy-plane circumferences\n");
-fprintf(fileptr, "# column 19 = ratio of yz/xy-plane circumferences\n");
-fprintf(fileptr, "# column 20 = area\n");
-fprintf(fileptr, "# column 21 = m_irreducible\n");
+fprintf(fileptr, "# column  9 = quadrupole_xx\n");
+fprintf(fileptr, "# column 10 = quadrupole_xy\n");
+fprintf(fileptr, "# column 11 = quadrupole_xz\n");
+fprintf(fileptr, "# column 12 = quadrupole_yy\n");
+fprintf(fileptr, "# column 13 = quadrupole_yz\n");
+fprintf(fileptr, "# column 14 = quadrupole_zz\n");
+fprintf(fileptr, "# column 15 = min x\n");
+fprintf(fileptr, "# column 16 = max x\n");
+fprintf(fileptr, "# column 17 = min y\n");
+fprintf(fileptr, "# column 18 = max y\n");
+fprintf(fileptr, "# column 19 = min z\n");
+fprintf(fileptr, "# column 20 = max z\n");
+fprintf(fileptr, "# column 21 = xy-plane circumference\n");
+fprintf(fileptr, "# column 22 = xz-plane circumference\n");
+fprintf(fileptr, "# column 23 = yz-plane circumference\n");
+fprintf(fileptr, "# column 24 = ratio of xz/xy-plane circumferences\n");
+fprintf(fileptr, "# column 25 = ratio of yz/xy-plane circumferences\n");
+fprintf(fileptr, "# column 26 = area\n");
+fprintf(fileptr, "# column 27 = irreducible mass\n");
+fprintf(fileptr, "# column 28 = areal radius\n");
+fprintf(fileptr, "# column 29 = [not implemented yet] (outer) expansion Theta_(l)\n");
+fprintf(fileptr, "# column 30 = [not implemented yet] inner expansion Theta_(n)\n");
+fprintf(fileptr, "# column 31 = [not implemented yet] product of inner and outer expansions\n");
+fprintf(fileptr, "# column 32 = [not implemented yet] mean curvature\n");
+fprintf(fileptr, "# column 33 = [not implemented yet] d/d(coordinate radius) of area\n");
+fprintf(fileptr, "# column 34 = [not implemented yet] d/d(coordinate radius) of (outer) expansion Theta_(l)\n");
+fprintf(fileptr, "# column 35 = [not implemented yet] d/d(coordinate radius) of inner expansion Theta_(n)\n");
+fprintf(fileptr, "# column 36 = [not implemented yet] d/d(coordinate radius) of product of inner and outer expansions\n");
+fprintf(fileptr, "# column 37 = [not implemented yet] d/d(coordinate radius) of mean curvature\n");
 fflush(fileptr);
 
 return fileptr;
@@ -420,6 +502,14 @@ fprintf(fileptr,
 	double(min_radius), double(max_radius), double(mean_radius));
 
 fprintf(fileptr,
+     //  quadrupole_{xx,xy,xz,yy,yz,zz}
+     //  ======  ======  ======  ======  ======  ======
+	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t",
+	double(quadrupole_xx), double(quadrupole_xy), double(quadrupole_xz),
+			       double(quadrupole_yy), double(quadrupole_yz),
+						      double(quadrupole_zz));
+
+fprintf(fileptr,
      //  {min,max}_x     {min,max}_y     {min,max}_z
      //  ==============  ==============  ==============
 	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t",
@@ -428,17 +518,41 @@ fprintf(fileptr,
 	double(min_z), double(max_z));
 
 fprintf(fileptr,
-     //  {xy,xz,yz}-plane         xz/xy  yz/xy   area    m_irreducible
-     //  circumferences          circumference   ======  ======
+     //  {xy,xz,yz}-plane         xz/xy  yz/xy
+     //  circumferences          circumference
      //                          ratios
-     //  ======================  ==============  ======  ======
-	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
+     //  ======================  ==============
+	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t",
 	double(circumference_xy),
 	double(circumference_xz),
 	double(circumference_yz),
 	double(circumference_xz / circumference_xy),
-	double(circumference_yz / circumference_xy),
-	double(area), double(m_irreducible));
+	double(circumference_yz / circumference_xy));
+
+fprintf(fileptr,
+     //  area            areal_radius
+     //  ======  irreducible_mass
+     //  ======  ======  ======
+	"%#.10g\t%#.10g\t%#.10g\t",
+	double(area), double(irreducible_mass), double(areal_radius));
+
+// these diagnostics aren't implemented yet :(
+fprintf(fileptr,
+     //  expansion       product_expansion
+     //  ======  inner_expansion
+     //  ======  ======  ======  mean_curvature
+     //  ======  ======  ======  ======
+	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t",
+	0.0, 0.0, 0.0, 0.0);
+
+// these diagnostics aren't implemented yet :(
+fprintf(fileptr,
+     //  dr_area         dr_inner_expansion
+     //  ======  dr_expansion    dr_product_expansion
+     //  ======  ======  ======  ======  dr_mean_curvature
+     //  ======  ======  ======  ======  ======
+	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
+	0.0, 0.0, 0.0, 0.0, 0.0);
 
 fflush(fileptr);
 }
diff --git a/src/driver/BH_diagnostics.hh b/src/driver/BH_diagnostics.hh
index 5e8dc37..0d368e8 100644
--- a/src/driver/BH_diagnostics.hh
+++ b/src/driver/BH_diagnostics.hh
@@ -4,6 +4,7 @@
 //
 // prerequisites:
 //	<stdio.h>
+//	"cctk.h"
 //
 
 // everything in this file is inside this namespace
@@ -30,33 +31,58 @@ struct	BH_diagnostics_info
 // Note that all the diagnostics are for the full apparent horizon, even
 // if we don't actually store all of it due to patch-system symmetries.
 //
+// All the "mean" quantities are computed via surface integrals using the
+// induced metric on the surface.
+//
 struct	BH_diagnostics
 	{
 public:
+	// mean x,y,z
 	fp centroid_x, centroid_y, centroid_z;
+
+	// these are quadrupole moments about the centroid, i.e.
+	// mean(xi*xj) - centroid_i*centroid_j
+	fp quadrupole_xx, quadrupole_xy, quadrupole_xz,
+			  quadrupole_yy, quadrupole_yz,
+					 quadrupole_zz;
+
+	// min,max,mean surface radius about local coordinate origin
 	fp min_radius, max_radius, mean_radius;
 
 	// xyz bounding box
-	fp min_x, max_x, min_y, max_y, min_z, max_z;
+	fp min_x, max_x,
+	   min_y, max_y,
+	   min_z, max_z;
 
-	fp circumference_xy, circumference_xz, circumference_yz;
-	fp area;
-	fp m_irreducible;
+	// proper circumference
+	// (computed using induced metric along these local-coordinate planes)
+	fp circumference_xy,
+	   circumference_xz,
+	   circumference_yz;
+
+	// surface area (computed using induced metric)
+	// and quantities derived from it
+	fp area, irreducible_mass, areal_radius;
 
 public:
 	// position of diagnostics in buffer and number of diagnostics
 	enum	{
 		posn__centroid_x = 0, posn__centroid_y, posn__centroid_z,
+		posn__quadrupole_xx, posn__quadrupole_xy, posn__quadrupole_xz,
+				     posn__quadrupole_yy, posn__quadrupole_yz,
+							  posn__quadrupole_zz,
 		posn__min_radius, posn__max_radius, posn__mean_radius,
 
 		posn__min_x, posn__max_x,
 		posn__min_y, posn__max_y,
 		posn__min_z, posn__max_z,
 
-		posn__circumference_xy, posn__circumference_xz,
-					posn__circumference_yz,
-		posn__area,
-		posn__m_irreducible,
+		posn__circumference_xy,
+		posn__circumference_xz,
+		posn__circumference_yz,
+
+		posn__area, posn__irreducible_mass, posn__areal_radius,
+
 		N_buffer // no comma	// size of buffer
 		};
 
diff --git a/src/driver/README b/src/driver/README
index 66712cd..f3de3e7 100644
--- a/src/driver/README
+++ b/src/driver/README
@@ -25,7 +25,11 @@ mask.cc			# sees CCTK_ARGUMENTS
 
 announce.cc		# sees CCTK_ARGUMENTS
 	this is called from the scheduler to announce apparent horizon
-	info to other thorns
+	info to other thorns via aliased function calls
+
+spherical_surface.cc	# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
+	this is called from the scheduler to store apparent horizon
+	info in the SphericalSurface variables
 
 aliased_functions.cc	# sees CCTK_ARGUMENTS, CCTK_FUNCTIONS
 	this is called from other thorns via the flesh aliased-function
diff --git a/src/driver/announce.cc b/src/driver/announce.cc
index f604d1c..8116c6a 100644
--- a/src/driver/announce.cc
+++ b/src/driver/announce.cc
@@ -84,11 +84,14 @@ if (state.AH_data_array[hn] == NULL)
 // is there anyone to announce it to?
 if (CCTK_IsFunctionAliased("SetDriftCorrectPosition"))
    then {
+	assert(state.AH_data_array[hn] != NULL);
 	const struct AH_data& AH_data = *state.AH_data_array[hn];
+
 	const struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
 	const CCTK_REAL xx = BH_diagnostics.centroid_x;
 	const CCTK_REAL yy = BH_diagnostics.centroid_y;
 	const CCTK_REAL zz = BH_diagnostics.centroid_z;
+
 	if (verbose_info.print_physics_details)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "horizon %d centroid (%g,%g,%g) --> DriftCorrect",
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index 65cb95b..1233b50 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -309,6 +309,13 @@ struct	AH_data
 	struct BH_diagnostics BH_diagnostics;
 	FILE *BH_diagnostics_fileptr;
 
+	bool store_info_in_SS_vars;	// should we store this horizon
+					// in the SphericalSurface variables?
+	int SS_surface_number;		// if  store_info_in_SS_vars ,
+					//   this is the SphericalSurface
+					//   surface number to store into;
+					// otherwise this is ignored
+
 	// interprocessor-communication buffers
 	// for this horizon's BH diagnostics and (optionally) horizon shape
 	struct horizon_buffers horizon_buffers;
@@ -395,6 +402,11 @@ extern "C"
 extern "C"
   void AHFinderDirect_announce(CCTK_ARGUMENTS);
 
+// spherical_surface.cc
+// ... called from Cactus Scheduler
+extern "C"
+  void AHFinderDirect_store_SS_info(CCTK_ARGUMENTS);
+
 // mask.cc
 // ... called from Cactus Scheduler
 extern "C"
diff --git a/src/driver/make.code.defn b/src/driver/make.code.defn
index 737cb3d..9cbbddd 100644
--- a/src/driver/make.code.defn
+++ b/src/driver/make.code.defn
@@ -7,7 +7,7 @@ SRCS = state.cc \
        initial_guess.cc Newton.cc \
        io.cc misc-driver.cc \
        BH_diagnostics.cc horizon_sequence.cc \
-       mask.cc announce.cc aliased_functions.cc
+       mask.cc announce.cc spherical_surface.cc aliased_functions.cc
 
 # Subdirectories containing source files
 SUBDIRS =
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 3ec6ca1..13a404f 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -570,6 +570,25 @@ if (strlen(surface_interpolator_name) > 0)
 	AH_data.found_flag = false;
 	AH_data.h_files_written = false;
 	AH_data.BH_diagnostics_fileptr = NULL;
+
+	AH_data.store_info_in_SS_vars = (which_surface_to_store_info[hn] >= 0);
+	if (AH_data.store_info_in_SS_vars)
+	   then {
+		AH_data.SS_surface_number = which_surface_to_store_info[hn];
+		if (AH_data.SS_surface_number
+		    >= /* SphericalSurface:: */ nsurfaces)
+		   then CCTK_VWarn(FATAL_ERROR,
+				   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   AHFinderDirect_setup():\n"
+"        which_surface_to_store_info[%d] = %d\n"
+"        is outside the legal range [0,SphericalSurface::nsurfaces=%d)!\n"
+				   ,
+				   hn, AH_data.SS_surface_number,
+				   int(/* SphericalSurface:: */ nsurfaces));
+								/*NOTREACHED*/
+		}
+	   else AH_data.SS_surface_number = 0;	// dummy value
 	}
 	  }
 }
diff --git a/src/driver/spherical_surface.cc b/src/driver/spherical_surface.cc
new file mode 100644
index 0000000..20617bd
--- /dev/null
+++ b/src/driver/spherical_surface.cc
@@ -0,0 +1,257 @@
+// spherical_surface.cc -- interface with SphericalSurface thorn
+// $Header$
+//
+// <<<access to persistent data>>>
+// AHFinderDirect_store_SS_info - ... SphericalSurface information
+/// store_diagnostic_info - store BH_diagnostics info in SphericalSurface vars
+/// store_horizon_shape - store horizon shape in SphericalSurface vars
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "config.h"
+#include "stdc.h"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
+#include "driver.hh"
+
+// all the code in this file is inside this namespace
+namespace AHFinderDirect
+	  {
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+void store_diagnostic_info(CCTK_ARGUMENTS,
+			   const patch_system& ps,
+			   const struct BH_diagnostics& BH_diagnostics,
+			   const int sn);
+void store_horizon_shape(CCTK_ARGUMENTS,
+			 const patch_system& ps,
+			 const int sn);
+	  }
+
+//******************************************************************************
+
+//
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler, to store any desired
+// apparent horizon info to the SphericalSurface variables.
+//
+// Cactus parameters used:
+// SphericalSurface::nsurfaces = (in)
+// 
+extern "C"
+  void AHFinderDirect_store_SS_info(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+const int N_surfaces = /* SphericalSurface:: */ nsurfaces;
+const struct verbose_info& verbose_info = state.verbose_info;
+
+	for (int hn = 1; hn <= N_horizons; ++ hn)
+	{
+	assert(state.AH_data_array[hn] != NULL);
+	const struct AH_data& AH_data = *state.AH_data_array[hn];
+	assert(AH_data.ps_ptr != NULL);
+	const patch_system& ps = *AH_data.ps_ptr;
+
+	if (! AH_data.store_info_in_SS_vars)
+	   then continue;				// *** LOOP CONTROL ***
+	const int sn = AH_data.SS_surface_number;
+	assert(sn >= 0);
+	assert(sn < N_surfaces);
+
+	if (! state.find_now(cctk_iteration))
+	   then {
+		// we didn't try to find this (or any!) horizon
+		// at this time step
+		/* SphericalSurface:: */ sf_valid[sn] = 0;
+		continue;				// *** LOOP CONTROL ***
+		}
+
+	if (! AH_data.found_flag)
+	   then {
+		// we tried to find this horizon, but failed
+		/* SphericalSurface:: */ sf_valid[sn] = -1;
+		continue;				// *** LOOP CONTROL ***
+		}
+
+	// get to here ==> we found this horizon
+	/* SphericalSurface:: */ sf_valid[sn] = 1;
+
+	if (verbose_info.print_algorithm_highlights)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "storing horizon %d info in SphericalSurface surface %d",
+			   hn, sn);
+
+	const struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
+	store_diagnostic_info(CCTK_PASS_CTOC, ps, BH_diagnostics, sn);
+
+	store_horizon_shape(CCTK_PASS_CTOC, ps, sn);
+	}
+}
+
+//******************************************************************************
+
+//
+// Assuming that we found this horizon, this function stores various
+// diagnostic info about it in the corresponding SphericalSurface variables.
+//
+// Arguments:
+// sn = (in) The SphericalSurface surface number to store the information in.
+//
+// Cactus variables:
+// sf_* = (out) The SphericalSurface variables.
+//
+namespace {
+void store_diagnostic_info(CCTK_ARGUMENTS,
+			   const patch_system& ps,
+			   const struct BH_diagnostics& BH_diagnostics,
+			   const int sn)
+{
+DECLARE_CCTK_ARGUMENTS
+
+/* SphericalSurface:: */ sf_origin_x[sn] = ps.origin_x();
+/* SphericalSurface:: */ sf_origin_y[sn] = ps.origin_y();
+/* SphericalSurface:: */ sf_origin_z[sn] = ps.origin_z();
+
+/* SphericalSurface:: */ sf_mean_radius  [sn] = BH_diagnostics.areal_radius;
+/* SphericalSurface:: */ sf_min_radius   [sn] = BH_diagnostics.min_radius;
+/* SphericalSurface:: */ sf_max_radius   [sn] = BH_diagnostics.max_radius;
+
+/* SphericalSurface:: */ sf_centroid_x   [sn] = BH_diagnostics.centroid_x;
+/* SphericalSurface:: */ sf_centroid_y   [sn] = BH_diagnostics.centroid_y;
+/* SphericalSurface:: */ sf_centroid_z   [sn] = BH_diagnostics.centroid_z;
+/* SphericalSurface:: */ sf_quadrupole_xx[sn] = BH_diagnostics.quadrupole_xx;
+/* SphericalSurface:: */ sf_quadrupole_xy[sn] = BH_diagnostics.quadrupole_xy;
+/* SphericalSurface:: */ sf_quadrupole_xz[sn] = BH_diagnostics.quadrupole_xz;
+/* SphericalSurface:: */ sf_quadrupole_yy[sn] = BH_diagnostics.quadrupole_yy;
+/* SphericalSurface:: */ sf_quadrupole_yz[sn] = BH_diagnostics.quadrupole_yz;
+/* SphericalSurface:: */ sf_quadrupole_zz[sn] = BH_diagnostics.quadrupole_zz;
+
+/* SphericalSurface:: */ sf_min_x        [sn] = BH_diagnostics.min_x;
+/* SphericalSurface:: */ sf_max_x        [sn] = BH_diagnostics.max_x;
+/* SphericalSurface:: */ sf_min_y        [sn] = BH_diagnostics.min_y;
+/* SphericalSurface:: */ sf_max_y        [sn] = BH_diagnostics.max_y;
+/* SphericalSurface:: */ sf_min_z        [sn] = BH_diagnostics.min_z;
+/* SphericalSurface:: */ sf_max_z        [sn] = BH_diagnostics.max_z;
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function stores our information about a specified horizon to the
+// SphericalSurface variables.
+//
+// Arguments:
+// sn = (in) The SphericalSurface surface number to store the information in.
+//
+// Cactus variables:
+// sf_maxn{theta,phi} = (in) The SphericalSurface array dimensions for
+//			     the 2-D array indexing.
+// sf_n{theta,phi} = (in) The SphericalSurface array dimensions for the
+//			  part of the 2-D array that's actually used.
+// sf_radius = (out) The SphericalSurface radius.
+//
+// FIXME: The present implementation is quite inefficient, as it calls
+//        patch_system::radius_in_local_xyz_direction() (which does a
+//	  separate interpolator call) for each point.  It would be more
+//	  efficient to have a new  patch_system::  API which operated
+//	  on a whole array of points at once, to amortize the interpolator
+//	  overhead.
+//
+namespace {
+void store_horizon_shape(CCTK_ARGUMENTS,
+			 const patch_system& ps,
+			 const int sn)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+const double origin_theta = /* SphericalSurface:: */ sf_origin_theta[sn];
+const double origin_phi   = /* SphericalSurface:: */ sf_origin_phi  [sn];
+const double delta_theta  = /* SphericalSurface:: */ sf_delta_theta [sn];
+const double delta_phi    = /* SphericalSurface:: */ sf_delta_phi   [sn];
+
+const int     N_theta  = /* SphericalSurface:: */    ntheta[sn];
+const int     N_phi    = /* SphericalSurface:: */    nphi  [sn];
+const int max_N_theta  = /* SphericalSurface:: */ maxntheta;
+const int max_N_phi    = /* SphericalSurface:: */ maxnphi;
+
+// we want Fortran loop nesting here for cache efficiency in storing
+// to the SphericalSurface::sf_radius array (see comment below)
+	for (int i_phi = 0 ; i_phi < N_phi ; ++i_phi)
+	{
+	const double phi = origin_phi + i_phi*delta_phi;
+	const double sin_phi = sin(phi);
+	const double cos_phi = cos(phi);
+
+		for (int i_theta = 0 ; i_theta < N_theta ; ++i_theta)
+		{
+		const double theta = origin_theta + i_theta*delta_theta;
+		const double sin_theta = sin(theta);
+		const double cos_theta = cos(theta);
+
+		const double local_x = sin_theta * cos_phi;
+		const double local_y = sin_theta * sin_phi;
+		const double local_z = cos_theta;
+
+		const double r = ps.radius_in_local_xyz_direction
+					(gfns::gfn__h,
+					 local_x, local_y, local_z);
+
+		// SphericalSurface::sf_radius is actually stored as
+		// a 3-D contiguous array, with indices
+		//    theta     (contiguous, i.e. stride=1)
+		//    phi
+		//    surface   (largest stride)
+		const int sub = i_theta + max_N_theta * (i_phi + max_N_phi*sn);
+		/* SphericalSurface:: */ sf_radius[sub] = r;
+		}
+	}
+}
+	  }
+
+//******************************************************************************
+
+	  }	// namespace AHFinderDirect
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index e0b7ba8..4d7ae34 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -256,8 +256,17 @@ return expansion_success;				// *** NORMAL RETURN ***
 //******************************************************************************
 
 //
-// This function sets up the global xyz positions of the grid points
-// in the gridfns global_[xyz].  These will be used by interplate_geometry().
+// This function sets up the xyz-position gridfns:
+// * global_[xyz] (used by  interplate_geometry() )
+// * global_{xx,xy,xz,yy,yz,zz} (used by higher-level code to compute
+//                               quadrupole moments of horizons)
+//
+// Bugs:
+// * We initialize the global_{xx,xy,xz,yy,yz,zz} gridfns every time
+//   this function is called, i.e. at each horizon-finding Newton
+//   iteration, even though they're only needed at the end of the
+//   horizon-finding process.  In practice the extra cost is small,
+//   though, so it's probably not worth fixing this...
 //
 namespace {
 void setup_xyz_posns(patch_system& ps, bool print_msg_flag)
@@ -290,6 +299,20 @@ if (print_msg_flag)
 		p.gridfn(gfns::gfn__global_x, irho,isigma) = global_x;
 		p.gridfn(gfns::gfn__global_y, irho,isigma) = global_y;
 		p.gridfn(gfns::gfn__global_z, irho,isigma) = global_z;
+
+		const fp global_xx = global_x * global_x;
+		const fp global_xy = global_x * global_y;
+		const fp global_xz = global_x * global_z;
+		const fp global_yy = global_y * global_y;
+		const fp global_yz = global_y * global_z;
+		const fp global_zz = global_z * global_z;
+
+		p.gridfn(gfns::gfn__global_xx, irho,isigma) = global_xx;
+		p.gridfn(gfns::gfn__global_xy, irho,isigma) = global_xy;
+		p.gridfn(gfns::gfn__global_xz, irho,isigma) = global_xz;
+		p.gridfn(gfns::gfn__global_yy, irho,isigma) = global_yy;
+		p.gridfn(gfns::gfn__global_yz, irho,isigma) = global_yz;
+		p.gridfn(gfns::gfn__global_zz, irho,isigma) = global_zz;
 		}
 		}
 	}
diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh
index d7a0afd..ee80c90 100644
--- a/src/gr/gfns.hh
+++ b/src/gr/gfns.hh
@@ -39,6 +39,13 @@ enum	{
 	gfn__global_y,				// no access macro
 	gfn__global_z,				// no access macro
 
+	gfn__global_xx,				// no access macro
+	gfn__global_xy,				// no access macro
+	gfn__global_xz,				// no access macro
+	gfn__global_yy,				// no access macro
+	gfn__global_yz,				// no access macro
+	gfn__global_zz,				// no access macro
+
 	gfn__g_dd_11,
 	gfn__g_dd_12,
 	gfn__g_dd_13,