Make "Erik" branch the default trunk for AHFinderDirect.

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

18 files changed, 1898 insertions, 1025 deletions

diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc
index 125177c..dde7d5c 100644
--- a/src/driver/BH_diagnostics.cc
+++ b/src/driver/BH_diagnostics.cc
@@ -12,15 +12,22 @@
 // print - print a line or two summarizing the diagnostics
 // setup_output_file - create/open output file, write header describing fields
 // output - write a (long) line of all the diagnostics
+// store - copy the surface and the diagnostics into the SphericalSurface
+//         variables
+// save - copy the surface and the diagnostics into the Cactus variables
+// load - set the surface and the diagnostics from the Cactus variables
 //
 
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
 
+#include <vector>
+
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -28,6 +35,7 @@
 #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"
@@ -50,7 +58,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 //******************************************************************************
@@ -67,17 +74,28 @@ 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),
-	  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),
+	: origin_x(0.0), origin_y(0.0), origin_z(0.0),
+	  centroid_x(0.0), centroid_y(0.0), centroid_z(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), irreducible_mass(0.0), areal_radius(0.0) // no comma
+	  area(1.0),
+          expansion(0.0),
+          inner_expansion(0.0),
+          product_expansion(0.0),
+          mean_curvature(0.0),
+          area_gradient(0.0),
+          expansion_gradient(0.0),
+          inner_expansion_gradient(0.0),
+          product_expansion_gradient(0.0),
+          mean_curvature_gradient(0.0),
+          mean_curvature_minimum(0.0),
+          mean_curvature_maximum(0.0),
+          mean_curvature_integral(0.0)
 { }
 
 //******************************************************************************
@@ -90,6 +108,10 @@ BH_diagnostics::BH_diagnostics()
 void BH_diagnostics::copy_to_buffer(CCTK_REAL buffer[N_buffer])
 	const
 {
+buffer[posn__origin_x] = this->origin_x;
+buffer[posn__origin_y] = this->origin_y;
+buffer[posn__origin_z] = this->origin_z;
+
 buffer[posn__centroid_x] = centroid_x;
 buffer[posn__centroid_y] = centroid_y;
 buffer[posn__centroid_z] = centroid_z;
@@ -112,13 +134,24 @@ buffer[posn__max_y] = max_y;
 buffer[posn__min_z] = min_z;
 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__irreducible_mass] = irreducible_mass;
-buffer[posn__areal_radius]  = areal_radius;
+buffer[posn__circumference_xy]  = circumference_xy;
+buffer[posn__circumference_xz]  = circumference_xz;
+buffer[posn__circumference_yz]  = circumference_yz;
+buffer[posn__area]              = area;
+buffer[posn__expansion]         = expansion;
+buffer[posn__inner_expansion]   = inner_expansion;
+buffer[posn__product_expansion] = product_expansion;
+buffer[posn__mean_curvature]    = mean_curvature;
+
+buffer[posn__area_gradient]              = area_gradient;
+buffer[posn__expansion_gradient]         = expansion_gradient;
+buffer[posn__inner_expansion_gradient]   = inner_expansion_gradient;
+buffer[posn__product_expansion_gradient] = product_expansion_gradient;
+buffer[posn__mean_curvature_gradient]    = mean_curvature_gradient;
+
+buffer[posn__mean_curvature_minimum]    = mean_curvature_minimum;
+buffer[posn__mean_curvature_maximum]    = mean_curvature_maximum;
+buffer[posn__mean_curvature_integral]   = mean_curvature_integral;
 }
 
 //******************************************************************************
@@ -128,6 +161,10 @@ buffer[posn__areal_radius]  = areal_radius;
 //
 void BH_diagnostics::copy_from_buffer(const CCTK_REAL buffer[N_buffer])
 {
+this->origin_x = buffer[posn__origin_x];
+this->origin_y = buffer[posn__origin_y];
+this->origin_z = buffer[posn__origin_z];
+
 centroid_x = buffer[posn__centroid_x];
 centroid_y = buffer[posn__centroid_y];
 centroid_z = buffer[posn__centroid_z];
@@ -150,13 +187,24 @@ max_y = buffer[posn__max_y];
 min_z = buffer[posn__min_z];
 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];
-   irreducible_mass = buffer[posn__irreducible_mass];
-    areal_radius = buffer[posn__areal_radius];
+ circumference_xy = buffer[posn__circumference_xy];
+ circumference_xz = buffer[posn__circumference_xz];
+ circumference_yz = buffer[posn__circumference_yz];
+             area = buffer[posn__area];
+        expansion = buffer[posn__expansion];
+  inner_expansion = buffer[posn__inner_expansion];
+product_expansion = buffer[posn__product_expansion];
+   mean_curvature = buffer[posn__mean_curvature];
+
+             area_gradient = buffer[posn__area_gradient];
+        expansion_gradient = buffer[posn__expansion_gradient];
+  inner_expansion_gradient = buffer[posn__inner_expansion_gradient];
+product_expansion_gradient = buffer[posn__product_expansion_gradient];
+   mean_curvature_gradient = buffer[posn__mean_curvature_gradient];
+
+   mean_curvature_minimum  = buffer[posn__mean_curvature_minimum];
+   mean_curvature_maximum  = buffer[posn__mean_curvature_maximum];
+   mean_curvature_integral = buffer[posn__mean_curvature_integral];
 }
 
 //******************************************************************************
@@ -165,13 +213,12 @@ circumference_yz = buffer[posn__circumference_yz];
 
 //
 // Given that an apparent horizon has been found, this function computes
-// the black hole diagnostics.
+// various black hole diagnostics.
 //
 // Inputs (gridfns)
-// h					# ghosted
-// one					# nominal
-// global_[xyz]				# nominal
-// global_{xx,xy,xz,yy,yz,zz}		# nominal
+// h		# ghosted
+// one		# nominal
+// global_[xyz]	# nominal
 //
 // Bugs:
 // The computation is rather inefficient -- we make many passes over the
@@ -179,6 +226,16 @@ circumference_yz = buffer[posn__circumference_yz];
 //
 void BH_diagnostics::compute
 	(const patch_system& ps,
+         const fp the_area,
+         const fp mean_expansion,
+         const fp mean_inner_expansion,
+         const fp mean_product_expansion,
+         const fp mean_mean_curvature,
+         const fp the_area_gradient,
+         const fp mean_expansion_gradient,
+         const fp mean_inner_expansion_gradient,
+         const fp mean_product_expansion_gradient,
+         const fp mean_mean_curvature_gradient,
 	 const struct BH_diagnostics_info& BH_diagnostics_info)
 {
 //
@@ -194,19 +251,20 @@ max_radius = h_norms.max_abs_value();
 // xyz bounding box of horizon
 //
 
-// compute bounding box of nominal grid (for stored part of the horizon only)
+// compute bounding box of nominal grid
+// ... this is only the stored part of the horizon if there are symmetries
 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();
 
@@ -262,23 +320,31 @@ 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);
+                                        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);
+                                        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);
+                                        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);
+                                        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);
+                                        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);
+                                        gfns::gfn__global_zz, true, true, true,
+                                        BH_diagnostics_info.integral_method);
+
+
+//
+// originds
+//
+this->origin_x = ps.origin_x();
+this->origin_y = ps.origin_y();
+this->origin_z = ps.origin_z();
 
 
 //
@@ -290,32 +356,52 @@ centroid_z = integral_z / integral_one;
 
 
 //
-// quadrupoles (taken about centroid position)
+// quadrupoles
 //
-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;
+quadrupole_xx = integral_xx / integral_one;
+quadrupole_xy = integral_xy / integral_one;
+quadrupole_xz = integral_xz / integral_one;
+quadrupole_yy = integral_yy / integral_one;
+quadrupole_yz = integral_yz / integral_one;
+quadrupole_zz = integral_zz / integral_one;
 
 
 //
-// mean radius of horizon
+// area, mean radius, and mass
 //
 mean_radius = integral_h / integral_one;
 
 
 //
-// surface area and quantities derived from it
+// expansion
+//
+area              = the_area;
+expansion         = mean_expansion;
+inner_expansion   = mean_inner_expansion;
+product_expansion = mean_product_expansion;
+mean_curvature    = mean_mean_curvature;
+
+area_gradient              = the_area_gradient;
+expansion_gradient         = mean_expansion_gradient;
+inner_expansion_gradient   = mean_inner_expansion_gradient;
+product_expansion_gradient = mean_product_expansion_gradient;
+mean_curvature_gradient    = mean_mean_curvature_gradient;
+
+//
+// minimum, maximum and the integral of the mean curvature
 //
-area             = integral_one;
-irreducible_mass = sqrt(area / (16.0*PI));
-areal_radius     = sqrt(area / ( 4.0*PI));
+jtutil::norm<fp> mean_curvature_norms;
+ps.gridfn_norms(gfns::gfn__mean_curvature, mean_curvature_norms);
+mean_curvature_minimum = mean_curvature_norms.min_value();
+mean_curvature_maximum = mean_curvature_norms.max_value();
+mean_curvature_integral = surface_integral(ps,
+                                           gfns::gfn__mean_curvature,
+                                           true, true, true,
+                                           BH_diagnostics_info.integral_method);
 
 
 //
-// proper circumferences
+// circumferences
 //
 circumference_xy
   = ps.circumference("xy", gfns::gfn__h,
@@ -373,15 +459,16 @@ return ps.integrate_gridfn
 void BH_diagnostics::print(int N_horizons, int hn)
 	const
 {
+const fp m_irreducible = 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 irreducible_mass=%.10g",
+	   "AH %d/%d: area=%.10g m_irreducible=%.10g",
 	   hn, N_horizons,
-	   double(area), double(irreducible_mass));
+	   double(area), double(m_irreducible));
 }
 
 //******************************************************************************
@@ -419,10 +506,13 @@ snprintf(file_name_buffer, IO_info::file_name_buffer_size,
 	 directory, IO_info.BH_diagnostics_base_file_name,
 	 hn, IO_info.BH_diagnostics_file_name_extension);
 
-const char *const file_open_mode = IO_TruncateOutputFiles(state.cgi.GH)
-				   ? "w" : "a";
+const char *openMode;
+if (IO_TruncateOutputFiles(state.cgi.GH) == 1)
+   openMode = "w";
+else
+   openMode = "a";
 
-FILE *fileptr = fopen(file_name_buffer, file_open_mode);
+FILE *fileptr = fopen(file_name_buffer, openMode);
 if (fileptr == NULL)
    then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -460,17 +550,20 @@ 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 27 = m_irreducible\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");
+fprintf(fileptr, "# column 29 = expansion Theta_(l)\n");
+fprintf(fileptr, "# column 30 = inner expansion Theta_(n)\n");
+fprintf(fileptr, "# column 31 = product of the expansions\n");
+fprintf(fileptr, "# column 32 = mean curvature\n");
+fprintf(fileptr, "# column 33 = gradient of the areal radius\n");
+fprintf(fileptr, "# column 34 = gradient of the expansion Theta_(l)\n");
+fprintf(fileptr, "# column 35 = gradient of the inner expansion Theta_(n)\n");
+fprintf(fileptr, "# column 36 = gradient of the product of the expansions\n");
+fprintf(fileptr, "# column 37 = gradient of the mean curvature\n");
+fprintf(fileptr, "# column 38 = minimum  of the mean curvature\n");
+fprintf(fileptr, "# column 39 = maximum  of the mean curvature\n");
+fprintf(fileptr, "# column 40 = integral of the mean curvature\n");
 fflush(fileptr);
 
 return fileptr;
@@ -493,22 +586,25 @@ void BH_diagnostics::output(FILE*fileptr, const struct IO_info& IO_info)
 assert(fileptr != NULL);
 
 fprintf(fileptr,
-     //  cctk_iteration      min radius      mean radius
-     //  ==  cctk_time       ======  max radius
-     //  ==  ==  centroid_[xyz]      ======  ======
-     //  ==  ==  ==========  ======  ======  ======
-	"%d\t%f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t",
+     //  cctk_iteration        min radius      mean radius
+     //  ==  cctk_time         ======  max radius
+     //  ==  ====  centroid_[xyz]      ======  ======
+     //  ==  ====  ==========  ======  ======  ======
+	"%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t",
 	IO_info.time_iteration, double(IO_info.time),
 	double(centroid_x), double(centroid_y), double(centroid_z),
 	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));
+	double(quadrupole_xx - centroid_x*centroid_x),
+        double(quadrupole_xy - centroid_x*centroid_y),
+	double(quadrupole_xz - centroid_x*centroid_z),
+        double(quadrupole_yy - centroid_y*centroid_y),
+	double(quadrupole_yz - centroid_y*centroid_z),
+        double(quadrupole_zz - centroid_z*centroid_z));
 
 fprintf(fileptr,
      //  {min,max}_x     {min,max}_y     {min,max}_z
@@ -530,35 +626,248 @@ fprintf(fileptr,
 	double(circumference_xz / circumference_xy),
 	double(circumference_yz / circumference_xy));
 
+const fp m_irreducible = sqrt(area / (16*PI));;
+const fp areal_radius = sqrt(area / (4*PI));
+const fp areal_radius_gradient = sqrt(1 / (16*PI*area)) * area_gradient;
 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);
+     //  area    m_irre- areal   expan-  inner   prod.   mean    areal   expan-  inner   prod.   mean    mean    mean    mean
+     //          ducible radius  sion    expan-  of the  curva-  radius  sion    expan-  of the  curva-	 curva-  curva-  curva-
+     //                                  sion    expan-  ture    grad.   grad.   sion    exp.s   ture    ture    ture    ture
+     //                                          sions                           grad.   grad.   grad.   min.    max.    integ.
+     //  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======
+     //                                                                                                
+     //  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======  ======
+	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
+	double(area), double(m_irreducible), double(areal_radius),
+        double(expansion), double(inner_expansion), double(product_expansion), double(mean_curvature),
+        double(areal_radius_gradient),
+        double(expansion_gradient), double(inner_expansion_gradient), double(product_expansion_gradient), double(mean_curvature_gradient),
+	double(mean_curvature_minimum), double(mean_curvature_maximum), double(mean_curvature_integral));
 
 fflush(fileptr);
 }
 
 //******************************************************************************
+
+//
+// This function copies the BH-diagnostics into the SphericalSurface variables.
+//
+// Arguments:
+// CCTK_ARGUMENTS
+//
+void BH_diagnostics::store(CCTK_ARGUMENTS,
+                           const int horizon_number, const int surface_number)
+	const
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  assert (surface_number >= 0 && surface_number < nsurfaces);
+  
+  assert(state.AH_data_array[horizon_number] != NULL);
+  const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+  
+  // tell the world if we don't look for this horizon any more
+  const int my_dont_find_after = dont_find_after_individual[horizon_number];
+  const fp my_find_after_time = find_after_individual_time[horizon_number];
+  const fp my_dont_find_after_time = dont_find_after_individual_time[horizon_number];
+  if ((my_dont_find_after >= 0 and cctk_iteration > my_dont_find_after) or
+      (my_dont_find_after_time > my_find_after_time and cctk_time > my_dont_find_after_time))
+  {
+    assert (! AH_data.search_flag);
+    sf_active[surface_number] = 0;
+  }
+  
+  // only try to copy AH info if we've found AHs at this time level
+  if (! AH_data.search_flag) {
+    sf_valid[surface_number] = 0;
+    return;
+  }
+  
+  // did we actually *find* this horizon?
+  if (! AH_data.found_flag) {
+    sf_valid[surface_number] = -1;
+    return;
+  }
+  
+  sf_active       [surface_number] = 1;
+  sf_valid        [surface_number] = 1;
+  sf_origin_x     [surface_number] = this->origin_x;
+  sf_origin_y     [surface_number] = this->origin_y;
+  sf_origin_z     [surface_number] = this->origin_z;
+  sf_mean_radius  [surface_number] = mean_radius;
+  sf_min_radius   [surface_number] = min_radius;
+  sf_max_radius   [surface_number] = max_radius;
+  sf_centroid_x   [surface_number] = centroid_x;
+  sf_centroid_y   [surface_number] = centroid_y;
+  sf_centroid_z   [surface_number] = centroid_z;
+  sf_quadrupole_xx[surface_number] = quadrupole_xx - centroid_x*centroid_x;
+  sf_quadrupole_xy[surface_number] = quadrupole_xy - centroid_x*centroid_y;
+  sf_quadrupole_xz[surface_number] = quadrupole_xz - centroid_x*centroid_z;
+  sf_quadrupole_yy[surface_number] = quadrupole_yy - centroid_y*centroid_y;
+  sf_quadrupole_yz[surface_number] = quadrupole_yz - centroid_y*centroid_z;
+  sf_quadrupole_zz[surface_number] = quadrupole_zz - centroid_z*centroid_z;
+  sf_min_x        [surface_number] = min_x;
+  sf_max_x        [surface_number] = max_x;
+  sf_min_y        [surface_number] = min_y;
+  sf_max_y        [surface_number] = max_y;
+  sf_min_z        [surface_number] = min_z;
+  sf_max_z        [surface_number] = max_z;
+  sf_area         [surface_number] = area;
+  
+  std::vector<CCTK_REAL> xx, yy, zz, rr;
+  xx.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]);
+  yy.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]);
+  zz.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]);
+  rr.resize (sf_ntheta[surface_number] * sf_nphi[surface_number]);
+  size_t n;
+  n = 0;
+  for (int j=0; j<sf_nphi[surface_number]; ++j) {
+    for (int i=0; i<sf_ntheta[surface_number]; ++i) {
+      const CCTK_REAL theta
+        = sf_origin_theta[surface_number] + i * sf_delta_theta[surface_number];
+      const CCTK_REAL phi
+        = sf_origin_phi[surface_number] + j * sf_delta_phi[surface_number];
+      xx[n] = sf_origin_x[surface_number] + sin(theta) * cos(phi);
+      yy[n] = sf_origin_y[surface_number] + sin(theta) * sin(phi);
+      zz[n] = sf_origin_z[surface_number] + cos(theta);
+      ++n;
+    }
+  }
+  AHFinderDirect_radius_in_direction
+    (horizon_number, sf_ntheta[surface_number] * sf_nphi[surface_number],
+     &xx[0], &yy[0], &zz[0], &rr[0]);
+  n = 0;
+  for (int j=0; j<sf_nphi[surface_number]; ++j) {
+    for (int i=0; i<sf_ntheta[surface_number]; ++i) {
+      sf_radius[i + maxntheta * (j + maxnphi * surface_number)] = rr[n];
+      ++n;
+    }
+  }
+}
+
+//******************************************************************************
+
+//
+// This function copies the BH-diagnostics into the Cactus variables.
+//
+// Arguments:
+// CCTK_ARGUMENTS
+//
+void BH_diagnostics::save(CCTK_ARGUMENTS,
+                          const int horizon_number)
+	const
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  assert(state.AH_data_array[horizon_number] != NULL);
+  const struct verbose_info& verbose_info = state.verbose_info;
+  const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+  patch_system& ps = *AH_data.ps_ptr;
+  
+  if (AH_data.search_flag &&
+      AH_data.found_flag &&
+      cctk_iteration > ah_centroid_iteration[horizon_number-1])
+  {
+    ah_centroid_x_p[horizon_number-1] = ah_centroid_x[horizon_number-1];
+    ah_centroid_y_p[horizon_number-1] = ah_centroid_y[horizon_number-1];
+    ah_centroid_z_p[horizon_number-1] = ah_centroid_z[horizon_number-1];
+    ah_centroid_t_p[horizon_number-1] = ah_centroid_t[horizon_number-1];
+    ah_centroid_valid_p[horizon_number-1] = ah_centroid_valid[horizon_number-1];
+    ah_centroid_iteration_p[horizon_number-1] =
+      ah_centroid_iteration[horizon_number-1];
+  
+    ah_centroid_x[horizon_number-1] = centroid_x;
+    ah_centroid_y[horizon_number-1] = centroid_y;
+    ah_centroid_z[horizon_number-1] = centroid_z;
+    ah_centroid_t[horizon_number-1] = cctk_time;
+    ah_centroid_valid[horizon_number-1] = AH_data.found_flag;
+    ah_centroid_iteration[horizon_number-1] = cctk_iteration;
+  }
+  
+  ah_origin_x[horizon_number-1] = ps.origin_x();
+  ah_origin_y[horizon_number-1] = ps.origin_y();
+  ah_origin_z[horizon_number-1] = ps.origin_z();
+  
+  for (int pn = 0; pn < ps.N_patches(); ++pn) {
+    assert (pn < 6);
+    patch& p = ps.ith_patch(pn);
+    for (int i = 0, irho = p.min_irho(); irho <= p.max_irho(); ++i, ++irho) {
+      assert (i <= max_N_zones_per_right_angle);
+      for (int j = 0, isigma = p.min_isigma(); isigma <= p.max_isigma(); ++j, ++isigma) {
+        assert (j <= max_N_zones_per_right_angle);
+        ah_radius[i + (max_N_zones_per_right_angle+1) * (j + (max_N_zones_per_right_angle+1) * (pn + 6 * (horizon_number-1)))]
+          = p.ghosted_gridfn(gfns::gfn__h, irho,isigma);
+      }
+    }
+  }
+  
+  ah_initial_find_flag       [horizon_number-1] = AH_data.initial_find_flag;
+  ah_really_initial_find_flag[horizon_number-1] = AH_data.really_initial_find_flag;
+  ah_search_flag             [horizon_number-1] = AH_data.search_flag;
+  ah_found_flag              [horizon_number-1] = AH_data.found_flag;
+  if (verbose_info.print_algorithm_highlights) {
+    printf ("AHF BH_diagnostics::save[%d] initial_find_flag=%d\n",        horizon_number, (int) AH_data.initial_find_flag);
+    printf ("AHF BH_diagnostics::save[%d] really_initial_find_flag=%d\n", horizon_number, (int) AH_data.really_initial_find_flag);
+    printf ("AHF BH_diagnostics::save[%d] search_flag=%d\n",              horizon_number, (int) AH_data.search_flag);
+    printf ("AHF BH_diagnostics::save[%d] found_flag=%d\n",               horizon_number, (int) AH_data.found_flag);
+  }
+}
+
+//******************************************************************************
+
+//
+// This function copies the BH-diagnostics from the Cactus variables.
+//
+// Arguments:
+// CCTK_ARGUMENTS
+//
+void BH_diagnostics::load(CCTK_ARGUMENTS,
+                          const int horizon_number)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  assert(state.AH_data_array[horizon_number] != NULL);
+  const struct verbose_info& verbose_info = state.verbose_info;
+  struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+  patch_system& ps = *AH_data.ps_ptr;
+  
+  ps.origin_x(ah_origin_x[horizon_number-1]);
+  ps.origin_y(ah_origin_y[horizon_number-1]);
+  ps.origin_z(ah_origin_z[horizon_number-1]);
+  
+  for (int pn = 0; pn < ps.N_patches(); ++pn) {
+    assert (pn < 6);
+    patch& p = ps.ith_patch(pn);
+    for (int i = 0, irho = p.min_irho(); irho <= p.max_irho(); ++i, ++irho) {
+      assert (i <= max_N_zones_per_right_angle);
+      for (int j = 0, isigma = p.min_isigma(); isigma <= p.max_isigma(); ++j, ++isigma) {
+        assert (j <= max_N_zones_per_right_angle);
+        p.ghosted_gridfn(gfns::gfn__h, irho,isigma) =
+          ah_radius[i + (max_N_zones_per_right_angle+1) * (j + (max_N_zones_per_right_angle+1) * (pn + 6 * (horizon_number-1)))];
+      }
+    }
+  }
+  
+  // recover the full ghosted-grid horizon shape
+  // (we only save and load the nominal-grid shape)
+  ps.synchronize();
+  
+  AH_data.initial_find_flag        = ah_initial_find_flag       [horizon_number-1];
+  AH_data.really_initial_find_flag = ah_really_initial_find_flag[horizon_number-1];
+  AH_data.search_flag              = ah_search_flag             [horizon_number-1];
+  AH_data.found_flag               = ah_found_flag              [horizon_number-1];
+  if (verbose_info.print_algorithm_highlights) {
+    printf ("AHF BH_diagnostics::load[%d] initial_find_flag=%d\n",        horizon_number, (int) AH_data.initial_find_flag);
+    printf ("AHF BH_diagnostics::load[%d] really_initial_find_flag=%d\n", horizon_number, (int) AH_data.really_initial_find_flag);
+    printf ("AHF BH_diagnostics::load[%d] search_flag=%d\n",              horizon_number, (int) AH_data.search_flag);
+    printf ("AHF BH_diagnostics::load[%d] found_flag=%d\n",               horizon_number, (int) AH_data.found_flag);
+  }
+}
+
+//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
diff --git a/src/driver/BH_diagnostics.hh b/src/driver/BH_diagnostics.hh
index 0d368e8..6698b05 100644
--- a/src/driver/BH_diagnostics.hh
+++ b/src/driver/BH_diagnostics.hh
@@ -4,7 +4,7 @@
 //
 // prerequisites:
 //	<stdio.h>
-//	"cctk.h"
+//	"cctk_Arguments.h"
 //
 
 // everything in this file is inside this namespace
@@ -31,57 +31,63 @@ 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 origin_x, origin_y, origin_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 quadrupole_xx, quadrupole_xy, quadrupole_xz;
+        fp quadrupole_yy, quadrupole_yz, quadrupole_zz;
 	fp min_radius, max_radius, mean_radius;
 
 	// xyz bounding box
-	fp min_x, max_x,
-	   min_y, max_y,
-	   min_z, max_z;
-
-	// 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;
+	fp min_x, max_x, min_y, max_y, min_z, max_z;
+
+        fp circumference_xy, circumference_xz, circumference_yz;
+        fp area;
+        fp expansion;
+        fp inner_expansion;
+        fp product_expansion;
+        fp mean_curvature;
+
+        fp area_gradient;
+        fp expansion_gradient;
+        fp inner_expansion_gradient;
+        fp product_expansion_gradient;
+        fp mean_curvature_gradient;
+        fp mean_curvature_minimum;
+        fp mean_curvature_maximum;
+        fp mean_curvature_integral;
 
 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__origin_x = 0, posn__origin_y, posn__origin_z,
+		posn__centroid_x, 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__irreducible_mass, posn__areal_radius,
+		posn__circumference_xy, posn__circumference_xz,
+					posn__circumference_yz,
+		posn__area,
+                posn__expansion,
+                posn__inner_expansion,
+                posn__product_expansion,
+                posn__mean_curvature,
+
+		posn__area_gradient,
+                posn__expansion_gradient,
+                posn__inner_expansion_gradient,
+                posn__product_expansion_gradient,
+                posn__mean_curvature_gradient,
+                posn__mean_curvature_minimum,
+                posn__mean_curvature_maximum,
+                posn__mean_curvature_integral,
 
 		N_buffer // no comma	// size of buffer
 		};
@@ -93,6 +99,16 @@ public:
 public:
 	// compute diagnostics (assuming that apparent horizon has been found)
 	void compute(const patch_system& ps,
+                     const fp area,
+                     const fp mean_expansion,
+                     const fp mean_inner_expansion,
+                     const fp mean_product_expansion,
+                     const fp mean_mean_curvature,
+                     const fp area_gradient,
+                     const fp mean_expansion_gradient,
+                     const fp mean_inner_expansion_gradient,
+                     const fp mean_product_expansion_gradient,
+                     const fp mean_mean_curvature_gradient,
 		     const struct BH_diagnostics_info& BH_diagnostics_info);
 
 	// print (CCTK_VInfo()) a line or two summarizing diagnostics
@@ -112,6 +128,17 @@ public:
 	void output(FILE* fileptr, const struct IO_info& IO_info)
 		const;
 
+	// store the diagnostics in the Cactus variables
+        void store(CCTK_ARGUMENTS, int horizon_number, int surface_number)
+                const;
+
+	// save the horizon shape in grid arrays
+        void save(CCTK_ARGUMENTS, int horizon_number)
+                const;
+
+	// load the horizon shape from grid arrays
+        void load(CCTK_ARGUMENTS, int horizon_number);
+
 	// constructor initializes all diagnostics to 0.0
 	BH_diagnostics();
 
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 2a54d17..62bab31 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -10,15 +10,19 @@
 /// Newton_step - take the Newton step, scaling it down if it's too large
 //
 
+#include <iostream>
+
 #include <stdio.h>
 #include <assert.h>
 #include <limits.h>
 #include <float.h>
 #include <math.h>
+#include <string.h>
 
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -26,6 +30,7 @@
 #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"
@@ -48,7 +53,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -67,6 +71,7 @@ bool broadcast_status(const cGH* GH,
 		      struct iteration_status_buffers& isb);
 void broadcast_horizon_data(const cGH* GH,
 			    bool broadcast_flag, bool broadcast_horizon_shape,
+                            struct AH_data& AH_data,
 			    struct BH_diagnostics& BH_diagnostics,
 			    patch_system& ps,
 			    struct horizon_buffers& horizon_buffers);
@@ -106,22 +111,22 @@ void Newton(const cGH* GH,
 	    const struct verbose_info& verbose_info,
 	    struct iteration_status_buffers& isb)
 {
+cGH const * const cctkGH = GH;
+DECLARE_CCTK_ARGUMENTS;
+DECLARE_CCTK_PARAMETERS;
+
 const bool my_active_flag = hs.has_genuine_horizons();
-const int N_horizons = hs.N_horizons();
 
 // print out which horizons we're finding on this processor
-if (verbose_info.print_physics_details)
-   then {
-	if (hs.has_genuine_horizons())
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "proc %d: searching for horizon%s %s/%d",
-			   my_proc,
-			   (hs.my_N_horizons() > 1 ? "s" : ""),
-			   hs.sequence_string(","), N_horizons);
-	   else CCTK_VInfo(CCTK_THORNSTRING,
-			   "proc %d: dummy horizon(s) only",
-			   my_proc);
-	}
+if (hs.has_genuine_horizons())
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "proc %d: searching for horizon%s %s/%d",
+		   my_proc,
+		   (hs.my_N_horizons() > 1 ? "s" : ""),
+		   hs.sequence_string(","), N_horizons);
+   else CCTK_VInfo(CCTK_THORNSTRING,
+		   "proc %d: dummy horizon(s) only",
+		   my_proc);
 
     //
     // each pass through this loop finds a single horizon,
@@ -138,7 +143,10 @@ if (verbose_info.print_physics_details)
 		       "Newton_solve(): processor %d working on horizon %d",
 		       my_proc, hn);
 
-    const bool horizon_is_genuine               = hs.is_genuine();
+    // only try to find horizons every  find_every  time steps
+    const bool horizon_is_genuine =
+      hs.is_genuine() && AH_data_array[hn]->search_flag;
+    // this is only a pessimistic approximation
     const bool there_is_another_genuine_horizon = hs.is_next_genuine();
     if (verbose_info.print_algorithm_details)
        then {
@@ -152,12 +160,262 @@ if (verbose_info.print_physics_details)
 
     struct AH_data* AH_data_ptr
 	                  = horizon_is_genuine ? AH_data_array[hn]   : NULL;
+
     patch_system* const  ps_ptr = horizon_is_genuine
 				  ? AH_data_ptr->ps_ptr : NULL;
     Jacobian*     const Jac_ptr = horizon_is_genuine
 				  ? AH_data_ptr->Jac_ptr: NULL;
-    const fp add_to_expansion = horizon_is_genuine
-				? -AH_data_ptr->surface_expansion : 0.0;
+
+    if (horizon_is_genuine) {
+      // deal with dependent horizons
+      if (AH_data_ptr->depends_on > 0) {
+        assert (AH_data_ptr->depends_on != hn);
+        assert (AH_data_ptr->depends_on < hn);
+        // check that the other horizon is on the same processor!
+        AH_data *AH_other_ptr = AH_data_array[AH_data_ptr->depends_on];
+        assert (AH_other_ptr);
+        AH_data_ptr->compute_info.desired_value
+          = AH_other_ptr->compute_info.desired_value
+          * AH_data_ptr->desired_value_factor
+          + AH_data_ptr->desired_value_offset;
+        const int gnp = ps_ptr->ghosted_N_grid_points();
+        assert (AH_other_ptr->ps_ptr->ghosted_N_grid_points() == gnp);
+        memcpy (ps_ptr->ghosted_gridfn_data(gfns::gfn__h),
+                AH_other_ptr->ps_ptr->ghosted_gridfn_data(gfns::gfn__h),
+                gnp * sizeof(fp));
+        ps_ptr->origin_x (AH_other_ptr->ps_ptr->origin_x());
+        ps_ptr->origin_y (AH_other_ptr->ps_ptr->origin_y());
+        ps_ptr->origin_z (AH_other_ptr->ps_ptr->origin_z());
+      }
+    }
+
+    if (horizon_is_genuine) {
+      if (AH_data_ptr->move_origins
+          && AH_data_ptr->depends_on == 0
+          && AH_data_ptr->found_flag)
+      {
+        // move the origin to the centre
+        patch_system& ps = *ps_ptr;
+        fp cx=0, cy=0, cz=0;    // change of origin
+        fp sx=0, sy=0, sz=0;    // change of shape
+        if (! predict_origin_movement) {
+          size_t N=0;
+          for (int pn=0; pn<ps.N_patches(); ++pn) {
+            patch& p = ps.ith_patch(pn);
+            for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) {
+              for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) {
+                const fp rho = p.rho_of_irho (irho);
+                const fp sigma = p.sigma_of_isigma (isigma);
+                const fp r = p.ghosted_gridfn (gfns::gfn__h, irho, isigma);
+                fp x, y, z;
+                p.xyz_of_r_rho_sigma (r, rho, sigma, x, y, z);
+                cx += x;
+                cy += y;
+                cz += z;
+                ++N;
+              }
+            }
+          }
+          assert (N > 0);
+          cx /= N;
+          cy /= N;
+          cz /= N;
+          sx=cx; sy=cy; sz=cz;
+        } else {                // if predict_origin_movement
+          if (ah_centroid_valid[hn-1]) {
+            if (verbose_info.print_algorithm_highlights) {
+              std::cout << "AHF centroid x " << (ah_centroid_x[hn-1]) << " y " << (ah_centroid_y[hn-1]) << " z " << (ah_centroid_z[hn-1]) << " t " << (ah_centroid_t[hn-1]) << std::endl;
+            }
+            if (ah_centroid_valid_p[hn-1]) {
+              // have two previous origins: linear extrapolation
+              if (verbose_info.print_algorithm_highlights) {
+                std::cout << "AHF centroid_p x " << (ah_centroid_x_p[hn-1]) << " y " << (ah_centroid_y_p[hn-1]) << " z " << (ah_centroid_z_p[hn-1]) << " t " << (ah_centroid_t_p[hn-1]) << std::endl;
+              }
+              fp const dt   = ah_centroid_t[hn-1] - ah_centroid_t_p[hn-1];
+              fp const dt_n = cctk_time - ah_centroid_t  [hn-1];
+              fp const dt_p = cctk_time - ah_centroid_t_p[hn-1];
+              fp const timescale =
+                fabs (dt) + fabs (dt_p) + fabs (dt_n) +
+                0.001 * fabs (cctk_delta_time);
+              if (10 * fabs (dt  ) < timescale ||
+                  10 * fabs (dt_p) < timescale ||
+                  10 * fabs (dt_n) < timescale)
+              {
+                // the times are too similar
+                if (verbose_info.print_algorithm_highlights) {
+                  std::cout << "AHF toosim" << std::endl;
+                }
+                cx = ah_centroid_x[hn-1];
+                cy = ah_centroid_y[hn-1];
+                cz = ah_centroid_z[hn-1];
+              } else {
+                fp const fact_n = + dt_p / dt;
+                fp const fact_p = - dt_n / dt;
+                if (fabs (fact_n) > 5 || fabs (fact_p) > 5) {
+                  // don't trust a large extrapolation
+                  if (verbose_info.print_algorithm_highlights) {
+                    std::cout << "AHF notrust" << std::endl;
+                  }
+                  cx = ah_centroid_x[hn-1];
+                  cy = ah_centroid_y[hn-1];
+                  cz = ah_centroid_z[hn-1];
+                } else {
+                  if (verbose_info.print_algorithm_highlights) {
+                    std::cout << "AHF extrap fn " << fact_n << " fp " << fact_p << std::endl;
+                  }
+                  cx = fact_n * ah_centroid_x[hn-1] + fact_p * ah_centroid_x_p[hn-1];
+                  cy = fact_n * ah_centroid_y[hn-1] + fact_p * ah_centroid_y_p[hn-1];
+                  cz = fact_n * ah_centroid_z[hn-1] + fact_p * ah_centroid_z_p[hn-1];
+                  if (verbose_info.print_algorithm_highlights) {
+                    std::cout << "AHF xp " << (ah_centroid_x_p[hn-1]) << " x " << (ah_centroid_x[hn-1]) << " cx " << cx << std::endl;
+                  }
+                }
+              }
+            } else {
+              // have one previous origin: constant extrapolation
+              if (verbose_info.print_algorithm_highlights) {
+                std::cout << "AHF const" << std::endl;
+              }
+              cx = ah_centroid_x[hn-1];
+              cy = ah_centroid_y[hn-1];
+              cz = ah_centroid_z[hn-1];
+            }
+            if (verbose_info.print_algorithm_highlights) {
+              std::cout << "AHF predicted position cx " << cx << " cy " << cy << " cz " << cz << std::endl;
+            }
+            cx -= ps.origin_x();
+            cy -= ps.origin_y();
+            cz -= ps.origin_z();
+          }
+          sx = ah_centroid_x[hn-1] - ps.origin_x();
+          sy = ah_centroid_y[hn-1] - ps.origin_y();
+          sz = ah_centroid_z[hn-1] - ps.origin_z();
+        } // if predict_origin_movement
+        switch (ps.type()) {
+        case patch_system::patch_system__full_sphere:
+          break;                // do nothing
+        case patch_system::patch_system__plus_z_hemisphere:
+          cz = 0; sz = 0; break;
+        case patch_system::patch_system__plus_xy_quadrant_mirrored:
+          cx = cy = 0; sx = sy = 0; break;
+        case patch_system::patch_system__plus_xy_quadrant_rotating:
+          cx = cy = 0; sx = sy = 0; break;
+        case patch_system::patch_system__plus_xz_quadrant_mirrored:
+          cx = cz = 0; sx = sz = 0; break;
+        case patch_system::patch_system__plus_xz_quadrant_rotating:
+          cx = cz = 0; sx = sz = 0; break;
+        case patch_system::patch_system__plus_xyz_octant_mirrored:
+          cx = cy = cz = 0; sx = sy = sz = 0; break;
+        case patch_system::patch_system__plus_xyz_octant_rotating:
+          cx = cy = cz = 0; sx = sy = sz = 0; break;
+        default: assert(0);
+        }
+        ps.origin_x (ps.origin_x() + cx);
+        ps.origin_y (ps.origin_y() + cy);
+        ps.origin_z (ps.origin_z() + cz);
+        if (reshape_while_moving) {
+          for (int pn=0; pn<ps.N_patches(); ++pn) {
+            patch& p = ps.ith_patch(pn);
+            for (int irho = p.min_irho(); irho <= p.max_irho(); ++irho) {
+              for (int isigma = p.min_isigma(); isigma <= p.max_isigma(); ++isigma) {
+                const fp rho = p.rho_of_irho (irho);
+                const fp sigma = p.sigma_of_isigma (isigma);
+                fp & r = p.ghosted_gridfn (gfns::gfn__h, irho, isigma);
+                fp x, y, z;
+                p.xyz_of_r_rho_sigma (r, rho, sigma, x, y, z);
+                fp const proj = (sx*x + sy*y + sz*z) / r;
+                r -= proj;
+              }
+            }
+          }
+          ps.synchronize();
+        }
+      }
+      
+      // modify the initial guess
+      jtutil::norm<fp> norms;
+      ps_ptr->ghosted_gridfn_norms (gfns::gfn__h, norms);
+      // smoothing:
+      ps_ptr->scale_ghosted_gridfn
+        (1.0 - AH_data_ptr->smoothing_factor, gfns::gfn__h);
+      ps_ptr->add_to_ghosted_gridfn
+        (AH_data_ptr->smoothing_factor * norms.mean(), gfns::gfn__h);
+      // enlarging:
+      ps_ptr->scale_ghosted_gridfn
+        (AH_data_ptr->shiftout_factor, gfns::gfn__h);
+    }
+
+    bool I_am_pretracking = horizon_is_genuine && AH_data_ptr->use_pretracking;
+    bool I_was_pretracking = false;
+    bool pretracking_have_upper_bound = false;
+    bool pretracking_have_lower_bound = false;
+    bool pretracking_was_successful = false;
+    fp const old_pretracking_value = I_am_pretracking ? AH_data_ptr->pretracking_value : 0.0;
+    fp pretracking_upper_bound;
+    fp pretracking_lower_bound;
+    bool pretracking_have_horizon_info;
+    fp pretracking_mean_expansion;
+    for (int pretracking_iter = 0;
+         I_am_pretracking
+           ? pretracking_iter < AH_data_ptr->pretracking_max_iterations
+           : true;
+         ++pretracking_iter)
+      {
+        bool found_this_horizon;
+        if (I_am_pretracking || I_was_pretracking) {
+          CCTK_VInfo(CCTK_THORNSTRING,
+                     "Pretracking: iteration %d", pretracking_iter);
+          if (pretracking_have_lower_bound) {
+            if (pretracking_have_upper_bound) {
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: looking for value %g in [%g,%g]",
+                         double(AH_data_ptr->pretracking_value),
+                           double(pretracking_lower_bound),
+                         double(pretracking_upper_bound));
+            } else {
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: looking for value %g in [%g,*]",
+                         double(AH_data_ptr->pretracking_value),
+                         double(pretracking_lower_bound));
+            }
+          } else {
+            if (pretracking_have_upper_bound) {
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: looking for value %g in [*,%g]",
+                         double(AH_data_ptr->pretracking_value),
+                           double(pretracking_upper_bound));
+            } else {
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: looking for value %g in [*,*]",
+                         double(AH_data_ptr->pretracking_value));
+            }
+          }
+          AH_data_ptr->compute_info.desired_value = AH_data_ptr->pretracking_value;
+          ps_ptr->ghosted_gridfn_copy (gfns::gfn__h, gfns::gfn__save_h);
+          pretracking_have_horizon_info = false;
+        }
+
+        if (! I_am_pretracking) {
+          if (horizon_is_genuine) {
+            if (! AH_data_ptr->initial_guess_info.reset_horizon_after_not_finding) {
+              // save the surface for possible backtracking later
+              ps_ptr->ghosted_gridfn_copy (gfns::gfn__h, gfns::gfn__save_h);
+            }
+          }
+        }
+
+	struct what_to_compute dummy_compute_info;
+	struct what_to_compute & compute_info =
+		horizon_is_genuine
+		? AH_data_ptr->compute_info
+		: dummy_compute_info;
+
+    if (horizon_is_genuine) {
+      if (ps_ptr->N_additional_points()) {
+        const int gnp = ps_ptr->ghosted_N_grid_points();
+        ps_ptr->ghosted_gridfn_data(gfns::gfn__h)[gnp] = 0;
+      }
+    }
 
     const int max_iterations
        = horizon_is_genuine
@@ -166,10 +424,16 @@ if (verbose_info.print_physics_details)
 	    : solver_info.max_Newton_iterations__subsequent)
 	 : INT_MAX;
 
+    int num_Theta_growth_iterations = 0;
+    fp previous_Theta_norm = 1.0e30;
+    int num_Theta_nonshrink_iterations = 0;
+    fp best_Theta_norm = 1.0e30;
+
 	//
 	// each pass through this loop does a single Newton iteration
 	// on the current horizon (which might be either genuine or dummy)
 	//
+        bool do_return = false;
 	for (int iteration = 1 ; ; ++iteration)
 	{
 	if (verbose_info.print_algorithm_debug)
@@ -185,20 +449,136 @@ if (verbose_info.print_physics_details)
 	if (horizon_is_genuine
 	    && solver_info.debugging_output_at_each_Newton_iteration)
 	   then output_gridfn(*ps_ptr, gfns::gfn__h,
+                              "h", GH,
 			      IO_info, IO_info.h_base_file_name,
+                              IO_info.h_min_digits,
 			      hn, verbose_info.print_algorithm_highlights,
 			      iteration);
 
+        // calculate the norms also for a surface a bit further out,
+        // so that we know how they vary in space.
+        // perform this calculation first, so that the real values
+        // do not have to be saved.
+        const fp epsilon = Jacobian_info.perturbation_amplitude;
+	jtutil::norm<fp> shifted_Theta_norms;
+        jtutil::norm<fp> shifted_expansion_Theta_norms;
+        jtutil::norm<fp> shifted_inner_expansion_Theta_norms;
+        jtutil::norm<fp> shifted_product_expansion_Theta_norms;
+        jtutil::norm<fp> shifted_mean_curvature_Theta_norms;
+        fp shifted_area;
+	jtutil::norm<fp> shifted2_Theta_norms;
+        jtutil::norm<fp> shifted2_expansion_Theta_norms;
+        jtutil::norm<fp> shifted2_inner_expansion_Theta_norms;
+        jtutil::norm<fp> shifted2_product_expansion_Theta_norms;
+        jtutil::norm<fp> shifted2_mean_curvature_Theta_norms;
+        fp shifted2_area;
+        if (solver_info.want_expansion_gradients) {
+          if (horizon_is_genuine) {
+            ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+            // ps_ptr->scale_ghosted_gridfn(1.0+epsilon, gfns::gfn__h);
+          }
+          const enum expansion_status raw_shifted_expansion_status
+            = expansion(ps_ptr, compute_info,
+                        cgi, gi,
+                        error_info, (iteration == 1),
+                        false,	// no, we don't want Jacobian coeffs
+                        false,
+                        &shifted_Theta_norms,
+                        &shifted_expansion_Theta_norms,
+                        &shifted_inner_expansion_Theta_norms,
+                        &shifted_product_expansion_Theta_norms,
+                        &shifted_mean_curvature_Theta_norms);
+          if (horizon_is_genuine) {
+            shifted_area = ps_ptr->integrate_gridfn
+              (gfns::gfn__one, true, true, true,
+               gfns::gfn__h,
+               gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+                                   gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+                                                       gfns::gfn__g_dd_33,
+               BH_diagnostics_info.integral_method);
+            ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+            // ps_ptr->scale_ghosted_gridfn(1.0/(1.0+epsilon), gfns::gfn__h);
+            
+            ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+            // ps_ptr->scale_ghosted_gridfn(1.0/(1.0+epsilon), gfns::gfn__h);
+          }
+          const enum expansion_status raw_shifted2_expansion_status
+            = expansion(ps_ptr, compute_info,
+                        cgi, gi,
+                        error_info, (iteration == 1),
+                        false,	// no, we don't want Jacobian coeffs
+                        false,
+                        &shifted2_Theta_norms,
+                        &shifted2_expansion_Theta_norms,
+                        &shifted2_inner_expansion_Theta_norms,
+                        &shifted2_product_expansion_Theta_norms,
+                        &shifted2_mean_curvature_Theta_norms);
+          if (horizon_is_genuine) {
+            shifted2_area = ps_ptr->integrate_gridfn
+              (gfns::gfn__one, true, true, true,
+               gfns::gfn__h,
+               gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+                                   gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+                                                       gfns::gfn__g_dd_33,
+               BH_diagnostics_info.integral_method);
+            ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+            // ps_ptr->scale_ghosted_gridfn(1.0+epsilon, gfns::gfn__h);
+          }
+        } // if want_expansion_gradients
+
+        // now calculate the real values.
 	jtutil::norm<fp> Theta_norms;
+        jtutil::norm<fp> expansion_Theta_norms;
+        jtutil::norm<fp> inner_expansion_Theta_norms;
+        jtutil::norm<fp> product_expansion_Theta_norms;
+        jtutil::norm<fp> mean_curvature_Theta_norms;
 	const enum expansion_status raw_expansion_status
-		= expansion(ps_ptr, add_to_expansion,
+               = expansion(ps_ptr, compute_info,
 			   cgi, gi,
 			   error_info, (iteration == 1),
 			   true,	// yes, we want Jacobian coeffs
 			   verbose_info.print_algorithm_details,
-			   &Theta_norms);
+                           &Theta_norms,
+                           &expansion_Theta_norms,
+                           &inner_expansion_Theta_norms,
+                           &product_expansion_Theta_norms,
+                           &mean_curvature_Theta_norms);
+        fp area;
+        if (horizon_is_genuine)
+	        {
+	        area = ps_ptr->integrate_gridfn
+			(gfns::gfn__one, true, true, true,
+			 gfns::gfn__h,
+		         gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+					     gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+							         gfns::gfn__g_dd_33,
+	                 BH_diagnostics_info.integral_method);
+	        }
 	const bool Theta_is_ok = (raw_expansion_status == expansion_success);
-	const bool norms_are_ok = horizon_is_genuine && Theta_is_ok;
+        const bool norms_are_ok = horizon_is_genuine && Theta_is_ok;
+
+        if (norms_are_ok) {
+          const fp this_norm = Theta_norms.infinity_norm();
+          if (this_norm > previous_Theta_norm) {
+            ++ num_Theta_growth_iterations;
+          } else {
+            num_Theta_growth_iterations = 0;
+          }
+          previous_Theta_norm = this_norm;
+          
+          if (this_norm >= best_Theta_norm) {
+            ++ num_Theta_nonshrink_iterations;
+          } else {
+            num_Theta_nonshrink_iterations = 0;
+            best_Theta_norm = this_norm;
+          }
+        }
+
+        if (I_am_pretracking && norms_are_ok) {
+          pretracking_mean_expansion = expansion_Theta_norms.mean();
+          pretracking_have_horizon_info = true;
+        }
+
 	if (verbose_info.print_algorithm_debug)
 	   then {
 		CCTK_VInfo(CCTK_THORNSTRING,
@@ -212,27 +592,64 @@ if (verbose_info.print_physics_details)
 		}
 	if (horizon_is_genuine && Theta_is_ok
 	    && solver_info.debugging_output_at_each_Newton_iteration)
-	   then output_gridfn(*ps_ptr, gfns::gfn__Theta,
+	   then {
+		output_gridfn(*ps_ptr, gfns::gfn__Theta,
+                              "Theta", GH,
 			      IO_info, IO_info.Theta_base_file_name,
+                              IO_info.h_min_digits,
+			      hn, verbose_info.print_algorithm_highlights,
+			      iteration);
+		output_gridfn(*ps_ptr, gfns::gfn__mean_curvature,
+                              "mean_curvature", GH,
+			      IO_info, IO_info.mean_curvature_base_file_name,
+                              IO_info.h_min_digits,
 			      hn, verbose_info.print_algorithm_highlights,
 			      iteration);
+		}
 
 
 	//
 	// have we found this horizon?
 	// if so, compute and output BH diagnostics
 	//
-	const bool found_this_horizon
-	   = norms_are_ok && (Theta_norms.infinity_norm()
-			      <= solver_info.Theta_norm_for_convergence);
+        found_this_horizon
+           = (norms_are_ok
+              && (I_was_pretracking
+                  ? pretracking_was_successful
+                  : Theta_norms.infinity_norm() <= solver_info.Theta_norm_for_convergence));
 	if (horizon_is_genuine)
 	   then AH_data_ptr->found_flag = found_this_horizon;
 
-	if (found_this_horizon)
+	if (found_this_horizon && ! I_am_pretracking)
 	   then {
 		struct BH_diagnostics& BH_diagnostics
 			= AH_data_ptr->BH_diagnostics;
-		BH_diagnostics.compute(*ps_ptr, BH_diagnostics_info);
+                const fp mean_expansion       = expansion_Theta_norms.mean();
+                const fp mean_inner_expansion = inner_expansion_Theta_norms.mean();
+                const fp mean_product_expansion = product_expansion_Theta_norms.mean();
+                const fp mean_mean_curvature  = mean_curvature_Theta_norms.mean();
+                // const fp area_gradient                 = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_area                               - area                              ) / epsilon;
+                // const fp mean_expansion_gradient       = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_expansion_Theta_norms.mean()       - expansion_Theta_norms.mean()      ) / epsilon;
+                // const fp mean_inner_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_inner_expansion_Theta_norms.mean() - inner_expansion_Theta_norms.mean()) / epsilon;
+                // const fp mean_product_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_product_expansion_Theta_norms.mean() - product_expansion_Theta_norms.mean()) / epsilon;
+                // const fp mean_mean_curvature_gradient  = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_mean_curvature_Theta_norms.mean()  - mean_curvature_Theta_norms.mean() ) / epsilon;
+                const fp area_gradient                 = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_area                               - shifted2_area                              ) / (2*epsilon);
+                const fp mean_expansion_gradient       = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_expansion_Theta_norms.mean()       - shifted2_expansion_Theta_norms.mean()      ) / (2*epsilon);
+                const fp mean_inner_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_inner_expansion_Theta_norms.mean() - shifted2_inner_expansion_Theta_norms.mean()) / (2*epsilon);
+                const fp mean_product_expansion_gradient = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_product_expansion_Theta_norms.mean() - shifted2_product_expansion_Theta_norms.mean()) / (2*epsilon);
+                const fp mean_mean_curvature_gradient  = ! solver_info.want_expansion_gradients ? 0.0 : (shifted_mean_curvature_Theta_norms.mean()  - shifted2_mean_curvature_Theta_norms.mean() ) / (2*epsilon);
+		BH_diagnostics.compute(*ps_ptr,
+                                       area,
+                                       mean_expansion,
+                                       mean_inner_expansion,
+                                       mean_product_expansion,
+                                       mean_mean_curvature,
+                                       area_gradient,
+                                       mean_expansion_gradient,
+                                       mean_inner_expansion_gradient,
+                                       mean_product_expansion_gradient,
+                                       mean_mean_curvature_gradient,
+                                       BH_diagnostics_info);
 		if (IO_info.output_BH_diagnostics)
 		   then {
 			if (AH_data_ptr->BH_diagnostics_fileptr == NULL)
@@ -251,8 +668,13 @@ if (verbose_info.print_physics_details)
 	// (if so, we'll give up on this horizon)
 	//
 	const bool expansion_is_too_large
-		= norms_are_ok && (Theta_norms.infinity_norm()
-				   > solver_info.max_allowable_Theta);
+		= norms_are_ok
+                  && (   Theta_norms.infinity_norm() > solver_info.max_allowable_Theta
+                      || (   solver_info.max_allowable_Theta_growth_iterations > 0
+                          && num_Theta_growth_iterations > solver_info.max_allowable_Theta_growth_iterations)
+                      || (   solver_info.max_allowable_Theta_nonshrink_iterations > 0
+                          && num_Theta_nonshrink_iterations > solver_info.max_allowable_Theta_nonshrink_iterations)
+                         );
 
 
 	//
@@ -260,10 +682,11 @@ if (verbose_info.print_physics_details)
 	// then pretend expansion() returned a "surface too large" error status
 	//
 	jtutil::norm<fp> h_norms;
-	if (horizon_is_genuine)
-	   then ps_ptr->ghosted_gridfn_norms(gfns::gfn__h, h_norms);
-	const fp mean_horizon_radius = horizon_is_genuine ? h_norms.mean()
-							  : 0.0;
+	if (horizon_is_genuine) {
+          ps_ptr->ghosted_gridfn_norms(gfns::gfn__h, h_norms);
+        }
+	const fp mean_horizon_radius
+		= horizon_is_genuine ? h_norms.mean() : 0.0;
 	const bool horizon_is_too_large
 		= (mean_horizon_radius > solver_info
 					 .max_allowable_horizon_radius[hn]);
@@ -279,7 +702,7 @@ if (verbose_info.print_physics_details)
 
 	// does *this* horizon need more iterations?
 	// i.e. has this horizon's Newton iteration not yet converged?
-	const bool this_horizon_needs_more_iterations
+        const bool this_horizon_needs_more_iterations
 	   = horizon_is_genuine && Theta_is_ok
 	     && !found_this_horizon
 	     && !expansion_is_too_large
@@ -290,7 +713,8 @@ if (verbose_info.print_physics_details)
 	// on this or a following horizon?
 	const bool I_need_more_iterations
 	   = this_horizon_needs_more_iterations
-	     || there_is_another_genuine_horizon;
+	     || there_is_another_genuine_horizon
+             || I_am_pretracking;
 
 	if (verbose_info.print_algorithm_details)
 	   then {
@@ -305,7 +729,6 @@ if (verbose_info.print_physics_details)
 			   int(I_need_more_iterations));
 		}
 
-
 	//
 	// broadcast iteration status from each active processor
 	// to all processors, and inclusive-or the "we need more iterations"
@@ -373,6 +796,7 @@ if (verbose_info.print_physics_details)
 		broadcast_horizon_data(GH,
 				       my_proc == found_proc,
 				       broadcast_horizon_shape,
+                                       found_AH_data,
 				       found_AH_data.BH_diagnostics,
 				       *found_AH_data.ps_ptr,
 				       found_AH_data.horizon_buffers);
@@ -387,8 +811,9 @@ if (verbose_info.print_physics_details)
 	//
 	// if we found our horizon, maybe output the horizon shape?
 	//
-	if (found_this_horizon)
+	if (found_this_horizon && ! I_am_pretracking)
 	   then {
+		// printf("will output h/Th/mc: %d/%d/%d\n", IO_info.output_h, IO_info.output_Theta, IO_info.output_mean_curvature); //xxxxxxxxxxxx
 		if (IO_info.output_h)
 		   then {
 			// if this is the first time we've output h for this
@@ -396,13 +821,24 @@ if (verbose_info.print_physics_details)
 			if (!AH_data_ptr->h_files_written)
 			   then setup_h_files(*ps_ptr, IO_info, hn);
 			output_gridfn(*ps_ptr, gfns::gfn__h,
+                                      "h", GH,
 				      IO_info, IO_info.h_base_file_name,
+                                      IO_info.h_min_digits,
 				      hn, verbose_info
 					  .print_algorithm_highlights);
 			}
 		if (IO_info.output_Theta)
 		   then output_gridfn(*ps_ptr, gfns::gfn__Theta,
+                                      "Theta", GH,
 				      IO_info, IO_info.Theta_base_file_name,
+                                      IO_info.h_min_digits,
+				      hn, verbose_info
+					  .print_algorithm_highlights);
+		if (IO_info.output_mean_curvature)
+		   then output_gridfn(*ps_ptr, gfns::gfn__mean_curvature,
+                                      "mean_curvature", GH,
+				      IO_info, IO_info.mean_curvature_base_file_name,
+                                      IO_info.h_min_digits,
 				      hn, verbose_info
 					  .print_algorithm_highlights);
 		}
@@ -417,7 +853,8 @@ if (verbose_info.print_physics_details)
 		if (verbose_info.print_algorithm_details)
 		   then CCTK_VInfo(CCTK_THORNSTRING,
 				   "==> all processors are finished!");
-		return;					// *** NORMAL RETURN ***
+                do_return = true;
+		break;					// *** LOOP EXIT ***
 		}
 	if ((N_procs == 1) && !this_horizon_needs_more_iterations)
 	   then {
@@ -441,7 +878,7 @@ if (verbose_info.print_physics_details)
 	     = expansion_Jacobian
 		 (this_horizon_needs_more_iterations ? ps_ptr : NULL,
 		  this_horizon_needs_more_iterations ? Jac_ptr : NULL,
-		  add_to_expansion,
+		  compute_info,
 		  cgi, gi, Jacobian_info,
 		  error_info, (iteration == 1),
 		  verbose_info.print_algorithm_details);
@@ -474,7 +911,7 @@ if (verbose_info.print_physics_details)
 					  verbose_info.print_algorithm_details);
 	if	((rcond >= 0.0) && (rcond < 100.0*FP_EPSILON))
 	   then {
-		CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "Newton_solve: Jacobian matrix is numerically singular!");
 		// give up on this horizon
 		break;				// *** LOOP CONTROL ***
@@ -491,7 +928,9 @@ if (verbose_info.print_physics_details)
 
 	if (solver_info.debugging_output_at_each_Newton_iteration)
 	   then output_gridfn(*ps_ptr, gfns::gfn__Delta_h,
+                              "Delta_h", GH,
 			      IO_info, IO_info.Delta_h_base_file_name,
+                              IO_info.h_min_digits,
 			      hn, verbose_info.print_algorithm_details,
 			      iteration);
 
@@ -507,6 +946,139 @@ if (verbose_info.print_physics_details)
 	// end of this Newton iteration
 	}
 
+        if (! I_am_pretracking) {
+          if (horizon_is_genuine) {
+            if (! AH_data_ptr->initial_guess_info.reset_horizon_after_not_finding) {
+              if (! found_this_horizon) {
+                // the surface failed; backtrack and continue
+                AH_data_ptr->ps_ptr->ghosted_gridfn_copy
+                  (gfns::gfn__save_h, gfns::gfn__h);
+              }
+            }
+          }
+        }
+
+        // exit
+        if (do_return) return;				// *** NORMAL RETURN ***
+
+        // break out of the pretracking loop if we are not pretracking
+        if (! I_am_pretracking) break;
+        
+        if (! found_this_horizon) {
+          // the surface failed; backtrack and continue
+          CCTK_VInfo(CCTK_THORNSTRING,
+                     "Pretracking: solving failed; backtracking");
+          AH_data_ptr->ps_ptr->ghosted_gridfn_copy (gfns::gfn__save_h, gfns::gfn__h);
+          if (pretracking_have_lower_bound) {
+            assert (AH_data_ptr->pretracking_value >= pretracking_lower_bound - 1.0e-10 * fabs(pretracking_lower_bound));
+          }
+          pretracking_have_lower_bound = true;
+          pretracking_lower_bound = AH_data_ptr->pretracking_value;
+          if (pretracking_have_upper_bound) {
+            AH_data_ptr->pretracking_delta = 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+            AH_data_ptr->pretracking_value = pretracking_lower_bound + 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+          } else {
+            AH_data_ptr->pretracking_delta = fabs(AH_data_ptr->pretracking_delta);
+            AH_data_ptr->pretracking_delta *= 2.0;
+            if (AH_data_ptr->pretracking_delta > AH_data_ptr->pretracking_maximum_delta) {
+              AH_data_ptr->pretracking_delta = AH_data_ptr->pretracking_maximum_delta;
+            }
+            AH_data_ptr->pretracking_value += AH_data_ptr->pretracking_delta;
+            if (AH_data_ptr->pretracking_value > AH_data_ptr->pretracking_maximum_value) {
+              AH_data_ptr->pretracking_value = AH_data_ptr->pretracking_maximum_value;
+            }
+          }
+          CCTK_VInfo(CCTK_THORNSTRING,
+                     "Pretracking: new value %g, delta %g",
+                     double(AH_data_ptr->pretracking_value),
+                     double(AH_data_ptr->pretracking_delta));
+          if (pretracking_lower_bound >= AH_data_ptr->pretracking_maximum_value * 0.9999999999) {
+            // give up
+            CCTK_VInfo(CCTK_THORNSTRING,
+                       "Pretracking: upper bound reached, giving up.");
+            I_am_pretracking = false;
+            I_was_pretracking = true;
+            pretracking_was_successful = false;
+            // restore old pretracking goal
+            AH_data_ptr->pretracking_value = old_pretracking_value;
+          } else if (AH_data_ptr->pretracking_delta < AH_data_ptr->pretracking_minimum_delta) {
+            // give up
+            CCTK_VInfo(CCTK_THORNSTRING,
+                       "Pretracking: step size too small, giving up.");
+            I_am_pretracking = false;
+            I_was_pretracking = true;
+            pretracking_was_successful = true;
+          }
+        } else {
+          // the surface was okay
+          // get mean expansion
+          assert (pretracking_have_horizon_info);
+          CCTK_VInfo(CCTK_THORNSTRING,
+                     "Pretracking: solving succeeded; expansion is now %g",
+                     double(pretracking_mean_expansion));
+//           if (fabs(AH_data_ptr->pretracking_value) > 2.0 * AH_data_ptr->pretracking_minimum_delta) {
+          if (AH_data_ptr->pretracking_value > AH_data_ptr->pretracking_minimum_value + 1.0e-10 * AH_data_ptr->pretracking_minimum_delta) {
+            // it is not yet a horizon
+            if (pretracking_have_upper_bound) {
+              assert (AH_data_ptr->pretracking_value <= pretracking_upper_bound + 1.0e-10 * fabs(pretracking_upper_bound));
+            }
+            pretracking_have_upper_bound = true;
+            pretracking_upper_bound = AH_data_ptr->pretracking_value;
+            if (pretracking_have_lower_bound) {
+#if 1
+              // TODO
+              // move lower bound further down
+              pretracking_lower_bound -= pretracking_upper_bound - pretracking_lower_bound;
+              if (pretracking_lower_bound < AH_data_ptr->pretracking_minimum_value) pretracking_lower_bound = AH_data_ptr->pretracking_minimum_value;
+#endif
+              AH_data_ptr->pretracking_delta = 0.5 * (pretracking_lower_bound - pretracking_upper_bound);
+              AH_data_ptr->pretracking_value = pretracking_lower_bound + 0.5 * (pretracking_upper_bound - pretracking_lower_bound);
+            } else {
+              AH_data_ptr->pretracking_delta = - fabs(AH_data_ptr->pretracking_delta);
+              AH_data_ptr->pretracking_delta *= 2.0;
+              if (- AH_data_ptr->pretracking_delta > AH_data_ptr->pretracking_maximum_delta) {
+                AH_data_ptr->pretracking_delta = - AH_data_ptr->pretracking_maximum_delta;
+              }
+              AH_data_ptr->pretracking_value += AH_data_ptr->pretracking_delta;
+              if (AH_data_ptr->pretracking_value < AH_data_ptr->pretracking_minimum_value) {
+                AH_data_ptr->pretracking_value = AH_data_ptr->pretracking_minimum_value;
+              }
+            }
+            CCTK_VInfo(CCTK_THORNSTRING,
+                       "Pretracking: new value radius %g, delta %g",
+                       double(AH_data_ptr->pretracking_value),
+                       double(AH_data_ptr->pretracking_delta));
+            if (pretracking_upper_bound <= AH_data_ptr->pretracking_minimum_value * 1.00000000001) {
+              // give up
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: lower bound reached, giving up.");
+              I_am_pretracking = false;
+              I_was_pretracking = true;
+              pretracking_was_successful = false;
+              // restore old pretracking goal
+              AH_data_ptr->pretracking_value = old_pretracking_value;
+            } else if (- AH_data_ptr->pretracking_delta < AH_data_ptr->pretracking_minimum_delta) {
+              // give up
+              CCTK_VInfo(CCTK_THORNSTRING,
+                         "Pretracking: step size too small, giving up.");
+              I_am_pretracking = false;
+              I_was_pretracking = true;
+              pretracking_was_successful = true;
+            }
+          } else {
+            // a true horizon was found; we are done
+            CCTK_VInfo(CCTK_THORNSTRING,
+                       "Pretracking: done.");
+            I_am_pretracking = false;
+            I_was_pretracking = true;
+            pretracking_was_successful = true;
+          }
+        } // if surface found
+
+    // end of pretracking loop
+    }
+    I_am_pretracking = false;
+
     // end of this horizon
     }
 
@@ -703,8 +1275,8 @@ return any_proc_needs_more_iterations;
 
 //
 // This function (which must be called on *every* processor) broadcasts
-// the BH diagnostics and optionally also the (ghosted) horizon shape,
-// from a specified processor to all processors.
+// the BH diagnostics and (ghosted) horizon shape from a specified processor
+// to all processors.
 //
 // The present implementation of this function uses the Cactus reduction
 // API.  If AHFinderDirect is ported to some other software environment,
@@ -715,15 +1287,11 @@ return any_proc_needs_more_iterations;
 // GH --> The Cactus grid hierarchy.
 // broadcast_flag = true on the broadcasting processor,
 //		    false on all other processors
-// broadcast_horizon_shape = true to broadcast the (ghosted) horizon shape
-//				  as well as the BH diagnostics
-//			     false to only broadcast the BH diagnostics
 // BH_diagnostics = On the broadcasting processor, this is the BH diagnostics
 //		    to broadcast; on all other processors, it's set to the
 //		    broadcast BH diagnostics.
-// ps = Used only if broadcast_horizon_shape; in this case...
-//	On the broadcasting processor,  gfn__h  is broadcast;
-//	on all other processors,  gfn__h  is set to the broadcast values.
+// ps = On the broadcasting processor,  gfn__h  is broadcast;
+//      on all other processors,  gfn__h  is set to the broadcast values.
 // horizon_buffers = Internal buffers for use in the broadcast;
 //		     if  N_buffer == 0  then we set N_buffer and allocate
 //		     the buffers.
@@ -731,6 +1299,7 @@ return any_proc_needs_more_iterations;
 namespace {
 void broadcast_horizon_data(const cGH* GH,
 			    bool broadcast_flag, bool broadcast_horizon_shape,
+                            struct AH_data& AH_data,
 			    struct BH_diagnostics& BH_diagnostics,
 			    patch_system& ps,
 			    struct horizon_buffers& horizon_buffers)
@@ -751,7 +1320,8 @@ if (horizon_buffers.N_buffer == 0)
 	// allocate the buffers
 	horizon_buffers.N_buffer
 		= BH_diagnostics::N_buffer
-		  + (broadcast_horizon_shape ? ps.N_grid_points() : 0);
+		  + (broadcast_horizon_shape ? ps.N_grid_points() : 0)
+                  + 4;
 	horizon_buffers.send_buffer
 		= new CCTK_REAL[horizon_buffers.N_buffer];
 	horizon_buffers.receive_buffer
@@ -762,9 +1332,9 @@ if (broadcast_flag)
    then {
 	// pack the data to be broadcast into the send buffer
 	BH_diagnostics.copy_to_buffer(horizon_buffers.send_buffer);
+        int posn = BH_diagnostics::N_buffer;
 	if (broadcast_horizon_shape)
 	   then {
-		int posn = BH_diagnostics::N_buffer;
 			for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 			{
 			const patch& p = ps.ith_patch(pn);
@@ -782,8 +1352,12 @@ if (broadcast_flag)
 				}
 				}
 			}
-		assert( posn == horizon_buffers.N_buffer );
 		}
+        horizon_buffers.send_buffer[posn++] = AH_data.initial_find_flag;
+        horizon_buffers.send_buffer[posn++] = AH_data.really_initial_find_flag;
+        horizon_buffers.send_buffer[posn++] = AH_data.search_flag;
+        horizon_buffers.send_buffer[posn++] = AH_data.found_flag;
+        assert( posn == horizon_buffers.N_buffer );
 	}
    else jtutil::zero_C_array(horizon_buffers.N_buffer,
 			     horizon_buffers.send_buffer);
@@ -814,9 +1388,12 @@ if (!broadcast_flag)
    then {
 	// unpack the data from the receive buffer
 	BH_diagnostics.copy_from_buffer(horizon_buffers.receive_buffer);
+        ps.origin_x(BH_diagnostics.origin_x);
+        ps.origin_y(BH_diagnostics.origin_y);
+        ps.origin_z(BH_diagnostics.origin_z);
+	int posn = BH_diagnostics::N_buffer;
 	if (broadcast_horizon_shape)
 	   then {
-		int posn = BH_diagnostics::N_buffer;
 			for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 			{
 			patch& p = ps.ith_patch(pn);
@@ -834,12 +1411,16 @@ if (!broadcast_flag)
 				}
 				}
 			}
-		assert( posn == horizon_buffers.N_buffer );
 
 		// recover the full ghosted-grid horizon shape
 		// (we only broadcast the nominal-grid shape)
 		ps.synchronize();
 		}
+        AH_data.initial_find_flag        = horizon_buffers.receive_buffer[posn++];
+        AH_data.really_initial_find_flag = horizon_buffers.receive_buffer[posn++];
+        AH_data.search_flag              = horizon_buffers.receive_buffer[posn++];
+        AH_data.found_flag               = horizon_buffers.receive_buffer[posn++];
+	assert( posn == horizon_buffers.N_buffer );
 	}
 }
 	  }
@@ -918,6 +1499,7 @@ const fp scale = (max_Delta_h <= max_allowable_Delta_h)
 
 if (verbose_info.print_algorithm_details)
    then {
+
 	if (scale == 1.0)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "h += Delta_h (rms-norm=%.1e, infinity-norm=%.1e)",
@@ -948,6 +1530,14 @@ if (verbose_info.print_algorithm_details)
 		}
 		}
 	}
+
+if (ps.N_additional_points())
+	{
+	const int np = ps.N_grid_points();
+	const int gnp = ps.ghosted_N_grid_points();
+	ps.ghosted_gridfn_data(gfns::gfn__h)[gnp]
+		-= scale * ps.gridfn_data(gfns::gfn__Delta_h)[np];
+	}
 }
 	  }
 
diff --git a/src/driver/README b/src/driver/README
index f3de3e7..a19bbb0 100644
--- a/src/driver/README
+++ b/src/driver/README
@@ -16,24 +16,16 @@ state.cc
 setup.cc		# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
 	this is called from the scheduler to setup our data structures
 
-find_horizons.cc	# sees CCTK_ARGUMENTS
+find_horizons.cc	# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
 	this is called from the scheduler to find the apparent horizon(s)
 	in a slice
 
-mask.cc			# sees CCTK_ARGUMENTS
-	this is called from the scheduler to set an excision mask or masks
-
-announce.cc		# sees CCTK_ARGUMENTS
+announce.cc		# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
 	this is called from the scheduler to announce apparent horizon
-	info to other thorns via aliased function calls
+	info to other thorns
 
-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
-	interface, for the other thorns to find out about apparent horizons
+mask.cc			# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
+	this is called from the scheduler to set an excision mask or masks
 
 initial_guess.cc
 	this sets up the initial guess(es) for the horizon position(s)
@@ -46,24 +38,6 @@ Newton.cc		# sees cGH (for interprocessor broadcasts)
 io.cc
 	I/O routines
 
-BH_diagnostics.{cc,hh}
-	this struct (= class with public data) computes and keeps track of
-	the BH diagnostics (horizon circumferences, area, irreducible mass,
-	etc)
-
-misc-driver.cc
-	some misc code that doesn't fit anywhere else :)
-
-Makefile.standalone
-	Makefile for standalone test drivers etc
-
-horizon_sequence.{cc,hh}
-	this class keeps track of the sequence of horizons which a given
-	processor should try to find (possibly including dummy ones to
-	preserve the multiprocessor synchronization)
-test_horizon_sequence.cc
-	standalone test driver for  horizon_sequence.{cc,hh}
-
 ellipsoid.maple
 	this is a Maple script to compute the intersection of a given
 	ray with an ellipsoid; this is used in setting up the initial guess
diff --git a/src/driver/README.parallel b/src/driver/README.parallel
index 67e3394..136c026 100644
--- a/src/driver/README.parallel
+++ b/src/driver/README.parallel
@@ -195,11 +195,3 @@ received from the broadcast, and for each flag which is true, all
 processors then know to participate in the broadcast of the BH diagnostics
 and (optionally) the horizon shape from the just-found-it processor
 to all processors.
-
-We only broadcast the horizon shapes if they're needed, i.e. if either
-or both of the following conditions hold:
-* we're setting an excision mask
-* we're storing horizon information in the SphericalSurface variables
-  for at least one horizon
-In either case, each processor needs to know the shapes of *all* horizons,
-so we do the broadcast.
diff --git a/src/driver/aliased_functions.cc b/src/driver/aliased_functions.cc
index 43c15ed..8b1218a 100644
--- a/src/driver/aliased_functions.cc
+++ b/src/driver/aliased_functions.cc
@@ -12,10 +12,12 @@
 #include <assert.h>
 #include <math.h>
 
+#include <vector>
+
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
-#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -23,6 +25,7 @@
 #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"
@@ -45,7 +48,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -57,13 +59,8 @@ extern struct state state;
 //******************************************************************************
 
 //
-// This function is called (via the Cactus flesh function-aliasing mechanism)
-// by other thorns to find out our local coordinate origin for a given AH.
-//
-// Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-//			 1 to N_horizons inclusive.
-// *p_origin_[xyz] = (out) The local coordinate origin is stored here.
+// This function is called (via the magic of function aliasing) by
+// other thorns to find out our local coordinate origin for a given AH.
 //
 // Results:
 // This function returns 0 for ok, or -1 if the horizon number is invalid.
@@ -77,7 +74,7 @@ const struct verbose_info& verbose_info = state.verbose_info;
 
 if (!  ((horizon_number >= 1) && (horizon_number <= state.N_horizons))  )
    then {
-	CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_local_coordinate_origin():\n"
 "        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
 		   ,
@@ -108,10 +105,6 @@ return 0;						// *** NORMAL RETURN ***
 // by other thorns to query whether or not the specified horizon was found
 // the last time we searched for it.
 //
-// Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-//			 1 to N_horizons inclusive.
-//
 // Results:
 // This function returns
 //  1 if the horizon was found
@@ -125,7 +118,7 @@ const struct verbose_info& verbose_info = state.verbose_info;
 
 if (!  ((horizon_number >= 1) && (horizon_number <= state.N_horizons))  )
    then {
-	CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_horizon_was_found():\n"
 "        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
 		   ,
@@ -202,26 +195,16 @@ if (AH_data.found_flag)
 // its local coordinate origin to a given (x,y,z) coordinate or coordinates.
 //
 // Arguments:
-// horizon_number = (in) The horizon to inquire about; must be in the range
-//			 1 to N_horizons inclusive.
-// N_points = (in) Number of (x,y,z) coordinate tuples (should be >= 0).
-// x[], y[], z[] = (in) These arrays give the (x,y,z) coordinates.
-// radius[] = (out) This array is set to the horizon radius values
-//		    (Euclidean distances from the local coordinate origin),
-//		    or to all -1.0 if we didn't find this horizon the
-//		    last time we looked for it.
+// horizon_number = must be in the range 1 to N_horizons
+// N_points = should be >= 0
+// x[], y[], z[] = these give the (x,y,z) coordinates
+// radius[] = this is set to the horizon radius values (Euclidean distance
+//	      from the local coordinate origin), or to all -1.0 if we didn't
+//	      find this horizon the last time we looked for it
 //
 // Results:
 // This function returns 0 for ok, or -1 if the horizon number is invalid.
 //
-// FIXME:
-// * The current implementation is quite inefficient: it does a full 2-D
-//   interpolation (to find the horizon radius) for each of the xyz points.
-//   It would be more efficient to batch the interpolations.
-// * If an (x,y,z) coordinate tuple coincides with this horizon's local
-//   coordinate origin, we return the dummy radius 1.0.  It would be
-//   better to return -1.0.
-//
 extern "C"
   CCTK_INT AHFinderDirect_radius_in_direction
     (CCTK_INT horizon_number,
@@ -229,37 +212,52 @@ extern "C"
      const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const z,
      CCTK_REAL* const radius)
 {
-const struct verbose_info& verbose_info = state.verbose_info;
-
-if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) )
-   then {
-	CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
-"AHFinderDirect_distance_outside_thorn():\n"
-"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
-		   ,
-		   int(horizon_number), state.N_horizons);
-	return -1;					// *** ERROR RETURN ***
-	}
-
-assert(state.AH_data_array[horizon_number] != NULL);
-const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
-
-	for (int point = 0 ; point < N_points ; ++point)
-	{
-	if (AH_data.found_flag)
-	   then {
-		assert(AH_data.ps_ptr != NULL);
-		const patch_system& ps = *AH_data.ps_ptr;
-
-		const fp local_x = x[point] - ps.origin_x();
-		const fp local_y = y[point] - ps.origin_y();
-		const fp local_z = z[point] - ps.origin_z();
-		radius[point] = ps.radius_in_local_xyz_direction
-				      (gfns::gfn__h,
-				       local_x, local_y, local_z);
-		}
-	   else radius[point] = -1.0;
-	}
+  const struct verbose_info& verbose_info = state.verbose_info;
+  
+  if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) ) {
+    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "AHFinderDirect_distance_outside_thorn():\n"
+               "        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
+               ,
+               int(horizon_number), state.N_horizons);
+    return -1;					// *** ERROR RETURN ***
+  }
+  
+  assert(state.AH_data_array[horizon_number] != NULL);
+  const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+  
+  if (AH_data.found_flag) {
+    
+    assert(AH_data.ps_ptr != NULL);
+    const patch_system& ps = *AH_data.ps_ptr;
+    
+    std::vector<fp> local_xs(N_points);
+    std::vector<fp> local_ys(N_points);
+    std::vector<fp> local_zs(N_points);
+    
+    for (int point = 0 ; point < N_points ; ++point) {
+      
+      local_xs.at(point) = x[point] - ps.origin_x();
+      local_ys.at(point) = y[point] - ps.origin_y();
+      local_zs.at(point) = z[point] - ps.origin_z();
+      
+    }
+    
+    ps.radii_in_local_xyz_directions (gfns::gfn__h,
+                                      N_points,
+                                      & local_xs.front(),
+                                      & local_ys.front(),
+                                      & local_zs.front(),
+                                      radius);
+    
+  } else {
+    // if not found
+    
+    for (int point = 0 ; point < N_points ; ++point) {
+      radius[point] = -1.0;
+    }
+    
+  } // if not found
 
 return 0;						// *** NORMAL RETURN ***
 }
diff --git a/src/driver/announce.cc b/src/driver/announce.cc
index 54b8972..5c0b62b 100644
--- a/src/driver/announce.cc
+++ b/src/driver/announce.cc
@@ -12,6 +12,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -19,6 +20,7 @@
 #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"
@@ -41,7 +43,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -63,35 +64,42 @@ extern "C"
   void AHFinderDirect_announce(CCTK_ARGUMENTS)
 {
 DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
 
 const struct verbose_info& verbose_info = state.verbose_info;
 
-// only try to announce AH info if we've found AHs at this time level
-if (! state.find_now(cctk_iteration))
+// which horizon to announce?
+const int hn = which_horizon_to_announce_centroid;
+if (hn == 0)
    then return;						// *** NO-OP RETURN ***
 
-// only try to anounce AH info if we've been asked to do so
-if (! state.announce_centroid_flag)
-   then return;						// *** NO-OP RETURN ***
+if (! ((hn >= 1) && (hn <= N_horizons)) )
+   then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   AHFinderDirect_announce():\n"
+"        invalid horizon number %d to announce\n"
+"        (valid range is [1,N_horizons=%d])!\n"
+		   ,
+		   hn, int(N_horizons));			/*NOTREACHED*/
 
-// which horizon to announce?
-const int hn = state.which_horizon_to_announce_centroid;
+assert(state.AH_data_array[hn] != NULL);
+const struct AH_data& AH_data = *state.AH_data_array[hn];
+
+// only try to announce AH info if we've found AHs at this time level
+if (! AH_data.search_flag)
+   then return;						// *** NO-OP RETURN ***
 
 // did we actually *find* this horizon?
-if (state.AH_data_array[hn] == NULL)
+if (! AH_data.found_flag)
    then return;						// *** NO-OP RETURN ***
 
 // 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",
@@ -102,4 +110,89 @@ if (CCTK_IsFunctionAliased("SetDriftCorrectPosition"))
 
 //******************************************************************************
 
+//
+// This function is called by the Cactus scheduler, to copy any
+// desired apparent horizon info to Cactus variables.
+//
+extern "C"
+  void AHFinderDirect_store(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+for (int hn = 1; hn <= N_horizons; ++ hn)
+  {
+
+  // Store in spherical surface
+  const int sn = which_surface_to_store_info[hn];
+  if (sn == -1)
+    then continue;
+
+  if (sn < 0 || sn >= nsurfaces)
+    then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   AHFinderDirect_store():\n"
+"        invalid surface number %d for horizon number %d\n"
+"        (valid range is [0,nsurfaces-1=%d])!\n"
+		   ,
+		   sn, hn,
+                   int(nsurfaces-1));			/*NOTREACHED*/
+  
+  const struct AH_data& AH_data = *state.AH_data_array[hn];
+  const struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
+  BH_diagnostics.store(cctkGH, hn, sn);
+
+  }
+}
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler, to copy any
+// desired apparent horizon info to Cactus variables.
+//
+extern "C"
+  void AHFinderDirect_save(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+for (int hn = 1; hn <= N_horizons; ++ hn)
+  {
+
+  const struct AH_data& AH_data = *state.AH_data_array[hn];
+  const struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
+
+  // Save in grid array
+  BH_diagnostics.save(cctkGH, hn);
+
+  }
+}
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler, to copy any
+// desired apparent horizon info from Cactus variables.
+//
+extern "C"
+  void AHFinderDirect_recover(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+for (int hn = 1; hn <= N_horizons; ++ hn)
+  {
+
+  struct AH_data& AH_data = *state.AH_data_array[hn];
+  struct BH_diagnostics& BH_diagnostics = AH_data.BH_diagnostics;
+
+  // Load from grid array
+  BH_diagnostics.load(cctkGH, hn);
+
+  }
+}
+
+//******************************************************************************
+
 	  }	// namespace AHFinderDirect
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index ae169fd..6569ba8 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -34,7 +34,6 @@ enum	method
 //
 enum	verbose_level
 	{
-	verbose_level__no_output,
 	verbose_level__physics_highlights,
 	verbose_level__physics_details,
 	verbose_level__algorithm_highlights,
@@ -57,17 +56,6 @@ enum	initial_guess_method
 	initial_guess__coord_ellipsoid // no comma
 	};
 
-//
-// this enum holds the decoded  horizon_file_format  parameter, i.e.
-// it specifies what format of input/output file(s) we should use
-// for h and H (and other angular grid functions)
-//
-enum	horizon_file_format
-	{
-	horizon_file_format__ASCII_gnuplot,
-	horizon_file_format__HDF5 // no comma
-	};
-
 //******************************************************************************
 
 //
@@ -77,6 +65,7 @@ enum	horizon_file_format
 struct	initial_guess_info
 	{
 	enum initial_guess_method method;
+        bool reset_horizon_after_not_finding;
 
 	// parameters for method == initial_guess__read_from_named_file
 	struct	{
@@ -120,9 +109,12 @@ struct	solver_info
 	fp max_allowable_Delta_h_over_h;
 	fp Theta_norm_for_convergence;
 	fp max_allowable_Theta;
+        fp max_allowable_Theta_growth_iterations;
+        fp max_allowable_Theta_nonshrink_iterations;
 	fp *max_allowable_horizon_radius;	// --> new[]-allocated array
 						//     of size  N_horizons+1 ,
 						//     subscripted by hn
+	bool want_expansion_gradients;
 	};
 
 //
@@ -140,11 +132,13 @@ struct	IO_info
 	//   should cut this down to something reasonable
 	enum { default_directory_permission = 0777 };
 
-	enum horizon_file_format horizon_file_format;
+        bool output_ASCII_files;
+        bool output_HDF5_files;
 	bool output_initial_guess;
 	int output_h_every, output_Theta_every;
+	int output_mean_curvature_every;
 	// based on the above, do we want to output things now (this time step)?
-	bool output_h, output_Theta;
+	bool output_h, output_Theta, output_mean_curvature;
 
 	bool output_BH_diagnostics;
 	const char* BH_diagnostics_directory;
@@ -157,7 +151,9 @@ struct	IO_info
 	const char* h_directory;
 	const char* h_base_file_name;
 	const char* Theta_base_file_name;
+	const char* mean_curvature_base_file_name;
 	const char* Delta_h_base_file_name;
+        int h_min_digits;
 
 	const char* Jacobian_base_file_name;
 
@@ -289,7 +285,26 @@ struct	AH_data
 	//
 	patch_system* ps_ptr;
 	Jacobian* Jac_ptr;
-	fp surface_expansion;
+        what_to_compute compute_info;
+
+        bool move_origins;
+
+        bool use_pretracking;
+        int pretracking_max_iterations;
+
+        fp pretracking_value;
+        fp pretracking_minimum_value;
+        fp pretracking_maximum_value;
+        fp pretracking_delta;
+        fp pretracking_minimum_delta;
+        fp pretracking_maximum_delta;
+
+        int depends_on;
+        fp desired_value_factor;
+        fp desired_value_offset;
+
+        fp shiftout_factor;
+        fp smoothing_factor;
 
 	// are we finding this horizon "from scratch"?
 	// ... true if this is the first time we've tried to find it,
@@ -299,10 +314,16 @@ struct	AH_data
 	//           (so we have that position as a very good initial guess)
 	// ... also false if we're not finding this horizon on this processor
 	bool initial_find_flag;
+        // is this really the first time in this simulation that we are
+        // trying to find a horizon?
+        // ... true initially if this is a genuine horizon,
+        //     false after the first time that a horizon has been found
+	bool really_initial_find_flag;
 
 	// used only if we're finding this horizon on this processor
 	struct initial_guess_info initial_guess_info;
 
+	bool search_flag;	// did we search for this horizon
 	bool found_flag;	// did we find this horizon (successfully)
 	bool h_files_written;	// have we written horizon-shape or similar
 				// files for this horizon yet?
@@ -310,13 +331,6 @@ 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;
@@ -329,20 +343,9 @@ struct	AH_data
 //
 struct	state
 	{
-	int find_every;
-	bool find_now(int cctk_iteration_in) const
-		{
-		return (find_every != 0)
-		       && ((cctk_iteration_in % find_every) == 0);
-		}
-
 	enum method method;
-
 	struct error_info error_info;
 	struct verbose_info verbose_info;
-
-	const char* metric_type;	// copy of ADMBase::metric_type
-
 	int timer_handle;
 	int N_procs;			// total number of processors
 	int my_proc;			// processor number of this processor
@@ -356,7 +359,6 @@ struct	state
 
 	struct cactus_grid_info cgi;
 	struct geometry_info gi;
-	bool test_all_Jacobian_compute_methods;
 	struct Jacobian_info Jac_info;
 	struct solver_info solver_info;
 	struct IO_info IO_info;
@@ -364,19 +366,6 @@ struct	state
 	struct BH_diagnostics_info BH_diagnostics_info;
 	struct mask_info mask_info;
 
-	bool announce_centroid_flag;	// should we announce horizon centroid?
-	int which_horizon_to_announce_centroid;	// if so, which horizon?
-
-	// should we always broadcast each horizon shape from the processor
-	// that finds it, to all processors, even if there's no apparent need
-	// for this information on other processors?
-	bool always_broadcast_horizon_shape;
-
-	// this is the inclusive-or of  AH_data.store_info_in_SS_vars
-	// over all horizons, i.e. it says if we will store horizon
-	// information in the SphericalSurface variables for *any* horizon
-	bool store_info_in_SS_vars_for_any_horizon;
-
 	// interprocessor-communication buffers for broadcasting
 	// Newton-iteration status from active processors to all processors
 	struct iteration_status_buffers isb;
@@ -402,10 +391,6 @@ struct	state
 extern "C"
   void AHFinderDirect_setup(CCTK_ARGUMENTS);
 
-// ... called from Cactus Scheduler
-extern "C"
-  void AHFinderDirect_setupupdate(CCTK_ARGUMENTS);
-
 // find_horizons.cc
 // ... called from Cactus Scheduler
 extern "C"
@@ -416,11 +401,6 @@ 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"
@@ -438,7 +418,8 @@ extern "C"
   CCTK_INT AHFinderDirect_radius_in_direction
     (CCTK_INT horizon_number,
      CCTK_INT N_points,
-     const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const z,
+     const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const 
+z,
      CCTK_REAL* const radius);
 
 // initial_guess.cc
@@ -469,18 +450,22 @@ void Newton(const cGH* GH,
 // io.cc
 void input_gridfn(patch_system& ps, int unknown_gfn,
 		  const struct IO_info& IO_info, const char base_file_name[],
+		  int min_digits,
 		  int hn, bool print_msg_flag, int AHF_iteration = 0);
 void input_gridfn__explicit_name(patch_system& ps, int unknown_gfn,
 				 const struct IO_info& IO_info,
 				 const char file_name[], bool print_msg_flag);
 void setup_h_files(patch_system& ps, const struct IO_info& IO_info, int hn);
 void output_gridfn(patch_system& ps, int unknown_gfn,
+                   const char gfn_name[], const cGH *cctkGH,
 		   const struct IO_info& IO_info, const char base_file_name[],
+                   int min_digits,
 		   int hn, bool print_msg_flag, int AHF_iteration = 0);
 void output_Jacobians(const patch_system& ps,
 		      const Jacobian* Jac_NP_ptr,
 		      const Jacobian* Jac_SD_FDdr_ptr,
 		      const struct IO_info& IO_info, const char base_file_name[],
+                      int min_digits,
 		      int hn, bool print_msg_flag, int AHF_iteration = 0);
 
 // misc-driver.cc
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index bf3f25a..45cc780 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -6,8 +6,6 @@
 // AHFinderDirect_find_horizons - top-level driver to find apparent horizons
 ///
 /// find_horizon - find a horizon
-/// print_summary_of_which_horizons_found
-///
 /// do_evaluate_expansions
 /// do_test_expansion_Jacobian
 ///
@@ -15,10 +13,12 @@
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
+#include <string.h>
 
-#include "util_String.h"
+#include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -26,6 +26,7 @@
 #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"
@@ -48,7 +49,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -63,8 +63,6 @@ extern struct state state;
 // ***** prototypes for functions local to this file
 //
 namespace {
-void print_summary_of_which_horizons_found
-	(int N_horizons, const struct AH_data* const* AH_data_array);
 void do_evaluate_expansions(int my_proc, int N_horizons,
 			    horizon_sequence& hs,
 			       struct AH_data* const AH_data_array[],
@@ -89,6 +87,36 @@ void do_test_expansion_Jacobians(int my_proc, int N_horizons,
 //******************************************************************************
 
 //
+// This function is called by the Cactus scheduler to import
+// the excision mask.
+//
+extern "C"
+  void AHFinderDirect_import_mask(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+for (int k=0; k<cctk_lsh[2]; ++k)
+for (int j=0; j<cctk_lsh[1]; ++j)
+for (int i=0; i<cctk_lsh[0]; ++i)
+{
+        const int ind = CCTK_GFINDEX3D(cctkGH,i,j,k);
+        // zero means: point can be used,
+        // non-zero means: point must be avoided
+        ahmask[ind] = 0;
+        if (use_mask)
+           // grid points with mask values of 1.0 can be used,
+           // values of 0.0 and 0.5 must be avoided.
+           // the excision boundary cannot be used because
+           // (a) it is inaccurate
+           // (b) it does not respect the symmetries e.g. in Kerr.
+           then ahmask[ind] = fabs(emask[ind] - 1.0) > 0.01;
+}
+}
+
+//******************************************************************************
+
+//
 // This function is called by the Cactus scheduler to find the apparent
 // horizon or horizons in the current slice.
 //
@@ -96,22 +124,44 @@ extern "C"
   void AHFinderDirect_find_horizons(CCTK_ARGUMENTS)
 {
 DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
 
-// only try to find horizons every  find_every  time steps
-if (! state.find_now(cctk_iteration))
-   then return;						// *** NO-OP RETURN ***
+// determine whether a horizon should be found at this iteration
+bool find_any = false;
+for (int hn = 1; hn <= state.my_hs->N_horizons(); ++ hn)
+{
+  // only try to find horizons every  find_every  time steps
+  const int my_find_after = find_after_individual[hn];
+  const int my_dont_find_after = dont_find_after_individual[hn];
+  const fp my_find_after_time = find_after_individual_time[hn];
+  const fp my_dont_find_after_time = dont_find_after_individual_time[hn];
+  const int my_find_every = (find_every_individual[hn] >= 0
+                             ? find_every_individual[hn]
+                             : find_every);
+  const bool find_this =    cctk_iteration >= my_find_after
+                         && (my_dont_find_after < 0
+                             ? true
+                             : cctk_iteration <= my_dont_find_after)
+                         && cctk_time >= my_find_after_time
+                         && (my_dont_find_after_time <= my_find_after_time
+                             ? true
+                             : cctk_time <= my_dont_find_after_time)
+                         && my_find_every > 0
+                         && cctk_iteration % my_find_every == 0
+                         && ! disable_horizon[hn];
+  struct AH_data& AH_data = *state.AH_data_array[hn];
+  AH_data.search_flag = find_this;
+  find_any = find_any || find_this;
+}
+if (! find_any) return;
 
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
 
 const int my_proc = state.my_proc;
-const int N_horizons = state.N_horizons;
 horizon_sequence& hs = *state.my_hs;
 const bool active_flag = hs.has_genuine_horizons();
-const bool broadcast_horizon_shape
-	= state.always_broadcast_horizon_shape
-	  || state.mask_info.set_mask_for_any_horizon
-	  || state.store_info_in_SS_vars_for_any_horizon;
+const bool broadcast_horizon_shape = true;
 
       struct cactus_grid_info&          cgi = state.cgi;
 const struct    geometry_info&           gi = state.gi;
@@ -122,13 +172,50 @@ const struct     verbose_info& verbose_info = state.verbose_info;
 
 // what are the semantics of the Cactus gxx variables? (these may
 // change from one call to another, so we have to re-check each time)
-if      (CCTK_Equals(state.metric_type, "physical"))
+if      (CCTK_Equals(metric_type, "physical"))
    then cgi.use_Cactus_conformal_metric = false;
-else if (CCTK_Equals(state.metric_type, "static conformal"))
+else if (CCTK_Equals(metric_type, "static conformal"))
    then cgi.use_Cactus_conformal_metric = (*conformal_state > 0);
 else	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"AHFinderDirect_find_horizons(): unknown ADMBase::metric_type=\"%s\"!",
-		   state.metric_type);				/*NOTREACHED*/
+"AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
+		   metric_type);				/*NOTREACHED*/
+
+// update parameters
+IO_info.output_ASCII_files = (output_ASCII_files != 0);
+IO_info.output_HDF5_files = (output_HDF5_files != 0);
+IO_info.output_initial_guess = (output_initial_guess != 0);
+IO_info.output_h_every     = output_h_every;
+IO_info.output_Theta_every = output_Theta_every;
+IO_info.output_mean_curvature_every = output_mean_curvature_every;
+IO_info.output_h     = false;	// dummy value
+IO_info.output_Theta = false;	// dummy value
+IO_info.output_mean_curvature = false;	// dummy value
+
+IO_info.output_BH_diagnostics              = (output_BH_diagnostics != 0);
+IO_info.BH_diagnostics_directory
+	= (strlen(BH_diagnostics_directory) == 0)
+	  ? /* IO:: */ out_dir
+	  : BH_diagnostics_directory;
+IO_info.BH_diagnostics_base_file_name      = BH_diagnostics_base_file_name;
+IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension;
+
+IO_info.output_ghost_zones_for_h  = (output_ghost_zones_for_h != 0);
+IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension;
+IO_info.HDF5_file_name_extension          = HDF5_file_name_extension;
+IO_info.h_directory
+	= (strlen(h_directory) == 0)
+	  ? /* IO:: */ out_dir
+	  : h_directory;
+IO_info.h_base_file_name         = h_base_file_name;
+IO_info.Theta_base_file_name     = Theta_base_file_name;
+IO_info.mean_curvature_base_file_name     = mean_curvature_base_file_name;
+IO_info.Delta_h_base_file_name   = Delta_h_base_file_name;
+IO_info.h_min_digits             = h_min_digits;
+IO_info.Jacobian_base_file_name  = Jacobian_base_file_name;
+IO_info.output_OpenDX_control_files  = (output_OpenDX_control_files != 0);
+IO_info.OpenDX_control_file_name_extension = OpenDX_control_file_name_extension;
+IO_info.time_iteration = 0;
+IO_info.time           = 0.0;
 
 // get the Cactus time step and decide if we want to output h and/or Theta now
 IO_info.time_iteration = cctk_iteration;
@@ -139,6 +226,9 @@ IO_info.output_h
 IO_info.output_Theta
    = (IO_info.output_Theta_every > 0)
      && ((IO_info.time_iteration % IO_info.output_Theta_every) == 0);
+IO_info.output_mean_curvature
+   = (IO_info.output_mean_curvature_every > 0)
+     && ((IO_info.time_iteration % IO_info.output_mean_curvature_every) == 0);
 
 // set initial guess for any (genuine) horizons that need it,
 // i.e. for any (genuine) horizons where we didn't find the horizon previously
@@ -146,21 +236,39 @@ IO_info.output_Theta
 	{
 	assert( state.AH_data_array[hn] != NULL );
 	struct AH_data& AH_data = *state.AH_data_array[hn];
+        if (verbose_info.print_algorithm_highlights) {
+          printf ("AHF find_horizons[%d] initial_find_flag=%d\n", hn, (int) AH_data.initial_find_flag);
+          printf ("AHF find_horizons[%d] really_initial_find_flag=%d\n", hn, (int) AH_data.really_initial_find_flag);
+          printf ("AHF find_horizons[%d] search_flag=%d\n", hn, (int) AH_data.search_flag);
+          printf ("AHF find_horizons[%d] found_flag=%d\n", hn, (int) AH_data.found_flag);
+        }
 	if (AH_data.found_flag)
-	   then AH_data.initial_find_flag = false;
+           then {
+                AH_data.initial_find_flag = false;
+                AH_data.really_initial_find_flag = false;
+                }
 	   else {
-		patch_system& ps = *AH_data.ps_ptr;
-		setup_initial_guess(ps,
-				    AH_data.initial_guess_info,
-				    IO_info,
-				    hn, N_horizons, verbose_info);
-		if (active_flag && IO_info.output_initial_guess)
-		   then output_gridfn(ps, gfns::gfn__h,
-				      IO_info, IO_info.h_base_file_name,
-				      hn, verbose_info
-					  .print_algorithm_highlights);
-		AH_data.initial_find_flag = true;
-		}
+                if (AH_data.really_initial_find_flag
+                    || AH_data.initial_guess_info.reset_horizon_after_not_finding)
+                   then {
+		        patch_system& ps = *AH_data.ps_ptr;
+                        if (verbose_info.print_algorithm_highlights) {
+                          printf ("AHF find_horizons[%d] setup_initial_guess\n", hn);
+                        }
+        		setup_initial_guess(ps,
+        				    AH_data.initial_guess_info,
+        				    IO_info,
+        				    hn, N_horizons, verbose_info);
+        		if (active_flag && IO_info.output_initial_guess)
+        		   then output_gridfn(ps, gfns::gfn__h,
+                                              "h", cgi.GH,
+        				      IO_info, IO_info.h_base_file_name,
+                                              IO_info.h_min_digits,
+        				      hn, verbose_info
+        				      .print_algorithm_highlights);
+        		AH_data.initial_find_flag = true;
+        		}
+                }
 	}
 
 //
@@ -180,7 +288,7 @@ case method__test_expansion_Jacobians:
 	do_test_expansion_Jacobians(my_proc, N_horizons,
 				    state.AH_data_array,
 				    cgi, gi, Jac_info,
-				    state.test_all_Jacobian_compute_methods,
+				    (test_all_Jacobian_compute_methods != 0),
 				    IO_info, error_info, verbose_info,
 				    state.timer_handle);
 	break;
@@ -198,10 +306,6 @@ case method__find_horizons:
 	       state.isb);
 	if (state.timer_handle >= 0)
 	   then CCTK_TimerStopI(state.timer_handle);
-	if (verbose_info.print_physics_highlights
-	    && !verbose_info.print_physics_details)
-	   then print_summary_of_which_horizons_found(N_horizons,
-						      state.AH_data_array);
 	break;
 	  }
 
@@ -223,53 +327,6 @@ if (state.timer_handle >= 0)
 }
 
 //******************************************************************************
-
-//
-// This function prints (using CCTK_VInfo()) a 1-line summary of
-// which AHs were / were not found.
-//
-namespace {
-void print_summary_of_which_horizons_found
-	(int N_horizons, const struct AH_data* const* AH_data_array)
-{
-const int hn_buffer_size = 10;		// buffer for single "%d" for hn
-char hn_buffer[hn_buffer_size];
-const int msg_buffer_size = 200;	// buffer for the entire line
-char msg_buffer[msg_buffer_size];
-
-msg_buffer[0] = '\0';
-Util_Strlcpy(msg_buffer, "found horizon(s) ", msg_buffer_size);;
-
-bool already_found_flag = false;	// running inclusive-or of  found_flag
-	for (int hn = 1 ; hn <= N_horizons ; ++hn)
-	{
-	assert(AH_data_array[hn] != NULL);
-	bool found_flag = AH_data_array[hn]->found_flag;
-	if (found_flag)
-	   then {
-		if (already_found_flag)
-		   then Util_Strlcat(msg_buffer, ",", msg_buffer_size);
-		already_found_flag = true;
-		snprintf(hn_buffer, hn_buffer_size, "%d", hn);
-		Util_Strlcat(msg_buffer, hn_buffer, msg_buffer_size);
-		}
-	}
-
-if (already_found_flag)
-   then {
-	// we found one or more horizons
-	Util_Strlcat(msg_buffer, " of ", msg_buffer_size);
-
-	snprintf(hn_buffer, hn_buffer_size, "%d", N_horizons);
-	Util_Strlcat(msg_buffer, hn_buffer, msg_buffer_size);
-	}
-   else Util_Strlcpy(msg_buffer, "no horizons found", msg_buffer_size);
-
-CCTK_VInfo(CCTK_THORNSTRING, msg_buffer);
-}
-	  }
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -317,7 +374,8 @@ if (active_flag)
 		if (timer_handle >= 0)
 		   then CCTK_TimerStartI(timer_handle);
 		jtutil::norm<fp> Theta_norms;
-		const bool Theta_ok = expansion(&ps, -AH_data.surface_expansion,
+		const bool Theta_ok = expansion(&ps,
+                                                AH_data.compute_info,
 						cgi, gi,
 						error_info, true,// initial eval
 						false,	// no Jacobian coeffs
@@ -328,7 +386,9 @@ if (active_flag)
 
 		if (IO_info.output_h)
 		   then output_gridfn(ps, gfns::gfn__h,
+                                      "h", cgi.GH,
 				      IO_info, IO_info.h_base_file_name,
+                                      IO_info.h_min_digits,
 				      hn, verbose_info.print_algorithm_details);
 
 		if (Theta_ok)
@@ -338,17 +398,28 @@ if (active_flag)
 			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
 			if (IO_info.output_Theta)
 			   then output_gridfn(ps, gfns::gfn__Theta,
+                                              "Theta", cgi.GH,
 					      IO_info, IO_info
 						       .Theta_base_file_name,
+                                              IO_info.h_min_digits,
+					      hn, verbose_info
+						  .print_algorithm_details);
+			if (IO_info.output_mean_curvature)
+			   then output_gridfn(ps, gfns::gfn__mean_curvature,
+                                              "mean_curvature", cgi.GH,
+					      IO_info, IO_info
+						       .mean_curvature_base_file_name,
+                                              IO_info.h_min_digits,
 					      hn, verbose_info
 						  .print_algorithm_details);
 			}
 		}
 	}
    else {
+                struct what_to_compute new_compute_info;
 		for (int i = 0 ; i < N_horizons ; ++i)
 		{
-		expansion(NULL, 0.0,		// dummy computation
+                expansion(NULL, new_compute_info,
 			  cgi, gi,
 			  error_info, true);	// initial evaluation
 		}
@@ -402,17 +473,22 @@ const int hn = 1;
 
 struct AH_data* const AH_data_ptr = active_flag ? AH_data_array[hn]   : NULL;
 patch_system*   const      ps_ptr = active_flag ? AH_data_ptr->ps_ptr : NULL;
-const fp add_to_expansion = active_flag ? -AH_data_ptr->surface_expansion : 0.0;
+struct what_to_compute dummy_compute_info;
+struct what_to_compute & compute_info =
+	active_flag
+	? AH_data_ptr->compute_info
+	: dummy_compute_info;
 
 //
 // numerical-perturbation Jacobian
 //
 Jacobian* Jac_NP_ptr = active_flag ? AH_data_ptr->Jac_ptr : NULL;
-expansion(ps_ptr, add_to_expansion,
+expansion(ps_ptr, compute_info,
 	  cgi, gi,
 	  error_info, true);		// initial evaluation
 Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
-expansion_Jacobian(ps_ptr, Jac_NP_ptr, add_to_expansion,
+expansion_Jacobian(ps_ptr, Jac_NP_ptr,
+                   compute_info,
 		   cgi, gi, Jac_info,
 		   error_info, true,	// initial evaluation
 		   print_msg_flag);
@@ -426,12 +502,13 @@ if (test_all_Jacobian_compute_methods)
 					 *ps_ptr,
 					 verbose_info.print_algorithm_details)
 			  : NULL;
-	expansion(ps_ptr, add_to_expansion,
+	expansion(ps_ptr, compute_info,
 		  cgi, gi,
 		  error_info, true,	// initial evaluation
 		  true);		// compute SD Jacobian coeffs
 	Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
-	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr, add_to_expansion,
+	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
+                           compute_info,
 			   cgi, gi, Jac_info,
 			   error_info, true,	// initial evaluation
 			   print_msg_flag);
@@ -441,6 +518,7 @@ if (active_flag)
    then output_Jacobians(*ps_ptr,
 			 Jac_NP_ptr, Jac_SD_FDdr_ptr,
 			 IO_info, IO_info.Jacobian_base_file_name,
+			 IO_info.h_min_digits,
 			 hn, print_msg_flag);
 }
 	  }
diff --git a/src/driver/horizon_sequence.cc b/src/driver/horizon_sequence.cc
index 42345a7..e33c997 100644
--- a/src/driver/horizon_sequence.cc
+++ b/src/driver/horizon_sequence.cc
@@ -57,19 +57,19 @@ delete[] my_hn_;
 //
 char* horizon_sequence::sequence_string(const char sep[]) const
 {
-const int hn_buffer_size = 10;		// buffer size for single "%d"
-char hn_buffer[hn_buffer_size];
+const int N_hn_buffer = 10;
+char hn_buffer[N_hn_buffer];
 
-const int buffer_size = 200;		// buffer size for entire string
-static char buffer[buffer_size];
+const int N_buffer = 100;
+static char buffer[N_buffer];
 
 buffer[0] = '\0';
 	for (int pos = 0 ; pos < my_N_horizons_ ; ++pos)
 	{
 	if (pos > 0)
-	   then Util_Strlcat(buffer, sep, buffer_size);
-	snprintf(hn_buffer, hn_buffer_size, "%d", my_hn_[pos]);
-	Util_Strlcat(buffer, hn_buffer, buffer_size);
+	   then Util_Strlcat(buffer, sep, N_buffer);
+	snprintf(hn_buffer, N_hn_buffer, "%d", my_hn_[pos]);
+	Util_Strlcat(buffer, hn_buffer, N_buffer);
 	}
 
 return buffer;
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
index 5f3fae4..315a47f 100644
--- a/src/driver/initial_guess.cc
+++ b/src/driver/initial_guess.cc
@@ -17,6 +17,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -24,6 +25,7 @@
 #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"
@@ -46,7 +48,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -94,7 +95,7 @@ case initial_guess__read_from_named_file:
 
 case initial_guess__read_from_h_file:
 	input_gridfn(ps, gfns::gfn__h,
-		     IO_info, IO_info.h_base_file_name,
+		     IO_info, IO_info.h_base_file_name, IO_info.h_min_digits,
 		     hn, verbose_info.print_algorithm_highlights);
 	break;
 
diff --git a/src/driver/io.cc b/src/driver/io.cc
index e7ea3e7..6202eaa 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -9,8 +9,8 @@
 // output_gridfn - write an angular grid function to an output file
 // output_Jacobians - write a Jacobian matrix or matrices to an output file
 //
-/// io_file_name - compute file name for angular-gridfn I/O file
-/// OpenDX_control_file_name - compute file name for OpenDX control file
+/// io_ASCII_file_name - compute file name for angular-gridfn I/O file
+/// io_HDF5_file_name - compute file name for angular-gridfn I/O file
 //
 
 #include <stdio.h>
@@ -21,6 +21,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -28,6 +29,7 @@
 #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"
@@ -50,7 +52,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 //******************************************************************************
@@ -64,13 +65,15 @@ namespace {
 void output_OpenDX_control_file(const patch_system& ps,
 				const struct IO_info& IO_info,
 				int hn);
-
-const char* io_file_name(const struct IO_info& IO_info,
-			 const char base_file_name[],
-			 int hn, int AHF_iteration = 0);
-const char* OpenDX_control_file_name(const struct IO_info& IO_info,
-				     const char base_file_name[],
-				     int hn);
+const char* io_ASCII_file_name(const struct IO_info& IO_info,
+                               const char base_file_name[],
+                               int min_digits,
+                               int hn, int AHF_iteration = 0);
+const char* io_HDF5_file_name(const struct IO_info& IO_info,
+                              const char base_file_name[],
+                              int min_digits,
+                              int hn, int AHF_iteration = 0);
+void create_h_directory(const struct IO_info& IO_info);
 	  }
 
 //******************************************************************************
@@ -86,10 +89,11 @@ const char* OpenDX_control_file_name(const struct IO_info& IO_info,
 //
 void input_gridfn(patch_system& ps, int unknown_gfn,
 		  const struct IO_info& IO_info, const char base_file_name[],
+		  int min_digits,
 		  int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
 {
-const char* file_name = io_file_name(IO_info, base_file_name,
-				     hn, AHF_iteration);
+const char* file_name = io_ASCII_file_name(IO_info, base_file_name, min_digits,
+                                           hn, AHF_iteration);
 
 input_gridfn__explicit_name(ps, unknown_gfn,
 			    IO_info, base_file_name, print_msg_flag);
@@ -115,27 +119,25 @@ if (print_msg_flag)
 			   "   reading initial guess from \"%s\"", file_name);
 	}
 
-switch	(IO_info.horizon_file_format)
-	{
-case horizon_file_format__ASCII_gnuplot:
+if (IO_info.output_ASCII_files)
+   then	{
 	ps.read_ghosted_gridfn(unknown_gfn,
 			       file_name,
 			       false);		// no ghost zones in data file
-	break;
+	}
 
-case horizon_file_format__HDF5:
+else if (IO_info.output_HDF5_files)
+   then {
 	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"input_gridfn__explicit_name(): HDF5 data files not implemented yet!");
-								/*NOTREACHED*/
+"input_gridfn__explicit_name(): reading from HDF5 data files not implemented yet!");
+	}
 
-default:
+else
+	{
 	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
 "   input_gridfn__explicit_name():\n"
-"        unknown IO_info.horizon_file_format=(int)%d!\n"
-"        (this should never happen!)"
-		   ,
-		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
+"        (this should never happen!)");				/*NOTREACHED*/
 	}
 }
 
@@ -151,6 +153,19 @@ default:
 void setup_h_files(patch_system& ps, const struct IO_info& IO_info, int hn)
 {
 // create the output directory (if it doesn't already exist)
+create_h_directory(IO_info);
+output_OpenDX_control_file(ps, IO_info, hn);
+}
+
+//******************************************************************************
+
+//
+// This function creates the h ooutput directory if it does not exist
+//
+namespace {
+void create_h_directory(const struct IO_info& IO_info)
+{
+// create the output directory (if it doesn't already exist)
 const int status = CCTK_CreateDirectory(IO_info.default_directory_permission,
 					IO_info.h_directory);
 if (status < 0)
@@ -162,9 +177,8 @@ if (status < 0)
 		   ,
 		   status,
 		   IO_info.h_directory);			/*NOTREACHED*/
-
-output_OpenDX_control_file(ps, IO_info, hn);
 }
+	  }
 
 //******************************************************************************
 
@@ -178,14 +192,18 @@ void output_OpenDX_control_file(const patch_system& ps,
 				const struct IO_info& IO_info,
 				int hn)
 {
-const char* file_name = OpenDX_control_file_name(IO_info,
-						 IO_info.h_base_file_name,
-						 hn);
-FILE *fileptr = fopen(file_name, "w");
+if (! IO_info.output_OpenDX_control_files) return;
+static char file_name_buffer[IO_info::file_name_buffer_size];
+snprintf(file_name_buffer, IO_info::file_name_buffer_size,
+	 "%s/%s.ah%d.%s",
+	 IO_info.h_directory, IO_info.h_base_file_name,
+         hn, IO_info.OpenDX_control_file_name_extension);
+
+FILE *fileptr = fopen(file_name_buffer, "w");
 if (fileptr == NULL)
    then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 		"output_OpenDX_control_file(): can't open output file \"%s\"!",
-		   file_name);					/*NOTREACHED*/
+		   file_name_buffer);				/*NOTREACHED*/
 
 fprintf(fileptr, "# list the size of each patch (N_rho x N_sigma)\n");
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
@@ -222,15 +240,20 @@ fclose(fileptr);
 // FIXME: if the gridfn is not h, we assume that it's nominal-grid.
 //
 void output_gridfn(patch_system& ps, int unknown_gfn,
+                   const char gfn_name[], const cGH *cctkGH,
 		   const struct IO_info& IO_info, const char base_file_name[],
+                   int min_digits,
 		   int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
 {
-const char* file_name = io_file_name(IO_info, base_file_name,
+
+if (IO_info.output_ASCII_files)
+   then	{
+	const char* file_name
+		= io_ASCII_file_name(IO_info, base_file_name, min_digits,
 				     hn, AHF_iteration);
+	// create the output directory (if it doesn't already exist)
+	create_h_directory(IO_info);
 
-switch	(IO_info.horizon_file_format)
-	{
-case horizon_file_format__ASCII_gnuplot:
 	switch	(unknown_gfn)
 		{
 	case gfns::gfn__h:
@@ -250,6 +273,12 @@ case horizon_file_format__ASCII_gnuplot:
 				   "writing Theta to \"%s\"", file_name);
 		ps.print_gridfn(unknown_gfn, file_name);
 		break;
+	case gfns::gfn__mean_curvature:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "writing mean curvature to \"%s\"", file_name);
+		ps.print_gridfn(unknown_gfn, file_name);
+		break;
 	case gfns::gfn__Delta_h:
 		if (print_msg_flag)
 		   then CCTK_VInfo(CCTK_THORNSTRING,
@@ -264,19 +293,52 @@ case horizon_file_format__ASCII_gnuplot:
 		ps.print_gridfn(unknown_gfn, file_name);
 		break;
 		}
-	break;
+	}
 
-case horizon_file_format__HDF5:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"output_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/
+if (IO_info.output_HDF5_files)
+   then	{
+	const char* file_name
+		= io_HDF5_file_name(IO_info, base_file_name, min_digits,
+				     hn, AHF_iteration);
+	// create the output directory (if it doesn't already exist)
+	create_h_directory(IO_info);
 
-default:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   output_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n"
-"                    (this should never happen!)"
-		   ,
-		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
+	switch	(unknown_gfn)
+		{
+	case gfns::gfn__h:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "writing h to \"%s\"", file_name);
+		ps.output_ghosted_gridfn_with_xyz
+		      (unknown_gfn, gfn_name, cctkGH,
+		       true, gfns::gfn__h,
+		       file_name,
+		       IO_info.output_ghost_zones_for_h); // should we include
+							  // ghost zones?
+		break;
+	case gfns::gfn__Theta:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "writing Theta to \"%s\"", file_name);
+		ps.output_gridfn(unknown_gfn, gfn_name, cctkGH,
+                                 file_name);
+		break;
+	case gfns::gfn__Delta_h:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "writing Delta_h to \"%s\"", file_name);
+		ps.output_gridfn(unknown_gfn, gfn_name, cctkGH,
+                                 file_name);
+		break;
+	default:
+		if (print_msg_flag)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+				   "writing gfn=%d to \"%s\"",
+				   unknown_gfn, file_name);
+		ps.output_gridfn(unknown_gfn, gfn_name, cctkGH,
+                                 file_name);
+		break;
+		}
 	}
 }
 
@@ -296,17 +358,21 @@ void output_Jacobians(const patch_system& ps,
                       const Jacobian* Jac_NP_ptr,
                       const Jacobian* Jac_SD_FDdr_ptr,
 		      const struct IO_info& IO_info, const char base_file_name[],
+                      int min_digits,
 		      int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
 {
 if (Jac_NP_ptr == NULL)
    then {
-	CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "output_Jacobians(): Jac_NP_ptr == NULL is not (yet) supported!");
 	return;						// *** ERROR RETURN ***
 	}
 
-const char* file_name = io_file_name(IO_info, base_file_name,
-				     hn, AHF_iteration);
+const char* file_name = io_ASCII_file_name(IO_info, base_file_name, min_digits,
+                                           hn, AHF_iteration);
+// create the output directory (if it doesn't already exist)
+create_h_directory(IO_info);
+
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "   writing %s to \"%s\"",
@@ -427,40 +493,25 @@ fclose(fileptr);
 // result points into a private static buffer; the usual caveats apply.
 //
 namespace {
-const char* io_file_name(const struct IO_info& IO_info,
-			 const char base_file_name[],
-			 int hn, int AHF_iteration /* = 0 */)
+const char* io_ASCII_file_name(const struct IO_info& IO_info,
+                               const char base_file_name[], int min_digits,
+                               int hn, int AHF_iteration /* = 0 */)
 {
 static char file_name_buffer[IO_info::file_name_buffer_size];
 
-const char* file_name_extension;
-switch	(IO_info.horizon_file_format)
-	{
-case horizon_file_format__ASCII_gnuplot:
-	file_name_extension = IO_info.ASCII_gnuplot_file_name_extension;
-	break;
-case horizon_file_format__HDF5:
-	file_name_extension = IO_info.HDF5_file_name_extension;
-	break;
-default:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   io_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n"
-"                   (this should never happen!)"
-		   ,
-		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
-	}
+const char* file_name_extension
+	= IO_info.ASCII_gnuplot_file_name_extension;
 
 if (AHF_iteration == 0)
    then snprintf(file_name_buffer, IO_info::file_name_buffer_size,
-		 "%s/%s.t%d.ah%d.%s",
+		 "%s/%s.t%0*d.ah%d.%s",
 		 IO_info.h_directory, base_file_name,
-		 IO_info.time_iteration, hn,
+		 min_digits, IO_info.time_iteration, hn,
 		 file_name_extension);
    else snprintf(file_name_buffer, IO_info::file_name_buffer_size,
-		 "%s/%s.t%d.ah%d.it%d.%s",
+		 "%s/%s.t%0*d.ah%d.it%d.%s",
 		 IO_info.h_directory, base_file_name,
-		 IO_info.time_iteration, hn, AHF_iteration,
+		 min_digits, IO_info.time_iteration, hn, AHF_iteration,
 		 file_name_extension);
 
 return file_name_buffer;
@@ -470,48 +521,42 @@ return file_name_buffer;
 //******************************************************************************
 
 //
-// This function encapsulates our file-naming conventions for OpenDX
-// control files for our regular angular-gridfn output files.
+// This function encapsulates our file-naming conventions for angular-gridfn
+// output files (those used for h, H, and other angular grid functions).
 //
 // Arguments:
 // base_file_name[] = from the parameter file
 // hn = the horizon number
+// AHF_iteration = the apparent horizon finder's internal iteration
+//		   number (>= 1) if this is an intermediate iterate,
+//		   or the default (0) if this is a final computed
+//		   horizon position
 //
 // Results:
 // This function returns (a pointer to) the file name.  The returned
 // result points into a private static buffer; the usual caveats apply.
 //
 namespace {
-const char* OpenDX_control_file_name(const struct IO_info& IO_info,
-				     const char base_file_name[],
-				     int hn)
+const char* io_HDF5_file_name(const struct IO_info& IO_info,
+                              const char base_file_name[], int min_digits,
+                              int hn, int AHF_iteration /* = 0 */)
 {
-const char* file_name_extension;
+static char file_name_buffer[IO_info::file_name_buffer_size];
 
-switch	(IO_info.horizon_file_format)
-	{
-case horizon_file_format__ASCII_gnuplot:
-	file_name_extension = IO_info.OpenDX_control_file_name_extension;
-	break;
-case horizon_file_format__HDF5:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   OpenDX_control_file_name(): don't know how to form OpenDX control file\n"
-"                               name for HDF output files!");	/*NOTREACHED*/
-default:
-	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   OpenDX_control_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n"
-"                               (this should never happen!)"
-		   ,
-		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
-	}
+const char* file_name_extension
+	= IO_info.HDF5_file_name_extension;
 
-static char file_name_buffer[IO_info::file_name_buffer_size];
-snprintf(file_name_buffer, IO_info::file_name_buffer_size,
-	 "%s/%s.ah%d.%s",
-	 IO_info.h_directory, base_file_name, hn,
-	 file_name_extension);
+if (AHF_iteration == 0)
+   then snprintf(file_name_buffer, IO_info::file_name_buffer_size,
+		 "%s/%s.ah%d.%s",
+		 IO_info.h_directory, base_file_name,
+		 hn,
+		 file_name_extension);
+   else snprintf(file_name_buffer, IO_info::file_name_buffer_size,
+		 "%s/%s.ah%d.it%d.%s",
+		 IO_info.h_directory, base_file_name,
+		 hn, AHF_iteration,
+		 file_name_extension);
 
 return file_name_buffer;
 }
diff --git a/src/driver/make.code.defn b/src/driver/make.code.defn
index 7f9ae10..9da0915 100644
--- a/src/driver/make.code.defn
+++ b/src/driver/make.code.defn
@@ -7,21 +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 spherical_surface.cc aliased_functions.cc
+       mask.cc announce.cc aliased_functions.cc
 
 # Subdirectories containing source files
 SUBDIRS =
-
-# disable automatic template instantiation on DEC Alphas cxx compiler
-ifeq ($(shell uname), OSF1)
-  ifeq ($(CXX), cxx)
-    CXXFLAGS += -nopt
-  endif
-endif
-
-# disable automagic template instantiation on SGI Irix CC compiler
-ifneq (,$(findstring IRIX,$(shell uname)))
-  ifeq ($(notdir $(CXX)), CC)
-    CXXFLAGS += -no_auto_include
-  endif
-endif
diff --git a/src/driver/mask.cc b/src/driver/mask.cc
index 19ce183..ce7ebf1 100644
--- a/src/driver/mask.cc
+++ b/src/driver/mask.cc
@@ -21,6 +21,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "SpaceMask.h"			// from thorn SpaceMask
 
@@ -30,6 +31,7 @@
 #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"
@@ -55,7 +57,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -170,6 +171,7 @@ extern "C"
   void AHFinderDirect_maybe_do_masks(CCTK_ARGUMENTS)
 {
 DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
 
 const struct verbose_info& verbose_info = state.verbose_info;
       struct    mask_info&    mask_info = state.mask_info;
@@ -377,23 +379,20 @@ if (verbose_info.print_algorithm_highlights)
 		}
 	}
 
-if (verbose_info.print_algorithm_details)
+if (verbose_info.print_algorithm_debug)
    then {
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "      this grid has x=[%g,%g] Delta_x=%g",
+		   "      grid on this processor has x=[%g,%g]",
 		   double(mgi.global_xyz_of_ijk(X_AXIS, 0)),
-		   double(mgi.global_xyz_of_ijk(X_AXIS, mgi.proc_gridfn_dims[X_AXIS]-1)),
-		   double(mgi.coord_delta[X_AXIS]));
+		   double(mgi.global_xyz_of_ijk(X_AXIS, mgi.proc_gridfn_dims[X_AXIS]-1)));
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "                    y=[%g,%g] Delta_y=%g",
+		   "                                 y=[%g,%g]",
 		   double(mgi.global_xyz_of_ijk(Y_AXIS, 0)),
-		   double(mgi.global_xyz_of_ijk(Y_AXIS, mgi.proc_gridfn_dims[Y_AXIS]-1)),
-		   double(mgi.coord_delta[Y_AXIS]));
+		   double(mgi.global_xyz_of_ijk(Y_AXIS, mgi.proc_gridfn_dims[Y_AXIS]-1)));
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "                    z=[%g,%g] Delta_z=%g",
+		   "                                 z=[%g,%g]",
 		   double(mgi.global_xyz_of_ijk(Z_AXIS, 0)),
-		   double(mgi.global_xyz_of_ijk(Z_AXIS, mgi.proc_gridfn_dims[Z_AXIS]-1)),
-		   double(mgi.coord_delta[Z_AXIS]));
+		   double(mgi.global_xyz_of_ijk(Z_AXIS, mgi.proc_gridfn_dims[Z_AXIS]-1)));
 	}
 
 
@@ -659,7 +658,7 @@ long buffer_count = 0;
 	//
 	// FIXME:
 	// it would be more efficient here to compute the
-	// radia of a whole batch of points at once
+	// radii of a whole batch of points at once
 	const fp r_horizon = ps.radius_in_local_xyz_direction
 				   (gfns::gfn__h,
 				    local_x, local_y, local_z);
diff --git a/src/driver/misc-driver.cc b/src/driver/misc-driver.cc
index 0297529..e3cc6b0 100644
--- a/src/driver/misc-driver.cc
+++ b/src/driver/misc-driver.cc
@@ -11,6 +11,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -18,6 +19,7 @@
 #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"
@@ -40,7 +42,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index d010d32..5217097 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -5,24 +5,23 @@
 // <<<access to persistent data>>>
 //
 // AHFinderDirect_setup - top-level driver to setup persistent data structures
-// AHFinderDirect_update - top-level driver to update from steerable parameters
-//
-/// set_mask_pars - set internal data structures from mask parameters
 ///
 /// decode_method - decode the  method  parameter
 /// decode_verbose_level - decode the  verbose_level  parameter
-/// decode_horizon_file_format - decode the  horizon_file_format  parameter
 ///
 /// allocate_horizons_to_processor - choose which horizons this proc will find
 ///
 /// choose_patch_system_type - choose patch system type
 ///
 
+#include <assert.h>
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
 #include <string.h>
 
+#include <vector>
+
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -36,6 +35,7 @@
 #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"
@@ -58,7 +58,6 @@
 // all the code in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************
 
@@ -67,17 +66,14 @@ using jtutil::error_exit;
 //
 
 namespace {
-void set_mask_pars(CCTK_ARGUMENTS);
-
 enum method
   decode_method(const char method_string[]);
 enum verbose_level
   decode_verbose_level(const char verbose_level_string[]);
-enum horizon_file_format
-  decode_horizon_file_format(const char horizon_file_format_string[]);
 
 int allocate_horizons_to_processor(int N_procs, int my_proc,
 				   int N_horizons, bool multiproc_flag,
+                                   const CCTK_INT depends_on[],
 				   horizon_sequence& my_hs,
 				   const struct verbose_info& verbose_info);
 
@@ -113,9 +109,36 @@ CCTK_VInfo(CCTK_THORNSTRING,
 
 
 //
+// check parameters
+//
+int need_zones = 0;
+for (int n=1; n<N_horizons; ++n) {
+  need_zones = jtutil::max(need_zones, N_zones_per_right_angle[n]);
+}
+if (need_zones > max_N_zones_per_right_angle) {
+  CCTK_VWarn (FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+              "AHFinderDirect_setup(): "
+              "The parameter max_N_zones_per_right_angle must be at least the maximum of all N_zones_per_right_angle[] parameters.  "
+              "Set max_N_zones_per_right_angle to %d or higher to continue.",
+              need_zones);      /*NOTREACHED*/
+ }
+for (int n=1; n<N_horizons; ++n) {
+  if (depends_on[n] != 0) {
+    assert (depends_on[n] >= 1 && depends_on[n] < n);
+    if (N_zones_per_right_angle[n] != N_zones_per_right_angle[depends_on[n]]) {
+      CCTK_VWarn (FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "AHFinderDirect_setup(): "
+                  "The parameter N_zones_per_right_angle must be the same for a horizon as for the horizon on which it depends.  "
+                  "Horizon %d depends on horizon %d, but they have different resolutions.",
+                  n, depends_on[n]); /*NOTREACHED*/
+    }
+  }
+}
+
+
+//
 // basic setup
 //
-state.find_every = find_every;
 state.method = decode_method(method);
 
 state.error_info.warn_level__point_outside__initial
@@ -144,9 +167,7 @@ verbose_info.print_algorithm_details
 verbose_info.print_algorithm_debug
    = (state.verbose_info.verbose_level >= verbose_level__algorithm_debug);
 
-state.metric_type = /* ADMBase:: */ metric_type;
-
-state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1;
+state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreate("finding apparent horizons") : -1;
 
 state.N_procs = CCTK_nProcs(cctkGH);
 state.my_proc = CCTK_MyProc(cctkGH);
@@ -173,6 +194,7 @@ if (cgi.coord_system_handle < 0)
 		   coordinate_system_name);			/*NOTREACHED*/
 cgi.use_Cactus_conformal_metric = false;	// dummy value, may change later
 
+cgi.mask_varindex    = Cactus_gridfn_varindex("AHFinderDirect::ahmask");
 cgi.g_dd_11_varindex = Cactus_gridfn_varindex("ADMBase::gxx");
 cgi.g_dd_12_varindex = Cactus_gridfn_varindex("ADMBase::gxy");
 cgi.g_dd_13_varindex = Cactus_gridfn_varindex("ADMBase::gxz");
@@ -222,9 +244,6 @@ gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0);
 //
 // Jacobian info
 //
-state.test_all_Jacobian_compute_methods
-	= (test_all_Jacobian_compute_methods != 0);
-
 struct Jacobian_info& Jac_info = state.Jac_info;
 Jac_info.Jacobian_compute_method
 	= decode_Jacobian_compute_method(Jacobian_compute_method);
@@ -243,6 +262,8 @@ solver_info.linear_solver_pars.ILUCG_pars.error_tolerance
 	= ILUCG__error_tolerance;
 solver_info.linear_solver_pars.ILUCG_pars.limit_CG_iterations
 	= (ILUCG__limit_CG_iterations != 0);
+solver_info.linear_solver_pars.UMFPACK_pars.N_II_iterations
+	= UMFPACK__N_II_iterations;
 solver_info.max_Newton_iterations__initial
 	= max_Newton_iterations__initial;
 solver_info.max_Newton_iterations__subsequent
@@ -250,9 +271,13 @@ solver_info.max_Newton_iterations__subsequent
 solver_info.max_allowable_Delta_h_over_h = max_allowable_Delta_h_over_h;
 solver_info.Theta_norm_for_convergence   = Theta_norm_for_convergence;
 solver_info.max_allowable_Theta          = max_allowable_Theta;
+solver_info.max_allowable_Theta_growth_iterations
+                                       = max_allowable_Theta_growth_iterations;
+solver_info.max_allowable_Theta_nonshrink_iterations
+                                    = max_allowable_Theta_nonshrink_iterations;
 // ... horizon numbers run from 1 to N_horizons inclusive
 //     so the array size is N_horizons+1
-solver_info.max_allowable_horizon_radius = new double[state.N_horizons+1];
+solver_info.max_allowable_horizon_radius = new fp[state.N_horizons+1];
 	  {
 	for (int hn = 0 ; hn <= N_horizons ; ++hn)
 	{
@@ -260,18 +285,22 @@ solver_info.max_allowable_horizon_radius = new double[state.N_horizons+1];
 		= max_allowable_horizon_radius[hn];
 	}
 	  }
+solver_info.want_expansion_gradients = want_expansion_gradients;
 
 
 //
 // I/O info
 //
 struct IO_info& IO_info = state.IO_info;
-IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
+IO_info.output_ASCII_files = (output_ASCII_files != 0);
+IO_info.output_HDF5_files = (output_HDF5_files != 0);
 IO_info.output_initial_guess = (output_initial_guess != 0);
 IO_info.output_h_every     = output_h_every;
 IO_info.output_Theta_every = output_Theta_every;
+IO_info.output_mean_curvature_every = output_mean_curvature_every;
 IO_info.output_h     = false;	// dummy value
 IO_info.output_Theta = false;	// dummy value
+IO_info.output_mean_curvature = false;	// dummy value
 
 IO_info.output_BH_diagnostics              = (output_BH_diagnostics != 0);
 IO_info.BH_diagnostics_directory
@@ -290,7 +319,9 @@ IO_info.h_directory
 	  : h_directory;
 IO_info.h_base_file_name         = h_base_file_name;
 IO_info.Theta_base_file_name     = Theta_base_file_name;
+IO_info.mean_curvature_base_file_name     = mean_curvature_base_file_name;
 IO_info.Delta_h_base_file_name   = Delta_h_base_file_name;
+IO_info.h_min_digits             = h_min_digits;
 IO_info.Jacobian_base_file_name  = Jacobian_base_file_name;
 IO_info.output_OpenDX_control_files  = (output_OpenDX_control_files != 0);
 IO_info.OpenDX_control_file_name_extension = OpenDX_control_file_name_extension;
@@ -304,8 +335,6 @@ IO_info.time           = 0.0;
 state.BH_diagnostics_info.integral_method
    = patch::decode_integration_method(integral_method);
 
-state.always_broadcast_horizon_shape = (always_broadcast_horizon_shape != 0);
-
 
 //
 // mask parameters
@@ -326,33 +355,43 @@ mask_info.set_mask_for_this_horizon = new bool[N_horizons+1];
 	}
 	  }
 if (mask_info.set_mask_for_any_horizon)
-   then set_mask_pars(CCTK_PASS_CTOC);
-
-
-//
-// announce parameters
-//
-state.announce_centroid_flag = (which_horizon_to_announce_centroid != 0);
-if (state.announce_centroid_flag)
    then {
-	state.which_horizon_to_announce_centroid
-		= which_horizon_to_announce_centroid;
-	if (    (state.which_horizon_to_announce_centroid < 1)
-	     || (state.which_horizon_to_announce_centroid > state.N_horizons) )
-	   then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   AHFinderDirect_setup():\n"
-"        invalid which_horizon_to_announce_centroid = %d!\n"
-"        (valid range is [1,N_horizons=%d])!\n"
-		   ,
-		   state.which_horizon_to_announce_centroid,
-		   state.N_horizons);				/*NOTREACHED*/
+	mask_info.radius_multiplier  = mask_radius_multiplier;
+	mask_info.radius_offset      = mask_radius_offset;
+	mask_info.buffer_thickness   = mask_buffer_thickness;
+	mask_info.mask_is_noshrink   = mask_is_noshrink;
+	mask_info.min_horizon_radius_points_for_mask
+				     = min_horizon_radius_points_for_mask;
+	mask_info.set_old_style_mask = (set_old_style_mask != 0);
+	mask_info.set_new_style_mask = (set_new_style_mask != 0);
+	if (mask_info.set_old_style_mask)
+	   then {
+		struct mask_info::old_style_mask_info& osmi
+			= mask_info.old_style_mask_info;
+		osmi.gridfn_name     = old_style_mask_gridfn_name;
+		osmi.gridfn_varindex = Cactus_gridfn_varindex(osmi.gridfn_name);
+		osmi.gridfn_dataptr  = NULL;	// dummy value; fixup later
+		osmi.inside_value  = old_style_mask_inside_value;
+		osmi.buffer_value  = old_style_mask_buffer_value;
+		osmi.outside_value = old_style_mask_outside_value;
+		}
+	if (mask_info.set_new_style_mask)
+	   then {
+		struct mask_info::new_style_mask_info& nsmi
+			= mask_info.new_style_mask_info;
+		nsmi.gridfn_name     = new_style_mask_gridfn_name;
+		nsmi.gridfn_varindex = Cactus_gridfn_varindex(nsmi.gridfn_name);
+		nsmi.gridfn_dataptr  = NULL;	// dummy value; fixup later
+		nsmi.bitfield_name   = new_style_mask_bitfield_name;
+		nsmi.bitfield_bitmask = 0;	// dummy value; fixup later
+		nsmi.inside_value    = new_style_mask_inside_value;
+		nsmi.buffer_value    = new_style_mask_buffer_value;
+		nsmi.outside_value   = new_style_mask_outside_value;
+		nsmi.inside_bitvalue = 0;	// dummy value; fixup later
+		nsmi.buffer_bitvalue = 0;	// dummy value; fixup later
+		nsmi.outside_bitvalue = 0;	// dummy value; fixup later
+		}
 	}
-   else state.which_horizon_to_announce_centroid = 0;	// dummy value; unused
-
-// initial value, will get each horizon's  store_info_in_SS_vars
-// flag inclusive-ored into it below
-state.store_info_in_SS_vars_for_any_horizon = false;
 
 
 //
@@ -371,6 +410,7 @@ const bool multiproc_flag = (state.method == method__find_horizons);
 state.N_active_procs
    = allocate_horizons_to_processor(state.N_procs, state.my_proc,
 				    state.N_horizons, multiproc_flag,
+                                    depends_on,
 				    hs,
 				    verbose_info);
 
@@ -474,6 +514,8 @@ if (strlen(surface_interpolator_name) > 0)
 			      true, verbose_info.print_algorithm_details);
 	patch_system& ps = *AH_data.ps_ptr;
 	if (genuine_flag)
+	   then ps.set_gridfn_to_constant(0.0, gfns::gfn__zero);
+	if (genuine_flag)
 	   then ps.set_gridfn_to_constant(1.0, gfns::gfn__one);
 
 	AH_data.Jac_ptr = genuine_flag
@@ -482,9 +524,67 @@ if (strlen(surface_interpolator_name) > 0)
 					 verbose_info.print_algorithm_details)
 			  : NULL;
 
-	AH_data.surface_expansion = surface_expansion[hn];
-
-	AH_data.initial_find_flag = genuine_flag;
+        AH_data.compute_info.surface_definition =
+          STRING_EQUAL(surface_definition[hn], "expansion")
+          ? definition_expansion
+          : STRING_EQUAL(surface_definition[hn], "inner expansion")
+          ? definition_inner_expansion
+          : STRING_EQUAL(surface_definition[hn], "mean curvature")
+          ? definition_mean_curvature
+          : STRING_EQUAL(surface_definition[hn], "expansion product")
+          ? definition_expansion_product
+          : (CCTK_WARN (0, "internal error"), definition_error);
+        AH_data.compute_info.surface_modification =
+          STRING_EQUAL(surface_modification[hn], "none")
+          ? modification_none
+          : STRING_EQUAL(surface_modification[hn], "radius")
+          ? modification_radius
+          : STRING_EQUAL(surface_modification[hn], "radius^2")
+          ? modification_radius2
+#if 0
+          : STRING_EQUAL(surface_modification[hn], "mean radius")
+          ? modification_mean_radius
+          : STRING_EQUAL(surface_modification[hn], "areal radius")
+          ? modification_areal_radius
+#endif
+          : (CCTK_WARN (0, "internal error"), modification_error);
+        AH_data.compute_info.surface_selection = 
+          STRING_EQUAL(surface_selection[hn], "definition")
+          ? selection_definition
+          : STRING_EQUAL(surface_selection[hn], "mean coordinate radius")
+          ? selection_mean_coordinate_radius
+          : STRING_EQUAL(surface_selection[hn], "areal radius")
+          ? selection_areal_radius
+          : STRING_EQUAL(surface_selection[hn], "expansion times mean coordinate radius")
+          ? selection_expansion_mean_coordinate_radius
+          : STRING_EQUAL(surface_selection[hn], "expansion times areal radius")
+          ? selection_expansion_areal_radius
+          : (CCTK_WARN (0, "internal error"), selection_error);
+	AH_data.compute_info.desired_value = desired_value[hn];
+
+        AH_data.move_origins               = move_origins;
+
+        AH_data.use_pretracking            = use_pretracking[hn];
+        AH_data.pretracking_max_iterations = pretracking_max_iterations[hn];
+
+        AH_data.pretracking_value          = pretracking_value[hn];
+        AH_data.pretracking_minimum_value  = pretracking_minimum_value[hn];
+        AH_data.pretracking_maximum_value  = pretracking_maximum_value[hn];
+        AH_data.pretracking_delta          = pretracking_delta[hn];
+        AH_data.pretracking_minimum_delta  = pretracking_minimum_delta[hn];
+        AH_data.pretracking_maximum_delta  = pretracking_maximum_delta[hn];
+
+        AH_data.depends_on           = depends_on[hn];
+        AH_data.desired_value_factor = desired_value_factor[hn];
+        AH_data.desired_value_offset = desired_value_offset[hn];
+
+        AH_data.shiftout_factor = shiftout_factor[hn];
+        AH_data.smoothing_factor = smoothing_factor[hn];
+
+	// AH_data.initial_find_flag = genuine_flag;
+	// AH_data.really_initial_find_flag = AH_data.initial_find_flag;
+	AH_data.initial_find_flag = true;
+	AH_data.really_initial_find_flag = AH_data.initial_find_flag;
 
 	if (genuine_flag)
 	   then {
@@ -493,6 +593,8 @@ if (strlen(surface_interpolator_name) > 0)
 			   "      setting initial guess parameters etc");
 		AH_data.initial_guess_info.method
 			= decode_initial_guess_method(initial_guess_method[hn]);
+                AH_data.initial_guess_info.reset_horizon_after_not_finding
+                        = reset_horizon_after_not_finding[hn];
 		// ... read from named file
 		AH_data.initial_guess_info.read_from_named_file_info.file_name
 			= initial_guess__read_from_named_file__file_name[hn];
@@ -542,106 +644,47 @@ if (strlen(surface_interpolator_name) > 0)
 			= initial_guess__coord_ellipsoid__z_radius[hn];
 		}
 
+	AH_data.search_flag = false;
 	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
-
-	state.store_info_in_SS_vars_for_any_horizon
-		|= AH_data.store_info_in_SS_vars;
 	}
 	  }
-}
 
-//******************************************************************************
-
-//
-// This function is called by the Cactus scheduler to update some of our
-// persistent data structures (stored in  struct state ) from parameters
-// which may have been steered.  It is also called by  AHFinderDirect_setup()
-// (above)  as part of the initial setup of these data structures.
-//
-extern "C"
-  void AHFinderDirect_update(CCTK_ARGUMENTS)
-{
-set_mask_pars(CCTK_PASS_CTOC);
+  // Initialise the centroids.  These values may later be overwritten
+  // when they are recovered from a checkpoint.  However, if new
+  // horizons are enabled during recovery, they are then correctly
+  // initialised.
+  for (int n = 0; n < N_horizons; ++ n) {
+    // Horizon centroids are not valid
+    ah_centroid_x[n] = 0.0;
+    ah_centroid_y[n] = 0.0;
+    ah_centroid_z[n] = 0.0;
+    ah_centroid_t[n] = 0.0;
+    ah_centroid_valid[n] = 0;
+    ah_centroid_iteration[n] = -1;
+    
+    ah_centroid_x_p[n] = 0.0;
+    ah_centroid_y_p[n] = 0.0;
+    ah_centroid_z_p[n] = 0.0;
+    ah_centroid_t_p[n] = 0.0;
+    ah_centroid_valid_p[n] = 0;
+    ah_centroid_iteration_p[n] = -1;
+    
+    struct AH_data& AH_data = *state.AH_data_array[n+1];
+    ah_initial_find_flag[n]        = AH_data.initial_find_flag;
+    ah_really_initial_find_flag[n] = AH_data.really_initial_find_flag;
+    ah_search_flag[n]              = AH_data.search_flag;
+    ah_found_flag[n]               = AH_data.found_flag;
+    if (verbose_info.print_algorithm_highlights) {
+      printf ("AHF setup %d initial_find_flag=%d\n",        n+1, (int) AH_data.initial_find_flag);
+      printf ("AHF setup %d really_initial_find_flag=%d\n", n+1, (int) AH_data.really_initial_find_flag);
+      printf ("AHF setup %d search_flag=%d\n",              n+1, (int) AH_data.search_flag);
+      printf ("AHF setup %d found_flag=%d\n",               n+1, (int) AH_data.found_flag);
+    }
+  }
 }
 
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function sets our internal data structures from the the mask
-// parameters.
-//
-namespace {
-void set_mask_pars(CCTK_ARGUMENTS)
-{
-DECLARE_CCTK_ARGUMENTS
-DECLARE_CCTK_PARAMETERS
-
-struct mask_info& mask_info = state.mask_info;
-mask_info.set_old_style_mask = (set_old_style_mask != 0);
-mask_info.set_new_style_mask = (set_new_style_mask != 0);
-if (mask_info.set_mask_for_any_horizon)
-   then {
-	mask_info.radius_multiplier  = mask_radius_multiplier;
-	mask_info.radius_offset      = mask_radius_offset;
-	mask_info.buffer_thickness   = mask_buffer_thickness;
-	mask_info.mask_is_noshrink   = mask_is_noshrink;
-	mask_info.min_horizon_radius_points_for_mask
-				     = min_horizon_radius_points_for_mask;
-	mask_info.set_old_style_mask = (set_old_style_mask != 0);
-	mask_info.set_new_style_mask = (set_new_style_mask != 0);
-	if (mask_info.set_old_style_mask)
-	   then {
-		struct mask_info::old_style_mask_info& osmi
-			= mask_info.old_style_mask_info;
-		osmi.gridfn_name     = old_style_mask_gridfn_name;
-		osmi.gridfn_varindex = Cactus_gridfn_varindex(osmi.gridfn_name);
-		osmi.gridfn_dataptr  = NULL;	// dummy value; fixup later
-		osmi.inside_value  = old_style_mask_inside_value;
-		osmi.buffer_value  = old_style_mask_buffer_value;
-		osmi.outside_value = old_style_mask_outside_value;
-		}
-	if (mask_info.set_new_style_mask)
-	   then {
-		struct mask_info::new_style_mask_info& nsmi
-			= mask_info.new_style_mask_info;
-		nsmi.gridfn_name     = new_style_mask_gridfn_name;
-		nsmi.gridfn_varindex = Cactus_gridfn_varindex(nsmi.gridfn_name);
-		nsmi.gridfn_dataptr  = NULL;	// dummy value; fixup later
-		nsmi.bitfield_name   = new_style_mask_bitfield_name;
-		nsmi.bitfield_bitmask = 0;	// dummy value; fixup later
-		nsmi.inside_value    = new_style_mask_inside_value;
-		nsmi.buffer_value    = new_style_mask_buffer_value;
-		nsmi.outside_value   = new_style_mask_outside_value;
-		nsmi.inside_bitvalue = 0;	// dummy value; fixup later
-		nsmi.buffer_bitvalue = 0;	// dummy value; fixup later
-		nsmi.outside_bitvalue = 0;	// dummy value; fixup later
-		}
-	}
-}
-	  }
 
 //******************************************************************************
 //******************************************************************************
@@ -677,9 +720,7 @@ namespace {
 enum verbose_level
   decode_verbose_level(const char verbose_level_string[])
 {
-if	(STRING_EQUAL(verbose_level_string, "no output"))
-   then return verbose_level__no_output;
-else if (STRING_EQUAL(verbose_level_string, "physics highlights"))
+if	(STRING_EQUAL(verbose_level_string, "physics highlights"))
    then return verbose_level__physics_highlights;
 else if (STRING_EQUAL(verbose_level_string, "physics details"))
    then return verbose_level__physics_details;
@@ -696,26 +737,6 @@ else	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 	  }
 
 //******************************************************************************
-
-//
-// This function decodes the  horizon_file_format  parameter (string) into an
-// internal enum for future use.
-//
-namespace {
-enum horizon_file_format
-  decode_horizon_file_format(const char horizon_file_format_string[])
-{
-if	(STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)"))
-   then return horizon_file_format__ASCII_gnuplot;
-else if (STRING_EQUAL(horizon_file_format_string, "HDF5"))
-   then return horizon_file_format__HDF5;
-else	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
-"decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!",
-		   horizon_file_format_string);			/*NOTREACHED*/
-}
-	  }
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
@@ -743,6 +764,7 @@ else	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
 namespace {
 int allocate_horizons_to_processor(int N_procs, int my_proc,
 				   int N_horizons, bool multiproc_flag,
+                                   const CCTK_INT depends_on[],
 				   horizon_sequence& my_hs,
 				   const struct verbose_info& verbose_info)
 {
@@ -753,15 +775,38 @@ const int N_active_procs = multiproc_flag ? jtutil::min(N_procs, N_horizons)
 // Implementation note:
 // We allocate the horizons to active processors in round-robin order.
 //
+std::vector<int> proc_of_horizon (N_horizons+1);
+	for (int hn = 1 ; hn <= N_horizons ; ++hn)
+        {
+        proc_of_horizon.at(hn) = -1;
+        }
+
 int proc = 0;
 	for (int hn = 1 ; hn <= N_horizons ; ++hn)
 	{
+        if (depends_on[hn] < 0 || depends_on[hn] > N_horizons) {
+          CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "horizon %d depends on a horizon with the illegal index %d",
+                     hn, int(depends_on[hn]));
+        } else if (depends_on[hn] == hn) {
+          CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "horizon %d depends on itself",
+                     hn);
+        } else if (depends_on[hn] > hn) {
+          CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "horizon %d depends on a horizon with a larger index %d",
+                     hn, int(depends_on[hn]));
+        }
+        const int this_horizons_proc
+          = depends_on[hn] == 0 ? proc : proc_of_horizon.at(depends_on[hn]);
+        assert (this_horizons_proc >= 0 && this_horizons_proc < N_procs);
+        proc_of_horizon.at(hn) = this_horizons_proc;
 	if (verbose_info.print_algorithm_highlights)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "   allocating horizon %d to processor #%d",
-			   hn, proc);
-	if (proc == my_proc)
-	   then my_hs.append_hn(hn);
+			   hn, this_horizons_proc);
+	if (this_horizons_proc == my_proc)
+           then my_hs.append_hn(hn);
 	if (++proc >= N_active_procs)
 	   then proc = 0;
 	}
@@ -875,6 +920,19 @@ else if (STRING_EQUAL(grid_domain, "octant"))
 		   "      ==> using \"+xz quadrant (mirrored)\" patch system");
 		return patch_system::patch_system__plus_xz_quadrant_mirrored;
 		}
+	else if ((origin_y == 0.0) && (origin_z == 0.0))
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      Cactus has mirrored octant mode");
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      but patch system origin_x=%g != 0",
+		   double(origin_z));
+		CCTK_VInfo(CCTK_THORNSTRING,
+		   "      ==> using \"+yz quadrant (mirrored)\" patch system");
+		//return patch_system::patch_system__plus_yz_quadrant_mirrored;
+                CCTK_WARN (0, "This patch system type is not implemented");
+		return patch_system::patch_system__full_sphere;
+		}
 	else	{
 		CCTK_VInfo(CCTK_THORNSTRING,
 		   "      Cactus has mirrored octant mode");
@@ -888,7 +946,7 @@ else if (STRING_EQUAL(grid_domain, "octant"))
 	}
 
 else	{
-	CCTK_VWarn(SERIOUS_WARNING, __LINE__, __FILE__, CCTK_THORNSTRING,
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
 "   choose_patch_system_type()\n"
 "        grid::domain = \"%s\" not supported!\n"
diff --git a/src/driver/spherical_surface.cc b/src/driver/spherical_surface.cc
deleted file mode 100644
index 96d3bc8..0000000
--- a/src/driver/spherical_surface.cc
+++ /dev/null
@@ -1,264 +0,0 @@
-// 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"
-
-#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
-	  {
-using jtutil::error_exit;
-
-//******************************************************************************
-
-//
-// ***** 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.mean_radius;
-/* SphericalSurface:: */ sf_min_radius   [sn] = BH_diagnostics.min_radius;
-/* SphericalSurface:: */ sf_max_radius   [sn] = BH_diagnostics.max_radius;
-
-/* SphericalSurface:: */ sf_area         [sn] = BH_diagnostics.area;
-
-/* 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.
-// * It would be cleaner to abstract out the (complicated) indexing
-//   calculation for the SphericalSurface shape array into a separate
-//   access structure/function, the way we do with  struct cactus_grid_info
-//   for grid functions.
-//
-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 local_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] = local_r;
-		}
-	}
-}
-	  }
-
-//******************************************************************************
-
-	  }	// namespace AHFinderDirect
diff --git a/src/driver/state.cc b/src/driver/state.cc
index dde1e23..9d495be 100644
--- a/src/driver/state.cc
+++ b/src/driver/state.cc
@@ -8,6 +8,7 @@
 #include "util_Table.h"
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 #include "config.h"
 #include "stdc.h"
@@ -15,6 +16,7 @@
 #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"
@@ -37,7 +39,6 @@
 // everything in this file is inside this namespace
 namespace AHFinderDirect
 	  {
-using jtutil::error_exit;
 
 //******************************************************************************