* 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

32 files changed, 709 insertions, 111 deletions

diff --git a/README b/README
index 7207e71..dac5714 100644
--- a/README
+++ b/README
@@ -120,23 +120,37 @@ Timothy A. Davis, and is subect to the UMFPACK License:
 Other Software Required
 =======================
 
-This thorn inherits from ADMBase, StaticConformal, SpaceMask, and IO
-(all in the CactusEinstein arrangement).
+This thorn inherits from a number of other thorns (strictly speaking,
+implementations):
+
+Implementation		Typically provided by Thorn
+--------------		---------------------------
+Grid			CactusBase/CartGrid3d
+IO			CactusBase/IOUtil
+ADMBase			CactusEinstein/ADMBase
+StaticConformal		CactusEinstein/StaticConformal
+SpaceMask		CactusEinstein/SpaceMask
+SphericalSurface	AEIThorns/SphericalSurface
+
+You can get CactusBase and CactusEinstein thorns from the same place
+you got Cactus itself.  You can get AEIThorns thorns by anonymous CVS
+checkout:
+  cvs -d cvs_anon@cvs.aei.mpg.de:/numrelcvs checkout AEIThorns/SphericalSurface
 
 By default this thorn uses various Cactus APIs which are supplied by
 other thorns:
 * This thorn uses CCTK_InterpGridArrays() for interpolating grid arrays.
-  This is supplied by a driver-specific thorn, such as PUGHInterp or
-  CarpetInterp.
+  This is supplied by a driver-specific global interpolation thorn,
+  such as PUGHInterp or CarpetInterp.
 * This thorn uses CCTK_ReduceLocArrayToArray1D() for interprocessor
   communication in the multiprocessor Newton solver.
 	[see src/driver/README.parallel for details]
-  This is supplied by a driver-specific thorn; at present PUGHReduce
-  it this for the PUGH driver.
-* This thorn uses CCTK_InterpLocalUniform() for interpatch and surface
-  interpolation, and it needs various interpolation options which at
-  present are (only) supported by the AEIThorns/AEILocalInterp local
-  interpolator.
+  This is supplied by a driver-specific thorn global reduction thorn,
+  such as PUGHReduce or CarpetReduce.
+* This thorn uses CCTK_InterpLocalUniform() for (processor-local)
+  interpatch and surface interpolation, and it needs various interpolation
+  options which at present are (only) supported by the
+  AEIThorns/AEILocalInterp local interpolator.
 
 This thorn is written in C++, so you'll need a C++ compiler -- in fact
 a fairly modern one -- to compile this thorn.  See the "Compilation Notes"
@@ -174,7 +188,7 @@ In particular:
   form.)
 * <assert.h> is used fairly heavily for sanity checks.
 * To avoid various portability problems, none of the C++ standard
-  template library (STL) is used.
+  template library (STL) is used. :(
 
 
 Compiler Notes
diff --git a/doc/documentation.tex b/doc/documentation.tex
index 5c3709a..616545d 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -106,7 +106,7 @@
 \def\defn#1{``#1''}			% definition of a term
 \def\arrangement#1{{\bf #1}}		% name of an arrangement
 \def\thorn#1{{\bf #1}}			% name of a thorn
-\def\cvsplace#1{{\bf #1}}		% name of a CVS repository/directory
+\def\cvsplace#1{{\bf #1}}		% name of a CVS repository/directory/tag
 \def\cf{\hbox{cf.\hbox{}}}
 \def\eg{\hbox{eg.\hbox{}}}
 \def\ie{\hbox{i.e.\hbox{}}}
@@ -254,13 +254,42 @@ here's what works and what doesn't:
 \item	\thorn{AHFinderDirect} works fine with the \thorn{PUGH} unigrid
 	driver and with the \thorn{Carpet} mesh-refinement driver.  So
 	far as I know it's never been tested with any other driver.
+\item	\thorn{AHFinderDirect} works fine in single- or multi-processor
+	Cactus runs.
 \item	Obviously, your Cactus configuration must include
 	\thorn{AHFinderDirect}, and your \code{ActiveThorns} parameter(s)
-	must activate it.  \thorn{AHFinderDirect} inherits from
-	\thorn{Grid}, \thorn{ADMBase}, \thorn{StaticConformal}, 
-	\thorn{SpaceMask}, and \thorn{IO}, so you'll need those
-	(or more precisely some thorns providing them),
+	must activate it.
+\item	\thorn{AHFinderDirect} inherits from the other thorns
+	(strictly speaking, implementations) listed in
+	table~\ref{AEIThorns/AHFinderDirect/tab-inherits-from},
+	so you'll need them (or more precisely some thorns providing them)
 	in your configuration and activated, too.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{table}[tbp]
+\begin{center}
+\begin{tabular}{l@{\qquad}l}
+\thorn{Implementation}		& Typically provided by \thorn{Thorn}	\\
+\hline %---------------------------------------------------------------
+\thorn{Grid}			& \thorn{CactusBase/CartGrid3d}		\\
+\thorn{IO}			& \thorn{CactusBase/IOUtil}		\\
+\thorn{ADMBase}			& \thorn{CactusEinstein/ADMBase}	\\
+\thorn{StaticConformal}		& \thorn{CactusEinstein/StaticConformal}\\
+\thorn{SpaceMask}		& \thorn{CactusEinstein/SpaceMask}	\\
+\thorn{SphericalSurface}	& \thorn{AEIThorns/SphericalSurface}	%%%\\
+\end{tabular}
+\end{center}
+\caption[Other Thorns from which \thorn{AHFinderDirect} Inherits]
+	{
+	This table lists all the other implementations from which
+	\thorn{AHFinderDirect} inherits, and the thorns which typically
+	provide these implementations.
+	}
+\label{AEIThorns/AHFinderDirect/tab-inherits-from}
+\end{table}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 \item	\verb|Grid::domain = "full"|, \verb|"bitant"|,
 	\verb|"quadrant"|, and \verb|"octant"| are supported.
 	Alas, at present rotating (or more precisely nonlocal)
@@ -280,8 +309,6 @@ here's what works and what doesn't:
 	in this API which at present are only supported by thorn
 	\thorn{AEILocalInterp} (so you must have this thorn compiled in
 	and activated)
-\item	\thorn{AHFinderDirect} works fine in single- or multi-processor
-	Cactus runs.
 \item	\thorn{AHFinderDirect} uses various Cactus reduction APIs to
 	coordinate multi-processor horizon finding, so (even if you're
 	only going to run on a single processor) you must have a reduction
@@ -297,18 +324,12 @@ here's what works and what doesn't:
 	files.  This is a bug, which I plan to fix soon (the right behavior
 	is/will be to {\em append\/} to the existing \verb|BH_diagnostics|
 	file).
-\item	\thorn{AHFinderDirect} can set an excision mask based on
-	each horizon's shape, using either the old-style (\verb|CCTK_REAL|)
-	mask (compatible with \thorn{AHFinder}) or the new-style
-	(\verb|CCTK_INT|) mask bit-fields defined by \thorn{SpaceMask}
-	(or even both simultaneously).
-\item	\thorn{AHFinderDirect} can notify \thorn{DriftCorrect} about
-	a selected horizon's centroid position, for use in adjusting
-	corotating shift vectors.  This uses the new function-aliasing
-	version of \thorn{DriftCorrect}, not the old version which
-	worked with an auxiliary thorn \thorn{AHFSetDCCentroid}.
 \end{itemize}
 
+\thorn{AHFinderDirect} can pass information to the rest of Cactus in
+several ways; these are described in detail in
+section~\ref{AHFinderDirect/sect-parameters/communicating-with-other-thorns}.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Obtaining and Compiling \thorn{AHFinderDirect}}
@@ -777,7 +798,8 @@ to data files:
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\subsection{Other I/O Parameters}
+\subsection{Parameters for the ``BH\_diagnostics'' Files}
+\label{AHFinderDirect/sect-parameters/BH-diagnostics-parameters}
 
 As well as the apparent horizon shape files, this thorn can also
 write files giving time series of various diagnostics.  These are
@@ -845,6 +867,8 @@ controlled by the following parameters:
 		\item	the proper area of the apparent horizon, $A$
 		\item	the irreducible mass of the apparent horizon,
 			$\sqrt{A/16\pi}$
+		\item	the areal radius of the apparent horizon,
+			$\sqrt{A/4\pi}$
 		\end{itemize}
 	\end{itemize}
 
@@ -1148,28 +1172,40 @@ to do this before \thorn{AHFinderDirect} tries to find the horizon(s).
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{Communicating with Other Thorns}
+\label{AHFinderDirect/sect-parameters/communicating-with-other-thorns}
 
 Besides the data files it writes, \thorn{AHFinderDirect} currently
 has three ways to communicate with other Cactus thorns:
 \begin{itemize}
-\item	It can set a mask grid function(s), which can be used
-	for excision or other purposes.  This is described in
-	section~\ref{AHFinderDirect/sect-parameters/mask-parameters}.
-\item	It can announce the horizon centroid to a some other
-	thorn (typicaly \thorn{DriftCorrect}) using a
-	function-aliasing mechanism; this is described below;
-	this is described in
-	section~\ref{AHFinderDirect/sect-parameters/communicating-with-other-thorns/announcing-centroid}.
-\item	It provides a set of aliased functions which any other
-	thorn(s) can call to find out the shape of a specified horizon;
-	this is described in
-	section~\ref{AHFinderDirect/sect-parameters/communicating-with-other-thorns/horizon-shape-functions}.
+\item	\thorn{AHFinderDirect} can set a mask grid function(s) based
+	on the (a) horizon's shape; other thorns may then use this for
+	excision or other purposes.  \thorn{AHFinderDirect} supports
+	both the old-style (\verb|CCTK_REAL|) mask
+	(compatible with \thorn{AHFinder}) or
+	the new-style (\verb|CCTK_INT|) mask bit-fields defined by
+	\thorn{SpaceMask}; you can even use both styles simultaneously.
+\item	\thorn{AHFinderDirect} can announce a selected horizon's centroid
+	position to another thorn (typically \thorn{DriftCorrect}, which
+	uses this to adjust its corotating shift vector).
+	This uses the new function-aliasing version of \thorn{DriftCorrect},
+	not the old version which worked with an auxiliary thorn
+	\thorn{AHFSetDCCentroid}.
+\item	\thorn{AHFinderDirect} provides a set of aliased functions
+	which any other thorn(s) can call to find out the shape of
+	a specified horizon.
+\item	\thorn{AHFinderDirect} can store information about the
+	horizon(s) it finds in the \thorn{SphericalSurface} variables
+	for other thorns to use.
 \end{itemize}
 
+\thorn{AHFinderDirect}'s mask features are described in
+section~\ref{AHFinderDirect/sect-parameters/mask-parameters};
+the other communication mechanisms are described in the following
+subsections.
+
 %%%%%%%%%%%%%%%%%%%%
 
 \subsubsection{Parameters for Announcing a Horizon Centroid to Other Thorns}
-\label{AHFinderDirect/sect-parameters/communicating-with-other-thorns/announcing-centroid}
 
 This thorn can optionally announce the centroid of a specified
 apparent horizon to another thorn (typically \thorn{DriftCorrect})
@@ -1188,7 +1224,6 @@ following parameter:
 %%%%%%%%%%%%%%%%%%%%
 
 \subsubsection{Aliased Functions to Provide Horizon-Shape Information}
-\label{AHFinderDirect/sect-parameters/communicating-with-other-thorns/horizon-shape-functions}
 
 \thorn{AHFinderDirect} provides the following aliased functions
 to allow other thorns to find out about the horizons.  Each function
@@ -1223,6 +1258,29 @@ CCTK_INT FUNCTION HorizonRadiusInDirection
     CCTK_REAL OUT ARRAY radius)
 \end{verbatim}
 
+%%%%%%%%%%%%%%%%%%%%
+
+\subsubsection{Storing Horizon-Shape Information in the \thorn{SphericalSurface} Variables}
+
+\thorn{SphericalSurface} (in the \thorn{AEIThorns} arrangement)
+defines a set of generic grid arrays which describe ``spherical surfaces''.
+\thorn{AHFinderDirect} can optionally store information about the
+horizons it finds in the \thorn{SphericalSurface} variables.
+This is controlled by the following parameters:
+\begin{description}
+\item[\code{which\_surface\_to\_store\_info[}$n$\code{]}]
+\mbox{}\\
+	This parameter should be set to the \thorn{SphericalSurface}
+	surface number into which information on a given \thorn{AHFinderDirect}
+	horizon should be stored, or to -1 to skip storing the information.
+	It defaults to -1 for each horizon.
+
+	At present, if multiple \thorn{AHFinderDirect} horizons
+	specify the same \thorn{SphericalSurface} surface, the
+	highest-numbered horizon will ``win'', \ie{} it will overwrite
+	the data from any lower-numbered horizons.
+\end{description}
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{Other Parameters}
@@ -1344,6 +1402,7 @@ CCTK_INT FUNCTION HorizonRadiusInDirection
 	in an Eddington-Finkelsteon slice of the unit-mass
 	Schwarzschild spacetime.
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{figure}[htbp]
 \begin{center}
 \includegraphics{AEIThorns_AHFinderDirect_Schw_EF_Theta_of_r.eps}
@@ -1358,6 +1417,7 @@ CCTK_INT FUNCTION HorizonRadiusInDirection
 	}
 \label{AHFinderDirect/fig-Schwarzschild-EF-Theta(r)}
 \end{figure}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \end{description}
 
@@ -1370,7 +1430,7 @@ is doing during a Cactus run: the \verb|BH_diagnostics| files and the
 \verb|CCTK_INFO| messages written to the Cactus standard output:
 
 The \verb|BH_diagnostics| files are described in detail in
-section~\ref{AHFinderDirect/sect-parameters/other-parameters}.
+section~\ref{AHFinderDirect/sect-parameters/BH-diagnostics-parameters}.
 These files are written and ``flushed'' at each time step, so they're
 always up-to-date.
 
@@ -1682,9 +1742,19 @@ interested in some other related thorns:
 \item[\thorn{TGRapparentHorizon2D}] (in the \arrangement{TAT} arrangement)
 	was written by Erik Schnetter, and is another apparent horizon
 	finder.  It uses methods very similar to this thorn, and (like
-	this thorn) is very fast and accurate.  It's described in detail
-	in the papers~\cite{AEIThorns/AHFinderDirect/Schnetter02a}
+	this thorn) is very fast and accurate.  However, it's no longer
+	under active development.  It's described in detail in the
+	papers~\cite{AEIThorns/AHFinderDirect/Schnetter02a}
 	and~\cite{AEIThorns/AHFinderDirect/Schnetter03a}.
+\item[\thorn{AHFinderDirect} (\cvsplace{Erik} branch)]
+      (in the \arrangement{AEIThorns} arrangement)\\
+	Erik Schnetter has added a number of new features to
+	\thorn{AHFinderDirect} on a CVS branch with the tag \cvsplace{Erik},
+	including horizon pretracking (to locate places where horizons
+	are about to form), and the ability to find constant-expansion
+	and constant-mean-curvature surfaces specified by their areal radius.
+	We hope to integrate these into the main \thorn{AHFinderDirect}
+	branch during the summer of 2004.
 \end{description}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -1692,8 +1762,13 @@ interested in some other related thorns:
 \section{Acknowledgments}
 
 I thank Peter Diener, Ian Hawke, and Erik Schnetter for many valuable
-conversations, to Thomas Radke for his work on the new interpolators,
-and the whole Cactus crew for a great infrastructure!
+conversations.  I think Thomas Radke for his work on the new interpolators.
+I thank the whole Cactus crew for a great infrastructure!
+
+Erik Schnetter originally implemented a number of improvements to
+this thorn, notably the \thorn{SphericalSurface} interface and the
+new features in the \cvsplace{Erik} branch.
+
 I thank the Alexander von Humboldt foundation and the AEI visitors'
 and postdoctoral fellowships programs for financial support.
 
diff --git a/interface.ccl b/interface.ccl
index 41900df..9ca619d 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,7 +2,8 @@
 # $Header$
 
 implements: AHFinderDirect
-inherits: Grid ADMBase StaticConformal SpaceMask IO
+# thorns from: CactusBase   CactusEinstein .................    AEIThorns
+inherits:      Grid IO      ADMBase StaticConformal SpaceMask   SphericalSurface
 
 # we use an include file provided by SpaceMask
 USES INCLUDE: SpaceMask.h
diff --git a/par/Kerr-centered.par b/par/Kerr-centered.par
index a327861..c2b21cf 100644
--- a/par/Kerr-centered.par
+++ b/par/Kerr-centered.par
@@ -37,6 +37,7 @@ IOUtil::parfile_write = "no"
 
 ########################################
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/Kerr-mask.par b/par/Kerr-mask.par
index f222b7a..a69d534 100644
--- a/par/Kerr-mask.par
+++ b/par/Kerr-mask.par
@@ -42,6 +42,7 @@ IOASCII::out2D_vars = "SpaceMask::emask"
 
 ########################################
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/Kerr-order3.par b/par/Kerr-order3.par
index 705761c..36610aa 100644
--- a/par/Kerr-order3.par
+++ b/par/Kerr-order3.par
@@ -38,6 +38,7 @@ IOUtil::parfile_write = "no"
 
 ########################################
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/Kerr-tiny.par b/par/Kerr-tiny.par
index f0ce3a0..0da4804 100644
--- a/par/Kerr-tiny.par
+++ b/par/Kerr-tiny.par
@@ -40,6 +40,7 @@ IOUtil::parfile_write = "no"
 
 ########################################
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 AHFinderDirect::h_base_file_name     = "Kerr-tiny.h"
 
diff --git a/par/Kerr.par b/par/Kerr.par
index 2a236c5..e26d4fe 100644
--- a/par/Kerr.par
+++ b/par/Kerr.par
@@ -37,6 +37,7 @@ IOUtil::parfile_write = "no"
 
 ########################################
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/bl-bad.par b/par/bl-bad.par
index 449dc94..b9c2d32 100644
--- a/par/bl-bad.par
+++ b/par/bl-bad.par
@@ -69,6 +69,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/bl-different-mask-min-radius-points=5.par b/par/bl-different-mask-min-radius-points=5.par
index e04df91..8c06566 100644
--- a/par/bl-different-mask-min-radius-points=5.par
+++ b/par/bl-different-mask-min-radius-points=5.par
@@ -73,6 +73,7 @@ IOASCII::out2D_vars = "SpaceMask::emask"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/bl-different-mask.par b/par/bl-different-mask.par
index 79720b9..dbd641c 100644
--- a/par/bl-different-mask.par
+++ b/par/bl-different-mask.par
@@ -73,6 +73,7 @@ IOASCII::out2D_vars = "SpaceMask::emask"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/bl-y.par b/par/bl-y.par
index ec4fd6b..6bd203f 100644
--- a/par/bl-y.par
+++ b/par/bl-y.par
@@ -67,6 +67,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/bl.par b/par/bl.par
index 7823b20..d304359 100644
--- a/par/bl.par
+++ b/par/bl.par
@@ -67,6 +67,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner-mask-noshrink.par b/par/misner-mask-noshrink.par
index dd9edd6..440f7cc 100644
--- a/par/misner-mask-noshrink.par
+++ b/par/misner-mask-noshrink.par
@@ -71,6 +71,7 @@ IOASCII::out3D_vars = "SpaceMask::emask"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner-mask-xnoshrink.par b/par/misner-mask-xnoshrink.par
index 381966f..602656d 100644
--- a/par/misner-mask-xnoshrink.par
+++ b/par/misner-mask-xnoshrink.par
@@ -71,6 +71,7 @@ IOASCII::out3D_vars = "SpaceMask::emask"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner-max-allowable-horizon-radius=0.75.par b/par/misner-max-allowable-horizon-radius=0.75.par
index 65e5751..4d3c9f2 100644
--- a/par/misner-max-allowable-horizon-radius=0.75.par
+++ b/par/misner-max-allowable-horizon-radius=0.75.par
@@ -62,6 +62,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner-max-allowable-horizon-radius=0.8.par b/par/misner-max-allowable-horizon-radius=0.8.par
index 341e14d..3925685 100644
--- a/par/misner-max-allowable-horizon-radius=0.8.par
+++ b/par/misner-max-allowable-horizon-radius=0.8.par
@@ -62,6 +62,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner-umfpack.par b/par/misner-umfpack.par
index ffb468d..079135d 100644
--- a/par/misner-umfpack.par
+++ b/par/misner-umfpack.par
@@ -62,6 +62,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/par/misner.par b/par/misner.par
index 87464f8..c733fc4 100644
--- a/par/misner.par
+++ b/par/misner.par
@@ -62,6 +62,7 @@ IOUtil::parfile_write = "no"
 # AHFinderDirect
 #
 
+ActiveThorns = "SphericalSurface"
 ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::print_timing_stats = "true"
diff --git a/param.ccl b/param.ccl
index dc9170d..0397ba6 100644
--- a/param.ccl
+++ b/param.ccl
@@ -26,6 +26,15 @@ USES KEYWORD metric_type
 shares: IO
 USES STRING out_dir
 
+# we need to look at SphericalSurface::nsurfaces to know how many
+# surfaces there are, so that we don't access a non-existing one
+shares: SphericalSurface
+USES INT nsurfaces
+USES INT maxntheta
+USES INT maxnphi
+USES INT    ntheta
+USES INT    nphi
+
 # all remaining parameters are private to this thorn
 private:
 
@@ -94,6 +103,21 @@ int which_horizon_to_announce_centroid \
 } 0
 
 #
+# This parameter controls whether/how we should store our horizon
+# information into the SphericalSurface variables.
+#
+# If this parameter is set to any value >= 0, that value must be
+# in the range  0 <= value < SphericalSurface::nsurfaces .  We check
+# this at run-time...
+#
+int which_surface_to_store_info[101] \
+  "into which SphericalSurface surface should we store our horizon information?"
+{
+-1  :: "don't store this horizon's information"
+0:* :: "store this horizon's information into the specified SphericalSurface surface"
+} -1
+
+#
 # This parameter controls how verbose this thorn is in printing
 # informational (non-error) messages describing what it's doing.
 #
diff --git a/schedule.ccl b/schedule.ccl
index 19eca06..1464fcb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -24,6 +24,12 @@ if (find_every != 0)
 			 lang: C
 			 } "maybe announce horizon position(s) to other thorns"
 		}
+	schedule AHFinderDirect_store_SS_info at CCTK_ANALYSIS \
+	   after AHFinderDirect_find_horizons
+		 {
+		 lang: C
+		 options: global
+		 } "maybe store horizon info in SphericalSurface variables"
 
 	#
 	# *** KLUDGE ***
@@ -31,7 +37,7 @@ if (find_every != 0)
 	# We would really like to give this routine the
 	#	options:global
 	# attribute, so it only runs on those time levels where
-	# we've found (or at least tried to fnid) horizons.  But
+	# we've found (or at least tried to find) horizons.  But
 	# this doesn't work, because we need a GH -- and the other
 	# thorns' routines we call need valid grid variables in
 	# it -- in the announcing code, and  options:global
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,
diff --git a/test/Kerr.par b/test/Kerr.par
index 050bdef..cecd6ce 100644
--- a/test/Kerr.par
+++ b/test/Kerr.par
@@ -38,12 +38,9 @@ IO::out_fileinfo = "none"
 
 ########################################
 
-ActiveThorns = "AEILocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+ActiveThorns = "SphericalSurface AEILocalInterp PUGHInterp PUGHReduce"
+ActiveThorns = "AHFinderDirect"
 
-AHFinderDirect::print_timing_stats = "true"
-
-AHFinderDirect::h_base_file_name              = "h"
-AHFinderDirect::BH_diagnostics_base_file_name = "BH_diagnostics"
 AHFinderDirect::Theta_norm_for_convergence = 1.0e-12
 
 AHFinderDirect::N_horizons = 1