move high-level driver routines from gr/ to new driver/ directory

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

12 files changed, 57 insertions, 1536 deletions

diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
deleted file mode 100644
index 009f5f0..0000000
--- a/src/gr/AHFinderDirect.hh
+++ /dev/null
@@ -1,156 +0,0 @@
-// AHFinderDirect.hh -- misc global-within-this-thorn stuff
-// $Header$
-
-//******************************************************************************
-
-//
-// number of spatial dimensions in the main Cactus grid
-// and in our trial-horizon-surface grid
-//
-enum	{ N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 };
-
-//
-// this enum holds the (a) decoded  Jacobian_method  parameter,
-// i.e. it specifies how we compute the (a) Jacobian matrix
-//
-enum	Jacobian_method
-	{
-	Jacobian__numerical_perturb,
-	Jacobian__symbolic_diff_with_FD_dr,
-	Jacobian__symbolic_diff // no comma
-	};
-
-//
-// this enum holds the (a) decoded  geometry_method  parameter,
-// i.e. it specifies how we compute the (a) geometry
-//
-enum	geometry_method
-	{
-	geometry__interpolate_from_Cactus_grid,
-	geometry__Schwarzschild_EF // no comma
-	};
-
-//
-// This structure holds information for computing the spacetime geometry.
-// This is normally done by interpolating $g_{ij}$ and $K_{ij}$ from the
-// usual Cactus grid, but can optionally instead by done by hard-wiring
-// the Schwarzschild geometry in Eddington-Finkelstein coordinates.
-//
-struct	geometry_info
-	{
-	enum geometry_method geometry_method;
-
-	// parameters for geometry_method = interpolate_from_Cactus_grid
-	int operator_handle;		// Cactus handle to interpolation op
-	int param_table_handle;		// Cactus handle to parameter table
-					// for the interpolator
-
-	// parameters for geometry_method = Schwarzschild_EF
-	fp geometry__Schwarzschild_EF__x_posn;	// x posn of Schwarzschild BH
-	fp geometry__Schwarzschild_EF__y_posn;	// y posn of Schwarzschild BH
-	fp geometry__Schwarzschild_EF__z_posn;	// z posn of Schwarzschild BH
-	fp geometry__Schwarzschild_EF__mass;	// mass of Schwarzschild BH
-	fp geometry__Schwarzschild_EF__epsilon;	// use z axis limits if
-						// (x^2+y^2)/r^2 <= this
-	fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD
-						// computation of partial_k g_ij
-	};
-
-//
-// This structure holds all the information we need about the Cactus grid
-// and gridfns outside the top-level driver.  At present we use the
-//  CCTK_InterpLocalUniform()  local interpolator to access the Cactus
-// gridfns.  Much of the information in this structure is specific to
-// that API, and (alas) will need changing when we eventually switch to
-// a "global" multiprocessor/distributed interpolator.
-//
-struct	cactus_grid_info
-	{
-	cGH *GH;			// --> Cactus grid hierarchy
-
-	// Cactus coordinate system
-	fp coord_origin[N_GRID_DIMS];	// (x,y,z) of grid posn (0,0,0)
-	fp coord_delta[N_GRID_DIMS];	// (x,y,z) grid spacing
-
-	// dimensions of gridfn data, viewed as a 3-D array
-	// n.b. storage ordering is Fortran,
-	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
-	CCTK_INT gridfn_dims[N_GRID_DIMS];
-
-	// --> Cactus gridfn data for grid posn (0,0,0)
-	const fp* g_dd_11_data;
-	const fp* g_dd_12_data;
-	const fp* g_dd_13_data;
-	const fp* g_dd_22_data;
-	const fp* g_dd_23_data;
-	const fp* g_dd_33_data;
-	const fp* K_dd_11_data;
-	const fp* K_dd_12_data;
-	const fp* K_dd_13_data;
-	const fp* K_dd_22_data;
-	const fp* K_dd_23_data;
-	const fp* K_dd_33_data;
-	};
-
-//
-// This struct holds parameters for computing the Jacobian matrix.
-//
-struct	Jacobian_info
-	{
-	enum Jacobian_method Jacobian_method;
-	fp perturbation_amplitude;
-	};
-
-//
-// This struct holds parameters for solving the H(h) = 0 equations.
-//
-struct	solver_info
-	{
-	// stuff for Newton_solve()
-	int max_Newton_iterations;
-	bool output_h_and_H_at_each_Newton_iteration;
-	const char *h_file_name;
-	const char *H_of_h_file_name;
-	fp H_norm_for_convergence;
-	fp Delta_h_norm_for_convergence;
-	bool final_H_update_if_exit_x_H_small;
-	};
-
-//******************************************************************************
-
-//
-// prototypes for functions visible outside their source files
-//
-
-// driver.cc
-extern "C"
-  void AHFinderDirect_driver(CCTK_ARGUMENTS);
-
-// horizon_function.cc
-void horizon_function(patch_system& ps,
-		      const struct cactus_grid_info& cgi,
-		      const struct geometry_info& gi,
-		      bool Jacobian_flag = false,
-		      bool msg_flag = false,
-		      jtutil::norm<fp>* H_norms_ptr = NULL);
-
-// horizon_Jacobian.cc
-void horizon_Jacobian(patch_system& ps,
-		      const struct cactus_grid_info& cgi,
-		      const struct geometry_info& gi,
-		      const struct Jacobian_info& Jacobian_info,
-		      Jacobian& Jac);
-
-// Schwarzschild_EF.cc
-void Schwarzschild_EF_geometry(patch_system& ps,
-			       const struct geometry_info& gi,
-			       bool msg_flag);
-
-// Newton.cc
-// return true for success, false for failure to converge
-bool Newton_solve(patch_system& ps,
-		  const struct cactus_grid_info& cgi,
-		  const struct geometry_info& gi,
-		  const struct Jacobian_info& Jacobian_info,
-		  const struct solver_info& solver_info,
-		  Jacobian& Jac);
diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc
deleted file mode 100644
index 79f1bab..0000000
--- a/src/gr/Newton.cc
+++ /dev/null
@@ -1,180 +0,0 @@
-// Newton.cc -- solve H(h) = 0 via Newton's method
-// $Id$
-//
-// <<<prototypes for functions local to this file>>>
-// Newton_solve - driver to solve H(h) = 0 via Newton's method
-//
-
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-#include <vector>
-
-#include "util_Table.h"
-#include "cctk.h"
-#include "cctk_Arguments.h"
-
-#include "stdc.h"
-#include "config.hh"
-#include "../jtutil/util.hh"
-#include "../jtutil/array.hh"
-#include "../jtutil/cpm_map.hh"
-#include "../jtutil/linear_map.hh"
-using jtutil::error_exit;
-
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
-
-#include "../elliptic/Jacobian.hh"
-
-#include "gfns.hh"
-#include "AHFinderDirect.hh"
-
-//******************************************************************************
-
-//
-// ***** prototypes for functions local to this file *****
-//
-
-namespace {
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function solves H(h) = 0 for h via Newton's method.
-//
-bool Newton_solve(patch_system& ps,
-		  const struct cactus_grid_info& cgi,
-		  const struct geometry_info& gi,
-		  const struct Jacobian_info& Jacobian_info,
-		  const struct solver_info& solver_info,
-		  Jacobian& Jac)
-{
-	for (int iteration = 1 ;
-	     iteration <= solver_info.max_Newton_iterations ;
-	     ++iteration)
-	{
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "Newton iteration %d", iteration);
-
-	jtutil::norm<fp> H_norms;
-	horizon_function(ps, cgi, gi, true, true, &H_norms);
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   H rms-norm=%.2e, infinity-norm=%.2e",
-		   H_norms.rms_norm(), H_norms.infinity_norm());
-
-	if (solver_info.output_h_and_H_at_each_Newton_iteration)
-	   then {
-		const char N_buffer = 100;
-		char file_name_buffer[100];
-
-		snprintf(file_name_buffer, N_buffer,
-			 "%s.it%d",
-			 solver_info.h_file_name, iteration);
-		CCTK_VInfo(CCTK_THORNSTRING,
-			   "   writing h to \"%s\"",
-			   file_name_buffer);
-		ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
-						 true, gfns::gfn__h,
-						 file_name_buffer,
-						 false);   // no ghost zones
-        
-		snprintf(file_name_buffer, N_buffer,
-			 "%s.it%d",
-			 solver_info.H_of_h_file_name, iteration);
-		CCTK_VInfo(CCTK_THORNSTRING,
-			   "   writing H(h) to \"%s\"",
-			   file_name_buffer);
-		ps.print_gridfn_with_xyz(gfns::gfn__H,
-					 true, gfns::gfn__h,
-					 file_name_buffer);
-		}
-
-	if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
-	   then {
-		CCTK_VInfo(CCTK_THORNSTRING,
-			   "==> finished (||H|| < %g)",
-			   double(solver_info.H_norm_for_convergence));
-		return true;				// *** NORMAL RETURN ***
-		}
-
-	// compute the Newton step
-	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac);
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "solving linear system for Delta_h (%d unknowns)",
-		   Jac.NN());
-	const fp rcond = Jac.solve_linear_system(gfns::gfn__H,
-						 gfns::gfn__Delta_h);
-	if	(rcond == 0.0)
-	   then {
-		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "Newton_solve: Jacobian matrix is numerically singular!");
-		return false;				// *** ERROR RETURN ***
-		}
-	if (rcond > 0.0)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "done solving linear system (condition number %.2e)",
-			   double(rcond));
-	   else CCTK_VInfo(CCTK_THORNSTRING,
-			   "done solving linear system");
-
-	// take the Newton step
-	jtutil::norm<fp> Delta_h_norms;
-		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 Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma);
-		Delta_h_norms.data(Delta_h);
-
-		p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h;
-		}
-		}
-		}
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "h += Delta_h (rms-norm=%.2e, infinity-norm=%.2e)",
-		   Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm());
-
-	if (Delta_h_norms.infinity_norm() <= solver_info
-					     .Delta_h_norm_for_convergence)
-	   then {
-		CCTK_VInfo(CCTK_THORNSTRING,
-			   "==> finished (||Delta_h|| < %g)",
-			   double(solver_info.Delta_h_norm_for_convergence));
-		if (solver_info.final_H_update_if_exit_x_H_small)
-		   then {
-			CCTK_VInfo(CCTK_THORNSTRING,
-				   "    doing final H(h) evaluation");
-			horizon_function(ps, cgi, gi, false, true);
-			}
-		return true;				// *** NORMAL RETURN ***
-		}
-	}
-
-
-CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-	   "Newton_solve: no convergence in %d iterations!",
-	   solver_info.max_Newton_iterations);
-if (solver_info.final_H_update_if_exit_x_H_small)
-   then {
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "    doing final H(h) evaluation");
-	horizon_function(ps, cgi, gi, false, true);
-	}
-return false;						// *** ERROR RETURN ***
-}
diff --git a/src/gr/README b/src/gr/README
index 32b834f..fdd3450 100644
--- a/src/gr/README
+++ b/src/gr/README
@@ -3,16 +3,10 @@ general relativity.
 
 The main source files in this directory are as follows:
 
-driver.cc
-	Top-level driver for the whole thorn.
-
 horizon_function.cc
 	Computes the LHS of the AH equation, $H = H(h)$.
 
-horizon_Jacobian.cc
-	Computes the Jacobian of the $H(h)$ function, $J[H(h)]$.
-
-AHFinderDirect.hh
+gr.hh
 	Overall header file for all the external routines in this directory.
 
 gfns.hh
@@ -28,6 +22,9 @@ cg.hh
 doit.maple
 	Top-level driver for all the Maple code.
 
+maple.log
+	This is a log file (transcript) of the Maple computation.
+
 setup_gr_gfas.maple
 gr_gfas.minc
 	These files define the relativity gridfns for the Maple code,
@@ -39,5 +36,6 @@ horizon.maple
 	These files compute the $H(h)$ function and the $J[H(h)]$
 	coefficients in Maple.
 
-maple.log
-	This is a log file (transcript) of the Maple computation.
+makefile
+	This drives the Maple code-generation process; the resulting
+	(Maple-generated) C code lives in ../gr.cg/
diff --git a/src/gr/Schwarzschild_EF.cc b/src/gr/Schwarzschild_EF.cc
index 1bce005..628e641 100644
--- a/src/gr/Schwarzschild_EF.cc
+++ b/src/gr/Schwarzschild_EF.cc
@@ -39,7 +39,7 @@ using jtutil::error_exit;
 #include "../elliptic/Jacobian.hh"
 
 #include "gfns.hh"
-#include "AHFinderDirect.hh"
+#include "gr.hh"
 
 //******************************************************************************
 
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
deleted file mode 100644
index 5da7b3f..0000000
--- a/src/gr/driver.cc
+++ /dev/null
@@ -1,623 +0,0 @@
-// driver.cc -- top level driver for finding apparent horizons
-// $Id$
-//
-// <<<prototypes for functions local to this file>>>
-// AHFinderDirect_driver - top-level driver
-/// decode_Jacobian_method - decode the Jacobian_method parameter
-/// decode_geometry_method - decode the geometry_method parameter
-///
-/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
-/// setup_ellipsoid - setup up a coordiante ellipsoid in h
-///
-/// print_Jacobians - print a pair of Jacobians
-///
-
-#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 "stdc.h"
-#include "config.hh"
-#include "../jtutil/util.hh"
-#include "../jtutil/array.hh"
-#include "../jtutil/cpm_map.hh"
-#include "../jtutil/linear_map.hh"
-using jtutil::error_exit;
-
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
-
-#include "../elliptic/Jacobian.hh"
-
-#include "gfns.hh"
-#include "AHFinderDirect.hh"
-
-//******************************************************************************
-
-//
-// ***** prototypes for functions local to this file *****
-//
-
-namespace {
-enum Jacobian_method
-  decode_Jacobian_method(const char Jacobian_method_string[]);
-enum geometry_method
-  decode_geometry_method(const char geometry_method_string[]);
-
-void setup_Kerr_horizon(patch_system& ps,
-			fp x_posn, fp y_posn, fp z_posn,
-			fp m, fp a,
-			bool Kerr_Schild_flag);
-void setup_ellipsoid(patch_system& ps,
-		     fp x_center, fp y_center, fp z_center,
-		     fp x_radius, fp y_radius, fp z_radius);
-
-void print_Jacobians(const patch_system& ps,
-		     const Jacobian* Jac_NP,
-		     bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr,
-		     const char file_name[]);
-	  };
-
-//******************************************************************************
-
-//
-// This function is the Cactus interface for the test driver.
-//
-extern "C"
-  void AHFinderDirect_driver(CCTK_ARGUMENTS)
-{
-DECLARE_CCTK_ARGUMENTS
-DECLARE_CCTK_PARAMETERS
-
-CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures");
-
-
-//
-// set up the geometry parameters and the geometry interpolator
-//
-struct geometry_info gi;
-gi.geometry_method = decode_geometry_method(geometry_method);
-
-// parameters for geometry_method = "interpolate from Cactus grid"
-CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
-gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
-if (gi.operator_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "couldn't find interpolator \"%s\"!",
-		   geometry_interpolator_name);		/*NOTREACHED*/
-gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
-if (gi.param_table_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "bad geometry-interpolator parameter(s) \"%s\"!",
-		   geometry_interpolator_pars);		/*NOTREACHED*/
-
-// parameters for geometry_method = "Schwarzschild/EF"
-gi.geometry__Schwarzschild_EF__mass     = geometry__Schwarzschild_EF__mass;
-gi.geometry__Schwarzschild_EF__x_posn   = geometry__Schwarzschild_EF__x_posn;
-gi.geometry__Schwarzschild_EF__y_posn   = geometry__Schwarzschild_EF__y_posn;
-gi.geometry__Schwarzschild_EF__z_posn   = geometry__Schwarzschild_EF__z_posn;
-gi.geometry__Schwarzschild_EF__epsilon  = geometry__Schwarzschild_EF__epsilon;
-gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz;
-
-
-//
-// set up the interpatch interpolator
-//
-CCTK_VInfo(CCTK_THORNSTRING, "   setting up interpatch interpolator");
-const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name);
-if (interp_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "couldn't find interpolator \"%s\"!",
-		   interpatch_interpolator_name);		/*NOTREACHED*/
-const int interp_param_table_handle
-	= Util_TableCreateFromString(interpatch_interpolator_pars);
-if (interp_param_table_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "bad interpatch-interpolator parameter(s) \"%s\"!",
-		   interpatch_interpolator_pars);		/*NOTREACHED*/
-
-
-//
-// set up the Cactus grid info
-//
-CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
-struct cactus_grid_info cgi;
-cgi.GH = cctkGH;
-cgi.coord_origin[0] = cctk_origin_space[0];
-cgi.coord_origin[1] = cctk_origin_space[1];
-cgi.coord_origin[2] = cctk_origin_space[2];
-cgi.coord_delta[0] = cctk_delta_space[0];
-cgi.coord_delta[1] = cctk_delta_space[1];
-cgi.coord_delta[2] = cctk_delta_space[2];
-cgi.gridfn_dims[0] = cctk_lsh[0];
-cgi.gridfn_dims[1] = cctk_lsh[1];
-cgi.gridfn_dims[2] = cctk_lsh[2];
-// n.b. The  cgi.[gK]_dd_??_data  are actually  const fp *  pointers,
-//	since we won't modify the 3-D gridfn data!  But  static_cast<...>
-//      won't change const modifiers, so we just cast to  fp*  and let
-//	the assignment take care of the const part...
-cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxx"));
-cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy"));
-cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz"));
-cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy"));
-cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz"));
-cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz"));
-cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx"));
-cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy"));
-cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz"));
-cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy"));
-cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz"));
-cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));
-
-
-//
-// create the patch system and initialize the xyz derivative coefficients
-//
-patch_system ps(origin_x, origin_y, origin_z,
-		patch_system::type_of_name(patch_system_type),
-		N_ghost_points, N_overlap_points, delta_drho_dsigma,
-		gfns::nominal_min_gfn, gfns::nominal_max_gfn,
-		gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
-		interp_handle, interp_param_table_handle);
-
-
-//
-// set up the initial guess for the apparent horizon shape
-//
-if	(STRING_EQUAL(initial_guess_method, "read from file"))
-   then {
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "reading initial guess h from \"%s\"",
-		   h_file_name);
-	ps.read_ghosted_gridfn(gfns::gfn__h,
-			       h_file_name,
-			       false);		// no ghost zones
-	}
-   else {
-	if	(STRING_EQUAL(initial_guess_method, "ellipsoid"))
-	   then setup_ellipsoid(ps,
-				initial_guess__ellipsoid__x_center,
-				initial_guess__ellipsoid__y_center,
-				initial_guess__ellipsoid__z_center,
-				initial_guess__ellipsoid__x_radius,
-				initial_guess__ellipsoid__y_radius,
-				initial_guess__ellipsoid__z_radius);
-	else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
-	   then setup_Kerr_horizon(ps,
-				   initial_guess__Kerr__x_posn,
-				   initial_guess__Kerr__y_posn,
-				   initial_guess__Kerr__z_posn,
-				   initial_guess__Kerr__mass,
-				   initial_guess__Kerr__spin,
-				   false);	// use Kerr coords
-	else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
-	   then setup_Kerr_horizon(ps,
-				   initial_guess__Kerr__x_posn,
-				   initial_guess__Kerr__y_posn,
-				   initial_guess__Kerr__z_posn,
-				   initial_guess__Kerr__mass,
-				   initial_guess__Kerr__spin,
-				   true);	// use Kerr-Schild coords
-	else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "unknown initial_guess_method=\"%s\"!",
-			   initial_guess_method);		/*NOTREACHED*/
-
-	// write initial guess back to this file
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing initial guess h to \"%s\"",
-		   h_file_name);
-	ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
-					 true, gfns::gfn__h,
-					 h_file_name,
-					 false);	// no ghost zones
-	}
-
-
-//
-// decode/copy parameters into structures
-//
-
-struct Jacobian_info Jacobian_info;
-Jacobian_info.Jacobian_method = decode_Jacobian_method(Jacobian_method);
-Jacobian_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
-
-struct solver_info solver_info;
-solver_info.max_Newton_iterations = max_Newton_iterations;
-solver_info.output_h_and_H_at_each_Newton_iteration
-	= (output_h_and_H_at_each_Newton_iteration != 0);
-solver_info.h_file_name = h_file_name;
-solver_info.H_of_h_file_name = H_of_h_file_name;
-solver_info.H_norm_for_convergence = H_norm_for_convergence;
-solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
-solver_info.final_H_update_if_exit_x_H_small
-	= (final_H_update_if_exit_x_H_small != 0);
-
-//
-// find the apparent horizon
-//
-if      (STRING_EQUAL(method, "horizon function"))
-   then {
-	jtutil::norm<fp> H_norms;
-
-	const int timer_handle = CCTK_TimerCreateI();
-	CCTK_TimerStartI(timer_handle);
-	  horizon_function(ps, cgi, gi, false, true, &H_norms);
-	CCTK_TimerStopI(timer_handle);
-	printf("timer stats for evaluating H(h):\n");
-	CCTK_TimerPrintDataI(timer_handle, -1);
-
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   H(h) rms-norm %.2e, infinity-norm %.2e",
-		   H_norms.rms_norm(), H_norms.infinity_norm());
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing H(h) to \"%s\"",
-		   H_of_h_file_name);
-	ps.print_gridfn_with_xyz(gfns::gfn__H,
-				 true, gfns::gfn__h,
-				 H_of_h_file_name);
-	}
-else if (    STRING_EQUAL(method, "Jacobian test")
-	  || STRING_EQUAL(method, "Jacobian test (NP only)")    )
-   then {
-	const bool NP_only = STRING_EQUAL(method, "Jacobian test (NP only)");
-	const bool do_non_NP_Jac = ! NP_only;
-
-	// numerical perturbation
-	Jacobian* const Jac_NP = & new_Jacobian(ps, Jacobian_storage_method);
-	horizon_function(ps, cgi, gi, true);
-	Jacobian_info.Jacobian_method = Jacobian__numerical_perturb;
-	horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_NP);
-
-	// symbolic differentiation with finite diff d/dr
-	Jacobian* Jac_SD_FDdr = NULL;
-	if (do_non_NP_Jac)
-	   then {
-		Jac_SD_FDdr = & new_Jacobian(ps, Jacobian_storage_method);
-		horizon_function(ps, cgi, gi, true);
-		Jacobian_info.Jacobian_method
-			= Jacobian__symbolic_diff_with_FD_dr;
-		horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_SD_FDdr);
-		}
-
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing Jacobian%s to \"%s\"",
-		   (NP_only ? "" : "s"), Jacobian_file_name);
-	print_Jacobians(ps,
-			Jac_NP,
-			do_non_NP_Jac, Jac_SD_FDdr,
-			Jacobian_file_name);
-	}
-else if (STRING_EQUAL(method, "Newton solve"))
-   then {
-	Jacobian& Jac = new_Jacobian(ps, Jacobian_storage_method);
-
-	const int timer_handle = CCTK_TimerCreateI();
-	CCTK_TimerStartI(timer_handle);
-	  Newton_solve(ps, cgi, gi, Jacobian_info, solver_info, Jac);
-	CCTK_TimerStopI(timer_handle);
-	printf("timer stats for finding the apparent horizon:\n");
-	CCTK_TimerPrintDataI(timer_handle, -1);
-
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing final horizon shape h to \"%s\"",
-		   h_file_name);
-	ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
-					 true, gfns::gfn__h,
-					 h_file_name,
-					 false);	// no ghost zones
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing H(h) to \"%s\"",
-		   H_of_h_file_name);
-	ps.print_gridfn_with_xyz(gfns::gfn__H,
-				 true, gfns::gfn__h,
-				 H_of_h_file_name);
-	}
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "unknown method=\"%s\"!",
-		   method);					/*NOTREACHED*/
-}
-
-//******************************************************************************
-
-//
-// This function decodes the  Jacobian_method  parameter (string) into
-// an internal enum for future use.
-//
-namespace {
-enum Jacobian_method
-  decode_Jacobian_method(const char Jacobian_method_string[])
-{
-if	(STRING_EQUAL(Jacobian_method_string,
-		      "numerical perturbation"))
-   then return Jacobian__numerical_perturb;
-else if (STRING_EQUAL(Jacobian_method_string,
-		      "symbolic differentiation with finite diff d/dr"))
-   then return Jacobian__symbolic_diff_with_FD_dr;
-else if (STRING_EQUAL(Jacobian_method_string,
-		      "symbolic differentiation"))
-   then return Jacobian__symbolic_diff;
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "unknown Jacobian_method_string=\"%s\"!",
-		   Jacobian_method_string);			/*NOTREACHED*/
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function decodes the  geometry_method  parameter (string) into
-// an internal enum for future use.
-//
-namespace {
-enum geometry_method
-  decode_geometry_method(const char geometry_method_string[])
-{
-if	(STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid"))
-   then return geometry__interpolate_from_Cactus_grid;
-else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF"))
-   then return geometry__Schwarzschild_EF;
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "unknown geometry_method_string=\"%s\"!",
-		   geometry_method_string);			/*NOTREACHED*/
-}
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function sets up the horizon of a Kerr black hole in Kerr or
-// Kerr-Schild coordinates, on the nominal grid, in the  h  gridfn.
-//
-// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
-// and the horizon is worked out on page 13.2 of my AHFinderDirect notes.
-//
-// Arguments:
-// [xyz]_posn = The position of the Kerr black hole.
-// (m,a) = Describe the Kerr black hole.  Note that my convention has
-//	   a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a).
-// Kerr_Schild_flag = false to use Kerr coordinates,
-//		      true to use Kerr-Schild coordinates
-//
-namespace {
-void setup_Kerr_horizon(patch_system& ps,
-			fp x_posn, fp y_posn, fp z_posn,
-			fp m, fp a,
-			bool Kerr_Schild_flag)
-{
-const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr";
-
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "setting up horizon for Kerr in %s coords",
-	   name);
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "   posn=(%g,%g,%g) mass=%g spin=J/m^2=%g",
-	   double(x_posn), double(y_posn), double(z_posn),
-	   double(m), double(a));
-
-// horizon in Kerr coordinates is coordinate sphere
-const fp r = m * (1.0 + sqrt(1.0 - a*a));
-
-// horizon in Kerr-Schild coordinates is coordinate ellipsoid
-const fp  z_radius = r;
-const fp xy_radius = Kerr_Schild_flag ? r * sqrt(1.0 + a*a*m*m/(r*r)) : r;
-
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "   horizon is coordinate %s",
-	   Kerr_Schild_flag ? "ellipsoid" : "sphere");
-setup_ellipsoid(ps,
-		x_posn, y_posn, z_posn,
-		xy_radius, xy_radius, z_radius);
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function sets up an ellipsoid in the gridfn h, using the
-// formulas in "ellipsoid.maple" and the Maple-generated C code in
-// "ellipsoid.c":
-//
-// ellipsoid has center (A,B,C), radius (a,b,c)
-// angular coordinate system has center (U,V,W)
-//
-// direction cosines wrt angular coordinate center are (xcos,ycos,zcos)
-// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
-//
-// then the equation of the ellipsoid is
-//	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
-//	-----------------  +  ----------------  +  -----------------  =  1
-//	        a^2                  b^2                   c^2
-//
-// to solve this, we introduce intermediate variables
-//	AU = A - U
-//	BV = B - V
-//	CW = C - W
-//
-namespace {
-void setup_ellipsoid(patch_system& ps,
-		     fp x_center, fp y_center, fp z_center,
-		     fp x_radius, fp y_radius, fp z_radius)
-{
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "setting h = ellipsoid: center=(%g,%g,%g)",
-	   double(x_center), double(y_center), double(z_center));
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "                       radius=(%g,%g,%g)",
-	   double(x_radius), double(y_radius), double(z_radius));
-
-	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 xcos, ycos, zcos;
-		p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos);
-
-		// set up variables used by Maple-generated code
-		const fp AU = x_center - ps.origin_x();
-		const fp BV = y_center - ps.origin_y();
-		const fp CW = z_center - ps.origin_z();
-		const fp a = x_radius;
-		const fp b = y_radius;
-		const fp c = z_radius;
-
-		// compute the solutions r_plus and r_minus
-		fp r_plus, r_minus;
-		#include "ellipsoid.c"
-
-		// exactly one of the solutions (call it r) should be positive
-		fp r;
-		if      ((r_plus > 0.0) && (r_minus < 0.0))
-		   then r = r_plus;
-		else if ((r_plus < 0.0) && (r_minus > 0.0))
-		   then r = r_minus;
-		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-				   "\n"
-"   expected exactly one r>0 solution to quadratic, got 0 or 2!\n"
-"   %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
-"   direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n"
-"   ==> r_plus=%g r_minus=%g\n"
-				   ,
-				   p.name(), irho, isigma,
-				   double(rho), double(sigma),
-				   double(xcos), double(ycos), double(zcos),
-				   double(r_plus), double(r_minus));
-		   						/*NOTREACHED*/
-
-		// r = horizon radius at this grid point
-		p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = r;
-		}
-		}
-	}
-}
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function prints one or two Jacobian matrices (and their difference
-// in the latter case) to a named output file.
-//
-// Arguments:
-// do_non_NP_Jac = true if we should also print *Jac_SD_FDdr
-//		   false to skip this (and print only *Jac_NP)
-//
-namespace {
-void print_Jacobians(const patch_system& ps,
-		     const Jacobian* Jac_NP,
-		     bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr,
-		     const char file_name[])
-{
-FILE *fileptr = fopen(file_name, "w");
-if (fileptr == NULL)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "print_Jacobians(): can't open output file \"%s\"!",
-		   file_name);					/*NOTREACHED*/
-
-fprintf(fileptr, "# column 1 = x patch name\n");
-fprintf(fileptr, "# column 2 = x patch number\n");
-fprintf(fileptr, "# column 3 = x irho\n");
-fprintf(fileptr, "# column 4 = x isigma\n");
-fprintf(fileptr, "# column 5 = x II\n");
-fprintf(fileptr, "# column 6 = y patch name\n");
-fprintf(fileptr, "# column 7 = y patch number\n");
-fprintf(fileptr, "# column 8 = y irho\n");
-fprintf(fileptr, "# column 9 = y isigma\n");
-fprintf(fileptr, "# column 10 = y JJ\n");
-fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n");
-if (do_non_NP_Jac)
-   then {
-	fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n");
-	fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n");
-	fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n");
-	}
-
-    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
-    {
-    patch& xp = ps.ith_patch(xpn);
-
-    for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
-    {
-    for (int x_isigma = xp.min_isigma() ;
-	 x_isigma <= xp.max_isigma() ;
-	 ++x_isigma)
-    {
-    const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
-
-	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
-	{
-	patch& yp = ps.ith_patch(ypn);
-
-	for (int y_irho = yp.min_irho() ;
-	     y_irho <= yp.max_irho() ;
-	     ++y_irho)
-	{
-	for (int y_isigma = yp.min_isigma() ;
-	     y_isigma <= yp.max_isigma() ;
-	     ++y_isigma)
-	{
-	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
-
-	if (! Jac_NP->is_explicitly_stored(II,JJ))
-	   then continue;			// skip sparse points
-
-	const fp NP = (*Jac_NP)(II,JJ);
-	const fp SD_FDdr = do_non_NP_Jac ? (*Jac_SD_FDdr)(II,JJ) : 0.0;
-
-	if ((NP == 0.0) && (SD_FDdr == 0.0))
-	   then continue;		// skip zero values (if == )
-
-	fprintf(fileptr,
-		"%s %d %d %d %d\t%s %d %d %d %d\t%.10g",
-		xp.name(), xpn, x_irho, x_isigma, II,
-		yp.name(), ypn, y_irho, y_isigma, JJ,
-		double(NP));
-
-	if (do_non_NP_Jac)
-	   then {
-		const fp abs_NP      = jtutil::abs(NP     );
-		const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
-		const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
-		const fp SD_FDdr_abs_error = SD_FDdr - NP;
-		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;
-
-		fprintf(fileptr,
-			"\t%.10g\t%e\t%e",
-			double(SD_FDdr),
-			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
-		}
-
-	fprintf(fileptr, "\n");
-	}
-	}
-	}
-    }
-    }
-    }
-
-fclose(fileptr);
-}
-	  }
diff --git a/src/gr/ellipsoid.c b/src/gr/ellipsoid.c
deleted file mode 100644
index 49a6457..0000000
--- a/src/gr/ellipsoid.c
+++ /dev/null
@@ -1,28 +0,0 @@
-fp t1, t2, t3, t5, t6, t7, t9, t10, t12, t28;
-fp t30, t33, t35, t36, t40, t42, t43, t48, t49, t52;
-fp t55;
-      t1 = a*a;
-      t2 = b*b;
-      t3 = t1*t2;
-      t5 = t3*zcos*CW;
-      t6 = c*c;
-      t7 = t1*t6;
-      t9 = t7*ycos*BV;
-      t10 = t2*t6;
-      t12 = t10*xcos*AU;
-      t28 = xcos*xcos;
-      t30 = CW*CW;
-      t33 = BV*BV;
-      t35 = t10*t28;
-      t36 = ycos*ycos;
-      t40 = AU*AU;
-      t42 = t7*t36;
-      t43 = zcos*zcos;
-      t48 = t3*t43;
-      t49 = -2.0*t1*zcos*CW*ycos*BV-2.0*t2*zcos*CW*xcos*AU-2.0*t6*ycos*BV*xcos*
-AU+t2*t28*t30+t6*t28*t33-t35+t1*t36*t30+t6*t36*t40-t42+t1*t43*t33+t2*t43*t40-
-t48;
-      t52 = sqrt(-t3*t6*t49);
-      t55 = 1/(t35+t42+t48);
-      r_plus = (t5+t9+t12+t52)*t55;
-      r_minus = (t5+t9+t12-t52)*t55;
diff --git a/src/gr/ellipsoid.maple b/src/gr/ellipsoid.maple
deleted file mode 100644
index 349cc8d..0000000
--- a/src/gr/ellipsoid.maple
+++ /dev/null
@@ -1,37 +0,0 @@
-# ellipsoid.maple -- compute equations for offset ellipsoid setup
-# $Header$
-
-#
-# This program finds the intersection of an ellipsoid with an offset ray.
-#
-# ellipsoid has center (A,B,C), radius (a,b,c)
-# angular coordinate system has center (U,V,W)
-#
-# direction cosines wrt angular coordinate center are (alpha,beta,gamma)
-# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
-# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
-#
-# then the equation of the ellipsoid is
-#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
-#	-----------------  +  ----------------  +  -----------------  =  1
-#	        a^2                  b^2                   c^2
-#
-# to solve this, we introduce intermediate variables
-#	AU = A - U
-#	BV = B - V
-#	CW = C - W
-#
-eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
-
-read "../maple/util.mm";
-read "../maple/codegen2.mm";
-
-[solve(eqn, r)];
-map(simplify, %);
-[r_plus = %[1], r_minus = %[2]];
-solnlist := [codegen[optimize](%)];
-
-ftruncate("ellipsoid.c");
-print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]),
-		    "fp", "ellipsoid.c");
-codegen[C](solnlist, filename="ellipsoid.c");
diff --git a/src/gr/ellipsoid.out b/src/gr/ellipsoid.out
deleted file mode 100644
index 4a0ee01..0000000
--- a/src/gr/ellipsoid.out
+++ /dev/null
@@ -1,38 +0,0 @@
-    |\^/|     Maple 7 (IBM INTEL LINUX)
-._|\|   |/|_. Copyright (c) 2001 by Waterloo Maple Inc.
- \  MAPLE  /  All rights reserved. Maple is a registered trademark of
- <____ ____>  Waterloo Maple Inc.
-      |       Type ? for help.
-# ellipsoid.maple -- compute equations for offset ellipsoid setup
-# $Header: /numrelcvs/ThornburgCVS/AHFinderDirect/src/gr/ellipsoid.maple,v 1.2 2002/04/22 15:10:32 jthorn Exp $
-> 
-#
-# ellipsoid has center (A,B,C), radius (a,b,c)
-# angular coordinate system has center (U,V,W)
-#
-# direction cosines wrt angular coordinate center are (alpha,beta,gamma)
-# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
-# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
-#
-# then the equation of the ellipsoid is
-#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
-#	-----------------  +  ----------------  +  -----------------  =  1
-#	        a^2                  b^2                   c^2
-#
-# to solve this, we introduce intermediate variables
-#	AU = A - U
-#	BV = B - V
-#	CW = C - W
-#
-> eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
-                              2                2                2
-                 (xcos r - AU)    (ycos r - BV)    (zcos r - CW)
-          eqn := -------------- + -------------- + -------------- = 1
-                        2                2                2
-                       a                b                c
-
-> 
-> read "../maple/util.mm";
-Error, unable to read `../maple/util.mm`
-> quit
-bytes used=129844, alloc=196572, time=0.03
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
deleted file mode 100644
index 24a59a7..0000000
--- a/src/gr/horizon_Jacobian.cc
+++ /dev/null
@@ -1,422 +0,0 @@
-// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
-// $Id$
-//
-// <<<prototypes for functions local to this file>>>
-//
-// horizon_Jacobian - top-level driver to compute the Jacobian
-///
-/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
-///
-/// horizon_Jacobian_SD - compute the Jacobian (symbolic differentiation)
-/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
-///
-
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-#include <vector>
-
-#include "util_Table.h"
-#include "cctk.h"
-#include "cctk_Arguments.h"
-
-#include "stdc.h"
-#include "config.hh"
-#include "../jtutil/util.hh"
-#include "../jtutil/array.hh"
-#include "../jtutil/cpm_map.hh"
-#include "../jtutil/linear_map.hh"
-using jtutil::error_exit;
-
-#include "../util/coords.hh"
-#include "../util/grid.hh"
-#include "../util/fd_grid.hh"
-#include "../util/patch.hh"
-#include "../util/patch_edge.hh"
-#include "../util/patch_interp.hh"
-#include "../util/ghost_zone.hh"
-#include "../util/patch_system.hh"
-
-#include "../elliptic/Jacobian.hh"
-
-#include "gfns.hh"
-#include "AHFinderDirect.hh"
-
-//******************************************************************************
-
-//
-// ***** prototypes for functions local to this file *****
-//
-
-namespace {
-void horizon_Jacobian_NP(patch_system& ps,
-			 const struct cactus_grid_info& cgi,
-			 const struct geometry_info& ii,
-			 const struct Jacobian_info& Jacobian_info,
-			 Jacobian& Jac);
-
-void horizon_Jacobian_SD(patch_system& ps,
-			 const struct cactus_grid_info& cgi,
-			 const struct geometry_info& gi,
-			 const struct Jacobian_info& Jacobian_info,
-			 Jacobian& Jac);
-void add_ghost_zone_Jacobian(const patch_system& ps,
-			     const patch& xp, const ghost_zone& xmgz,
-			     int x_II,
-			     int xm_irho, int xm_isigma,
-			     fp mol, Jacobian& Jac);
-	  }
-
-//******************************************************************************
-
-//
-// This function is a top-level driver to compute the Jacobian matrix
-// J[H(h)] of the LHS function H(h).  It just decodes the Jacobian_method
-// parameter and calls the appropriate subfunction.
-void horizon_Jacobian(patch_system& ps,
-		      const struct cactus_grid_info& cgi,
-		      const struct geometry_info& gi,
-		      const struct Jacobian_info& Jacobian_info,
-		      Jacobian& Jac)
-{
-switch	(Jacobian_info.Jacobian_method)
-	{
-case Jacobian__numerical_perturb:
-	horizon_Jacobian_NP(ps, cgi, gi, Jacobian_info, Jac);
-	break;
-case Jacobian__symbolic_diff_with_FD_dr:
-	horizon_Jacobian_SD(ps, cgi, gi, Jacobian_info, Jac);
-	break;
-case Jacobian__symbolic_diff:
-	// FIXME FIXME: need true symbolic diff for d/dr terms
-	horizon_Jacobian_SD(ps, cgi, gi, Jacobian_info, Jac);
-	break;
-default:
-	error_exit(PANIC_EXIT,
-"***** horizon_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n"
-"                          (this should never happen!)\n"
-,
-		   int(Jacobian_info.Jacobian_method));		/*NOTREACHED*/
-	}
-}
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function computes the Jacobian matrix of the LHS function H(h)
-// by numerical perturbation of the H(h) function.  The algorithm is as
-// follows:
-//
-// we assume that H = H(h) has already been evaluated
-// save_H = H
-//	for each point (y,JJ)
-//	{
-//	const fp save_h_y = h at y;
-//	h at y += perturbation_amplitude;
-//	evaluate H(h) (silently)
-//		for each point (x,II)
-//		{
-//		Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
-//		}
-//	h at y = save_h_y;
-//	}
-// H = save_H
-//
-namespace {
-void horizon_Jacobian_NP(patch_system& ps,
-			 const struct cactus_grid_info& cgi,
-			 const struct geometry_info& ii,
-			 const struct Jacobian_info& Jacobian_info,
-			 Jacobian& Jac)
-{
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian (numerical perturbation)");
-const fp epsilon = Jacobian_info.perturbation_amplitude;
-
-ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
-
-	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
-	{
-	patch& yp = ps.ith_patch(ypn);
-	CCTK_VInfo(CCTK_THORNSTRING, "      perturbing in %s patch", yp.name());
-
-	for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho)
-	{
-	for (int y_isigma = yp.min_isigma() ;
-	     y_isigma <= yp.max_isigma() ;
-	     ++y_isigma)
-	{
-	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
-
-	const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
-	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
-	horizon_function(ps, cgi, ii);
-		for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
-		{
-		patch& xp = ps.ith_patch(xpn);
-
-		for (int x_irho = xp.min_irho() ;
-		     x_irho <= xp.max_irho() ;
-		     ++x_irho)
-		{
-		for (int x_isigma = xp.min_isigma() ;
-		     x_isigma <= xp.max_isigma() ;
-		     ++x_isigma)
-		{
-		const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
-		const fp old_H = xp.gridfn(gfns::gfn__save_H, x_irho,x_isigma);
-		const fp new_H = xp.gridfn(gfns::gfn__H, x_irho,x_isigma);
-		Jac(II,JJ) = (new_H - old_H) / epsilon;
-		}
-		}
-		}
-	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y;
-	}
-	}
-   	} 
-
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
-}
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function computes the Jacobian matrix of the LHS function H(h)
-// by symbolic differentiation from the Jacobian coefficient (angular)
-// gridfns.  The computation is done a Jacobian row at a time, using
-// equation (25) of my 1996 apparent horizon finding paper, except that
-// the d/dr term is done by numerical perturbation.
-//
-// The overall algorithm is
-//	save_H = H
-//	h += perturbation_amplitude
-//	evaluate H(h) (silently)
-//	Jac = diagonal[ (H(h) - save_H)/perturbation_amplitude ]
-//	h -= perturbation_amplitude
-//	H = save_H
-//	Jac += Jacobian coeffs * molecule coeffs
-//
-// Inputs (angular gridfns, on ghosted grid):
-//	h				# shape of trial surface
-//	H				# H(h) assumed to already be computed
-//	partial_H_wrt_partial_d_h	# Jacobian coefficients
-//	partial_H_wrt_partial_dd_h	# (also assumed to already be computed)
-//
-// Outputs:
-//	The Jacobian matrix is stored in the Jacobian object Jac.
-//
-namespace {
-void horizon_Jacobian_SD(patch_system& ps,
-			 const struct cactus_grid_info& cgi,
-			 const struct geometry_info& gi,
-			 const struct Jacobian_info& Jacobian_info,
-			 Jacobian& Jac)
-{
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian (symbolic differentiation)");
-
-ps.compute_synchronize_Jacobian();
-Jac.zero_matrix();
-
-
-//
-// compute d/dr term in Jacobian by numerical perturbation
-//
-
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      computing d/dr terms by numerical perturbation");
-const fp epsilon = Jacobian_info.perturbation_amplitude;
-ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
-ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
-horizon_function(ps, cgi, gi);
-	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 int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma);
-		const fp old_H = p.gridfn(gfns::gfn__save_H, irho,isigma);
-		const fp new_H = p.gridfn(gfns::gfn__H, irho,isigma);
-		Jac(II,II) += (new_H - old_H) / epsilon;
-		}
-		}
-	}
-ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
-
-
-//
-// now the main symbolic-differentiation Jacobian computation
-//
-
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      computing main terms by symbolic differentiation");
-    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
-    {
-    patch& xp = ps.ith_patch(xpn);
-
-	for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
-	{
-	for (int x_isigma = xp.min_isigma() ;
-	x_isigma <= xp.max_isigma() ;
-	++x_isigma)
-	{
-	//
-	// compute the main Jacobian terms for this grid point, i.e.
-	//	  partial H(this point x, Jacobian row II)
-	//	---------------------------------------------
-	//	partial h(other points y, Jacobian column JJ)
-	//
-
-	// Jacobian row index
-	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
-
-	// Jacobian coefficients for this point
-	const fp Jacobian_coeff_rho
-	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1,
-		       x_irho, x_isigma);
-	const fp Jacobian_coeff_sigma
-	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2,
-		       x_irho, x_isigma);
-	const fp Jacobian_coeff_rho_rho
-	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11,
-		       x_irho, x_isigma);
-	const fp Jacobian_coeff_rho_sigma
-	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12,
-		       x_irho, x_isigma);
-	const fp Jacobian_coeff_sigma_sigma
-	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22,
-		       x_irho, x_isigma);
-
-	// partial_rho, partial_rho_rho
-	    for (int m_irho = xp.molecule_min_m() ;
-		 m_irho <= xp.molecule_max_m() ;
-		 ++m_irho)
-	    {
-	    const int xm_irho = x_irho + m_irho;
-	    const fp Jac_rho     = Jacobian_coeff_rho
-				   * xp.partial_rho_coeff(m_irho);
-	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
-				   * xp.partial_rho_rho_coeff(m_irho);
-	    const fp Jac_sum = Jac_rho + Jac_rho_rho;
-	    if (xp.is_in_nominal_grid(xm_irho, x_isigma))
-	       then Jac(II, xp,xm_irho,x_isigma) += Jac_sum;
-	       else add_ghost_zone_Jacobian
-			(ps, xp, xp.minmax_rho_ghost_zone(m_irho < 0),
-			 II, xm_irho, x_isigma,
-			 Jac_sum, Jac);
-	    }
-
-	// partial_sigma, partial_sigma_sigma
-	    for (int m_isigma = xp.molecule_min_m() ;
-		 m_isigma <= xp.molecule_max_m() ;
-		 ++m_isigma)
-	    {
-	    const int xm_isigma = x_isigma + m_isigma;
-	    const fp Jac_sigma       = Jacobian_coeff_sigma
-				       * xp.partial_sigma_coeff(m_isigma);
-	    const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
-				       * xp.partial_sigma_sigma_coeff(m_isigma);
-	    const fp Jac_sum = Jac_sigma + Jac_sigma_sigma;
-	    if (xp.is_in_nominal_grid(x_irho, xm_isigma))
-	       then Jac(II, xp,x_irho,xm_isigma) += Jac_sum;
-	       else add_ghost_zone_Jacobian
-			(ps, xp, xp.minmax_sigma_ghost_zone(m_isigma < 0),
-			 II, x_irho, xm_isigma,
-			 Jac_sum, Jac);
-	    }
-
-	// partial_rho_sigma
-	    for (int m_irho = xp.molecule_min_m() ;
-		 m_irho <= xp.molecule_max_m() ;
-		 ++m_irho)
-	    {
-	    for (int m_isigma = xp.molecule_min_m() ;
-		 m_isigma <= xp.molecule_max_m() ;
-		 ++m_isigma)
-	    {
-	    const int xm_irho   = x_irho   + m_irho;
-	    const int xm_isigma = x_isigma + m_isigma;
-	    const fp Jac_rho_sigma
-	       = Jacobian_coeff_rho_sigma
-		 * xp.partial_rho_sigma_coeff(m_irho, m_isigma);
-	    if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
-	       then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma;
-	       else add_ghost_zone_Jacobian
-			(ps, xp,
-		 xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0,
-						       xm_irho, xm_isigma),
-			 II, xm_irho, xm_isigma,
-			 Jac_rho_sigma, Jac);
-	    }
-	    }
-
-	}
-	}
-    }
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function adds the ghost-zone Jacobian dependency contributions
-// for a single ghost-zone point, to a Jacobian matrix.
-//
-// Arguments:
-// ps = The patch system.
-// xp = The patch containing the center point of the molecule.
-// xmgz = If the x+m point is in a ghost zone, this must be that ghost zone.
-//	  If the x+m point is not in a ghost zone, this argument is ignored.
-// x_II = The Jacobian row of the x point.
-// xm_(irho,isigma) = The coordinates (in xp) of the x+m point of the molecule.
-// mol = The molecule coefficient.
-// Jac = The Jacobian matrix.
-//
-namespace {
-void add_ghost_zone_Jacobian(const patch_system& ps,
-			     const patch& xp, const ghost_zone& xmgz,
-			     int x_II,
-			     int xm_irho, int xm_isigma,
-			     fp mol, Jacobian& Jac)
-{
-const patch_edge& xme = xmgz.my_edge();
-const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma);
-const int xm_ipar  = xme. ipar_of_irho_isigma(xm_irho, xm_isigma);
-
-// FIXME: this won't change from one call to another
-//        ==> it would be more efficient to reuse the same buffer
-//            across multiple calls on this function
-int global_min_ym, global_max_ym;
-ps.synchronize_Jacobian_global_minmax_ym(global_min_ym, global_max_ym);
-jtutil::array1d<fp> Jacobian_buffer(global_min_ym, global_max_ym);
-
-// on what other points y does this molecule point xm depend
-// via the patch_system::synchronize() operation?
-int y_iperp;
-int y_posn, min_ym, max_ym;
-const patch_edge& ye = ps.synchronize_Jacobian(xmgz,
-					       xm_iperp, xm_ipar,
-					       y_iperp,
-					       y_posn, min_ym, max_ym,
-					       Jacobian_buffer);
-const patch& yp = ye.my_patch();
-
-// add the Jacobian contributions from the ym points
-	for (int ym = min_ym ; ym <= max_ym ; ++ym)
-	{
-	const int y_ipar = y_posn + ym;
-	const int y_irho   = ye.  irho_of_iperp_ipar(y_iperp,y_ipar);
-	const int y_isigma = ye.isigma_of_iperp_ipar(y_iperp,y_ipar);
-	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma);
-	Jac(x_II,JJ) += mol*Jacobian_buffer(ym);
-	}
-}
-	  }
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 41eaf23..3749ea0 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -39,7 +39,7 @@ using jtutil::pow4;
 #include "../elliptic/Jacobian.hh"
 
 #include "gfns.hh"
-#include "AHFinderDirect.hh"
+#include "gr.hh"
 
 //******************************************************************************
 
@@ -48,19 +48,19 @@ using jtutil::pow4;
 //
 
 namespace {
-void setup_xyz_posns(patch_system& ps, bool msg_flag);
-void interpolate_geometry(patch_system& ps,
+void setup_xyz_posns(patch_system& ps, bool print_msg_flag);
+bool interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
-			  bool msg_flag);
+			  bool print_msg_flag);
 void Schwarzschild_EF_geometry(patch_system& ps,
 			       const struct cactus_grid_info& cgi,
 			       const struct geometry_info& gi,
-			       bool msg_flag);
+			       bool print_msg_flag);
 void compute_H(patch_system& ps,
 	       bool Jacobian_flag,
 	       jtutil::norm<fp>* H_norms_ptr,
-	       bool msg_flag);
+	       bool print_msg_flag);
 	  }
 
 //******************************************************************************
@@ -100,37 +100,43 @@ void compute_H(patch_system& ps,
 // Arguments:
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
-// msg_flag = true to print status messages,
-//	      false to skip this.
+// print_msg_flag = true to print status messages,
+//		    false to skip this.
 // H_norms_ptr = (out) If this pointer is non-NULL, the norm object it
 //		       points to is updated with all the H values in the
 //		       grid.  This norm object can then be used to compute
 //		       various (gridwise) norms of H.
 //
-void horizon_function(patch_system& ps,
+// Results:
+// This function returns an error flag: normally this is false, but if
+// the geometry interpolation would need data outside the Cactus grid,
+// or data from an excised region, then this function returns true.
+//
+bool horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_info& gi,
 		      bool Jacobian_flag = false,
-		      bool msg_flag = false,
+		      bool print_msg_flag = false,
 		      jtutil::norm<fp>* H_norms_ptr = NULL)
 {
-if (msg_flag)
+if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
 // fill in values of all ghosted gridfns in ghost zones
 ps.synchronize();
 
 // set up xyz positions of grid points
-setup_xyz_posns(ps, msg_flag);
+setup_xyz_posns(ps, print_msg_flag);
 
 // compute the "geometry" g_ij, K_ij, and partial_k g_ij
 switch	(gi.geometry_method)
 	{
 case geometry__interpolate_from_Cactus_grid:
-	interpolate_geometry(ps, cgi, gi, msg_flag);
+	if (interpolate_geometry(ps, cgi, gi, print_msg_flag))
+	   then return true;				// *** ERROR RETURN ***
 	break;
 case geometry__Schwarzschild_EF:
-	Schwarzschild_EF_geometry(ps, gi, msg_flag);
+	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
 	break;
 default:
 	error_exit(PANIC_EXIT,
@@ -142,7 +148,9 @@ default:
 
 // compute remaining gridfns --> $H$ and optionally Jacobian coefficients
 // by algebraic ops and angular finite differencing
-compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag);
+compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag);
+
+return false;						// *** NORMAL RETURN ***
 }
 
 //******************************************************************************
@@ -152,9 +160,9 @@ compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag);
 // in the gridfns global_[xyz].  These will be used by interplate_geometry().
 //
 namespace {
-void setup_xyz_posns(patch_system& ps, bool msg_flag)
+void setup_xyz_posns(patch_system& ps, bool print_msg_flag)
 {
-if (msg_flag)
+if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      xyz positions and derivative coefficients");
 
@@ -213,13 +221,18 @@ if (msg_flag)
 // On the first call this function also modifies the interpolator
 // parameter table.
 //
+// Results:
+// This function returns an error flag: normally this is false, but if
+// the interpolation would need data outside the Cactus grid, or data
+// from an excised region, then this function returns true.
+//
 namespace {
-void interpolate_geometry(patch_system& ps,
+bool interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
-			  bool msg_flag)
+			  bool print_msg_flag)
 {
-if (msg_flag)
+if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      interpolating g_ij and Kij from Cactus 3-D grid");
 
@@ -343,7 +356,7 @@ if (first_time)
 	first_time = false;
 
 	// store derivative info in interpolator parameter table
-	if (msg_flag)
+	if (print_msg_flag)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "         setting up interpolator derivative info");
 
@@ -370,7 +383,7 @@ if (first_time)
 			   status);				/*NOTREACHED*/
 	}
 
-if (msg_flag)
+if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "         calling interpolator (%d points)",
 		   N_interp_points);
@@ -425,13 +438,13 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
 	assert((out_of_range_end == -1) || (out_of_range_end == +1));
 	const char end = (out_of_range_end == -1) ? '-' : '+';
 
-	error_exit(ERROR_EXIT,
-"***** interpolate_geometry():\n"
-"        the point (%g,%g,%g) on the trial horizon surface\n"
-"        is outside the grid in the %c%c direction!\n"
-,
-		   global_x, global_y, global_z,
-		   end, axis);					/*NOTREACHED*/
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "*** the trial-horizon-surface point (%g,%g,%g)",
+		   global_x, global_y, global_z);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "*** is outside the grid in the %c%c direction!",
+		   end, axis);
+	return true;					// *** ERROR RETURN ***
 	}
 if (status < 0)
    then error_exit(ERROR_EXIT,
@@ -439,6 +452,8 @@ if (status < 0)
 "        CCTK_InterpLocalUniform() status=%d\n"
 ,
 		   status);					/*NOTREACHED*/
+
+return false;						// *** NORMAL RETURN ***
 }
 	  }
 
@@ -458,9 +473,9 @@ namespace {
 void compute_H(patch_system& ps,
 	       bool Jacobian_flag,
 	       jtutil::norm<fp>* H_norms_ptr,
-	       bool msg_flag)
+	       bool print_msg_flag)
 {
-if (msg_flag)
+if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index ef5576b..0fcbabb 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -2,11 +2,8 @@
 # $Header$
 
 # Source files in this directory
-SRCS = driver.cc \
-       horizon_function.cc \
-       Schwarzschild_EF.cc \
-       horizon_Jacobian.cc \
-       Newton.cc
+SRCS = horizon_function.cc \
+       Schwarzschild_EF.cc
 
 # Subdirectories containing source files
 SUBDIRS =
diff --git a/src/gr/makefile b/src/gr/makefile
index 0b980fb..d3618f9 100644
--- a/src/gr/makefile
+++ b/src/gr/makefile
@@ -1,5 +1,5 @@
 # Makefile for GR code
-# $Id: makefile,v 1.6 2002-04-22 11:02:27 jthorn Exp $
+# $Id: makefile,v 1.7 2002-09-03 15:14:03 jthorn Exp $
 
 #
 # Environment Variables:
@@ -11,7 +11,6 @@
 #	cg		==> run all the *.mm files to generate C code 
 #	mm		==> preprocess all *.maple files to produce *.mm files
 #	clean		==> delete *.mm (mpp Maple preprocessor output)
-#	ellipsoid	==> compute ellipsoid equations
 #
 
 ifneq ($(MAPLE_VERSION),)
@@ -33,7 +32,3 @@ mm	: $(patsubst %.maple, %.mm, $(wildcard *.maple))
 .PHONY	: clean
 clean	:
 	-rm *.mm
-
-.PHONY	: ellipsoid
-ellipsoid:
-	maple <ellipsoid.maple >ellipsoid.out