* finish Jacobian support in src/util

* test driver for this in src/util/test_patch_system.cc
  --> IT WORKS!!!!!
* also fix comments on #include prerequisites to say ../jtutil/ not jt/


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

15 files changed, 418 insertions, 464 deletions

diff --git a/src/patch/coords.hh b/src/patch/coords.hh
index c0ae3e8..8543f9c 100644
--- a/src/patch/coords.hh
+++ b/src/patch/coords.hh
@@ -9,8 +9,10 @@
 
 //
 // prerequisites:
-//    "jt/util.hh"
-//    fp.hh
+//    <math.h>
+//    "stdc.h"
+//    "config.hh"
+//    "../jtutil/util.hh"
 //
 
 //******************************************************************************
diff --git a/src/patch/fd_grid.hh b/src/patch/fd_grid.hh
index e9a88c7..697def4 100644
--- a/src/patch/fd_grid.hh
+++ b/src/patch/fd_grid.hh
@@ -13,15 +13,14 @@
 // prerequisites:
 //    <stdio.h>
 //    <assert.h>
-//    <math.h>		// for M_PI (used by degree/radian conversions)
-//    "jt/util.hh"	// jtutil:: stuff:
-//			//    how_many_in_range(),
-//			//    degrees_of_radians(), radians_of_degrees(),
-//    "jt/array.hh"
-//    "jt/linear_map.hh"
-//    fp.hh
-//    coords.hh
-//    grid.hh
+//    <math.h>
+//    "stdc.h"
+//    "config.hh"
+//    "../jtutil/util.hh"
+//    "../jtutil/array.hh"
+//    "../jtutil/linear_map.hh"
+//    "coords.hh"
+//    "grid.hh"
 //
 
 //******************************************************************************
diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index ab40101..0b7825e 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -451,11 +451,24 @@ interp_result_buffer_
 
 //
 // construct the patch_interp object to interpolate from the *other* patch
-//
+// ... tell it to use corners if and only if the adjacent ghost zones
+//     are symmetry; cf  min_ipar()  and  max_ipar()  functions above,
+//     and header comments in this file
+//
+const bool ok_to_use_min_par_corner
+	= min_par_adjacent_ghost_zone().is_symmetry()
+	  ? true
+	  : false;
+const bool ok_to_use_max_par_corner
+	= max_par_adjacent_ghost_zone().is_symmetry()
+	  ? true
+	  : false;
 other_patch_interp_ = new patch_interp(other_edge(),
 				       min_other_iperp_, max_other_iperp_,
 				       *min_ipar_used_, *max_ipar_used_,
 				       *other_par_,
+				       ok_to_use_min_par_corner,
+				       ok_to_use_max_par_corner,
 				       interp_handle, interp_par_table_handle);
 }
 
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 6b81ee2..99da6a9 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -12,17 +12,18 @@
 //	<stdio.h>
 //	<assert.h>
 //	<math.h>
-//	"jt/util.hh"
-//	"jt/array.hh"
-//	"jt/cpm_map.hh"
-//	"jt/linear_map.hh"
-//	fp.hh
-//	coords.hh
-//	grid.hh
-//	fd_grid.hh
-//	patch.hh
-//	patch_edge.hh
-//	patch_interp.hh
+//	"stdc.h"
+//	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/cpm_map.hh"
+//	"../jtutil/linear_map.hh"
+//	"coords.hh"
+//	"grid.hh"
+//	"fd_grid.hh"
+//	"patch.hh"
+//	"patch_edge.hh"
+//	"patch_interp.hh"
 //
 
 //*****************************************************************************
@@ -65,27 +66,37 @@
 //   corner in the example below.  In this case we could get the gridfn
 //   data by either of two distinct interpolation operations (presumably
 //   from two distinct patches), which would in general give slightly
-//   different results.  In some ideal world we might take the average
-//   of these, but at present we split the corner down its diagonal.
-//   For the points on the diagonal we make an arbitrary choice; at
-//   present this is that they belong to (and get their data via) the
-//   rho ghost zone.
-//	[For elliptic equations it's important to ensure that
-//	each grid point get a value from exactly one computation,
-//	and in particular we can't allow two different ghost
-//	zones to "own" (assign values) to the same point, since
-//	then the first value would be "lost" and wouldn't show
-//	up in the Jacobian matrix.]
+//   different results.  In some ideal world we might do a centered
+//   interpolation using data from both patches, but this would be
+//   complicated:
+//   - it would require a 2-D interpolation
+//   - it would require bookkeeping for interpolating from multiple
+//     patches within the same ghost zone, indeed for the same ghost
+//     zone point
+//   At present, we follow a simpler approach: we split the corner down
+//   its diagonal,
+//	[for the points on the diagonal we make an arbitrary choice;
+//	at present this is that they belong to (and get their data via)
+//	the rho ghost zone.]
+//   and off-center the interpolation as necessary so each ghost-zone
+//   point gets data solely from the neighboring patch on its own side.
 // * A corner between a symmetry and an interpatch ghost zone, for
 //   example the +x/-y or -x/+y corners in the example below.  In this
-//   case we need to first do a symmetry operation in the neighboring
-//   patch so we don't have to off-center the interpolation in the
-//   non-corner part of the interpatch ghost zone.  Then after the
-//   interpatch interpolation, we need to do a final symmetry operation
-//   in this patch to set up gridfn data in the corner.
-//
-// To handle all these cases, we use a 3-phase algorithm to synchronize
-// ghost zones:
+//   case we first do a symmetry operation in the neighboring patch,
+//   then a fully centered interpolation (using the data just obtained
+//   from a symmetry operation) to get data in the non-corner part of
+//   the interpatch ghost zone.  After the interpatch interpolation,
+//   we do a final symmetry operation to get gridfn data in the corner.
+//
+// In general, then, a ghost zone is rhomboid-shaped: iperp lies in a
+// fixed interval, while ipar lies in an interval which may depend on
+// iperp.  In general, this shape depends on the type (symmetry vs interpatch)
+// of the adjacent ghost zones.
+//
+
+//
+// To properly handle all the symmetry/interpatch cases described above,
+// we use a 3-phase algorithm to synchronize ghost zones:
 // Phase 1: Fill in gridfn data at all the non-corner points of symmetry
 //	    ghost zones, by using the symmetries to get this data from
 //	    its "home patch" nominal grids.
@@ -171,6 +182,8 @@
 //   adjacent edges! There are no interpolation regions inside the -x or
 //   -y boundaries, since no interpolation is needed across those boundaries
 //   of this patch.
+// The diagonal *** line shows the boundary between the +x and +y ghost
+// zones.
 //
 // Our 3-phase algorithm described above thus becomes:
 // Phase 1: Fill in gridfn values at points marked with "s-x" below or
@@ -187,21 +200,6 @@
 //	    (yz plane) or -y boundary (xz plane), as described by the
 //	    +z patch's -x or -y  symmetry_ghost_zone  object respectively.
 //
-// The diagonal *** line shows the boundary between the +x and +y ghost
-// zones; points there may be interpolated via either of the two possible
-// interpatch_ghost_zone  objects.
-//	[At present we require all points in a given ghost zone
-//	to be interpolated from the same neighboring patch and
-//	 patch_interp  object, so must arbitrarily choose one
-//	of the two neighbors for the diagonal points.  In theory
-//	it would be better to take the average of the two neighbors,
-//	but in practice this doesn't matter for horizon finding
-//	or other elliptic stuff (it would matter for stability in
-//	time evolutions using multipatch finite differencing).]
-//
-// In general, then, a ghost zone is trapezoid-shaped: iperp lies in a
-// fixed interval, while ipar lies in an interval which may depend on iperp.
-//
 
 //*****************************************************************************
 
@@ -230,7 +228,7 @@ class	ghost_zone
 	{
 public:
 	//
-	// ***** main client interface *****
+	// ***** main high-level client interface *****
 	//
 	// "synchronize" a ghost zone, i.e. update the ghost-zone values
 	// of the specified gridfns via the appropriate sequence of
@@ -255,6 +253,12 @@ public:
 	// FIXME: should these be moved out into a separate Jacobian
 	//        object/class?
 	//
+	// Note that this function just computes the Jacobian of this
+	// ghost zone's  synchronize()  operation -- it does *NOT* take
+	// into account the 3-phase synchronization algorithm described
+	// in the header comments for this file.  (That's done by
+	//  patch_system::synchronize_ghost_zones_Jacobian() .)
+	//
 	// n.b. terminology is
 	//	partial gridfn at x
 	//	-------------------
@@ -365,7 +369,7 @@ public:
 		{ return my_edge().max_ipar_with_corners(); }
 
 	// actual min/max ipar in the ghost zone at a particular iperp
-	// (only defined in derived classes)
+	// (may depend on type of the adjacent ghost zones)
 	virtual int min_ipar(int iperp) const = 0;
 	virtual int max_ipar(int iperp) const = 0;
 
@@ -463,7 +467,7 @@ class	symmetry_ghost_zone
 	{
 public:
 	//
-	// ***** main client interface *****
+	// ***** main high-level client interface *****
 	//
 	// "synchronize" a ghost zone, i.e. update the ghost-zone values
 	// of the specified gridfns via the appropriate symmetry operations
@@ -601,7 +605,7 @@ class	interpatch_ghost_zone
 	{
 public:
 	//
-	// ***** main client interface *****
+	// ***** main high-level client interface *****
 	//
 	// "synchronize" a ghost zone, i.e. update the ghost-zone
 	// values of the specified gridfns via the appropriate
@@ -666,6 +670,9 @@ public:
 	fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const
 		{
 		assert(Jacobian_buffer_ != NULL);
+		// we're actually only supposed to be called with y_ipar_m
+		// in range, but as a backup we validate this here and return
+		// 0.0 for out of range
 		if ((y_ipar_m >= min_y_ipar_m_) && (y_ipar_m <= max_y_ipar_m_))
 		   then {
 			const int y_iperp = Jacobian_y_iperp(x_iperp);
@@ -710,8 +717,9 @@ private:
 		const;
 
 	// min/max ipar of the ghost zone for specified iperp
-	// taking into account our "triangular" corners
-	// (cf. the example at the start of this file)
+	// with possibly "triangular" corners depending on the type
+	// (symmetry vs interpatch) of the adjacent ghost zones
+	// (cf. comments & example at the start of this file)
 	int min_ipar(int iperp) const;
 	int max_ipar(int iperp) const;
 
diff --git a/src/patch/grid.hh b/src/patch/grid.hh
index 5d63698..aadec38 100644
--- a/src/patch/grid.hh
+++ b/src/patch/grid.hh
@@ -9,14 +9,13 @@
 // prerequisites:
 //    <stdio.h>
 //    <assert.h>
-//    <math.h>		// for M_PI (used by degree/radian conversions)
-//    "jt/util.hh"	// jtutil:: stuff:
-//			//    how_many_in_range(),
-//			//    degrees_of_radians(), radians_of_degrees(),
-//    "jt/array.hh"
-//    "jt/linear_map.hh"
-//    fp.hh
-//    coords.hh
+//    <math.h>
+//    "stdc.h"
+//    "config.hh"
+//    "../jtutil/util.hh"
+//    "../jtutil/array.hh"
+//    "../jtutil/linear_map.hh"
+//    "coords.hh"
 //
 
 //*****************************************************************************
diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn
index 1929236..148b738 100644
--- a/src/patch/make.code.defn
+++ b/src/patch/make.code.defn
@@ -9,6 +9,7 @@ SRCS = coords.cc \
        patch.cc \
        ghost_zone.cc \
        patch_interp.cc \
+       patch_info.cc \
        patch_system.cc \
        test_patch_system.cc
 
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index 8035a8e..dfa4369 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -18,9 +18,6 @@
 // patch::edge_adjacent_to_patch
 // patch::assert_all_ghost_zones_fully_setup
 //
-// patch_info::grid_array_pars
-// patch_info::grid_pars
-//
 
 #include <cstdio>
 #include <cmath>
@@ -413,111 +410,3 @@ max_rho_ghost_zone().assert_fully_setup();
 min_sigma_ghost_zone().assert_fully_setup();
 max_sigma_ghost_zone().assert_fully_setup();
 }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function computes, and returns a reference to, a
-//  struct grid_arrays::grid_array_pars  from the info in a
-//  struct patch_info  and the additional information in the arguments.
-//
-// The result refers to an internal static buffer in this function; the
-// usual caveats about lifetimes/overwriting apply.
-//
-// Arguments:
-// N_ghost_points = Width in grid points of all ghost zones.
-// N_extend_points = Number of grid points to extend each patch past
-//		     "just touching" so as to overlap neighboring patches.
-//		     Thus patches overlap by
-//			N_overlap_points = 2*N_extend_points + 1
-//		     grid points.  For example, with N_extend_points == 2,
-//		     here are the grid points of two neighboring patches:
-//			x   x   x   x   x   X   X
-//                                      |
-//			        O   O   o   o   o   o   o
-//		     Here | marks the "just touching" boundary,
-//		     x and o the grid points before this extension,
-//		     and X and O the extra grid points added by this
-//		     extension.
-// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
-//
-const grid_arrays::grid_array_pars&
-  patch_info::grid_array_pars(int N_ghost_points, int N_extend_points,
-			      fp delta_drho_dsigma)
-	const
-{
-static
-  struct grid_arrays::grid_array_pars grid_array_pars_buffer;
-
-grid_array_pars_buffer.min_irho
-	= jtutil::round<fp>::to_integer(min_drho  /delta_drho_dsigma);
-grid_array_pars_buffer.min_isigma
-	= jtutil::round<fp>::to_integer(min_dsigma/delta_drho_dsigma);
-grid_array_pars_buffer.max_irho
-	= grid_array_pars_buffer.min_irho
-	  + jtutil::round<fp>::to_integer(
-		(max_drho  -min_drho  ) / delta_drho_dsigma
-					 );
-grid_array_pars_buffer.max_isigma
-	= grid_array_pars_buffer.min_isigma
-	  + jtutil::round<fp>::to_integer(
-		(max_dsigma-min_dsigma) / delta_drho_dsigma
-					 );
-grid_array_pars_buffer.min_irho   -= N_extend_points;
-grid_array_pars_buffer.min_isigma -= N_extend_points;
-grid_array_pars_buffer.max_irho   += N_extend_points;
-grid_array_pars_buffer.max_isigma += N_extend_points;
-
-grid_array_pars_buffer.min_rho_N_ghost_points = N_ghost_points;
-grid_array_pars_buffer.max_rho_N_ghost_points = N_ghost_points;
-grid_array_pars_buffer.min_sigma_N_ghost_points = N_ghost_points;
-grid_array_pars_buffer.max_sigma_N_ghost_points = N_ghost_points;
-
-return grid_array_pars_buffer;
-}
-
-//******************************************************************************
-//
-//
-// This function computes, and returns a reference to, a
-//  struct grid_arrays::grid_pars  from the info in a  struct patch_info
-// and the additional information in the arguments.
-//
-// The result refers to an internal static buffer in this function; the
-// usual caveats about lifetimes/overwriting apply.
-//
-// Arguments:
-// N_extend_points = Number of grid points to extend each patch past
-//		     "just touching" so as to overlap neighboring patches.
-//		     Thus patches overlap by  2*N_extend_points + 1  grid
-//		     points.  For example, with N_extend_points == 2, here
-//		     are the grid points of two neighboring patches:
-//			x   x   x   x   x   X   X
-//                                      |
-//			        O   O   o   o   o   o   o
-//		     Here | marks the "just touching" boundary,
-//		     x and o the grid points before this extension,
-//		     and X and O the extra grid points added by this
-//		     extension.
-// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
-//
-const grid::grid_pars& patch_info::grid_pars(int N_extend_points,
-					     fp delta_drho_dsigma)
-	const
-{
-static
-  struct grid::grid_pars grid_pars_buffer;
-
-const fp extend_drho_dsigma = fp(N_extend_points) * delta_drho_dsigma;
-
-grid_pars_buffer.  min_drho   = min_drho   - extend_drho_dsigma;
-grid_pars_buffer.delta_drho   = delta_drho_dsigma;
-grid_pars_buffer.  max_drho   = max_drho   + extend_drho_dsigma;
-grid_pars_buffer.  min_dsigma = min_dsigma - extend_drho_dsigma;
-grid_pars_buffer.delta_dsigma = delta_drho_dsigma;
-grid_pars_buffer.  max_dsigma = max_dsigma + extend_drho_dsigma;
-
-return grid_pars_buffer;
-}
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 6f812e5..b58882d 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -3,7 +3,6 @@
 //
 // ***** how patch boundaries are handled *****
 // patch - abstract base class to describe a coordinate/grid patch
-// patch_info - POD struct with summary info for patch construction
 //
 // z_patch - derived class for a +/- z patch
 // x_patch - derived class for a +/- x patch
@@ -15,10 +14,11 @@
 //	<stdio.h>
 //	<assert.h>
 //	<math.h>
-//	"jt/util.hh"
-//	"jt/array.hh"
-//	"jt/linear_map.hh"
+//	"stdc.h"
 //	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/linear_map.hh"
 //	"coords.hh"
 //	"grid.hh"
 //	"fd_grid.hh"
@@ -425,61 +425,6 @@ private:
 //*****************************************************************************
 
 //
-// This (POD, and hence static-initializable) struct gives a minimal
-// set of information which varies from one patch to another.
-//
-// The member functions allow a patch to be constructed from this struct
-// together with the additional information provided by their arguments.
-// This doesn't allow the most general possible type of patch (eg it
-// constrains all ghost zones to have the same width), but does cover
-// the common casees.
-//
-// Arguments for member functions:
-// N_ghost_points = Width in grid points of all ghost zones.
-// N_extend_points = Number of grid points to extend each patch past
-//		     "just touching" so as to overlap neighboring patches.
-//		     Thus patches overlap by
-//			N_overlap_points = 2*N_extend_points + 1
-//		     grid points.  For example, with N_extend_points == 2,
-//		     here are the grid points of two neighboring patches:
-//			x   x   x   x   x   X   X
-//                                      |
-//			        O   O   o   o   o   o   o
-//		     Here | marks the "just touching" boundary,
-//		     x and o the grid points before this extension,
-//		     and X and O the extra grid points added by this
-//		     extension.
-// delta_drho_dsigma = Grid spacing (both rho and sigma) in degrees.
-//
-struct	patch_info
-	{
-	const char* name;
-	bool is_plus;
-	char ctype;
-	fp min_drho, max_drho;
-	fp min_dsigma, max_dsigma;
-
-	// compute and return reference to  struct grid_arrays::grid_array_pars
-	// ... result refers to internal static buffer;
-	//     the usual caveats about lifetimes/overwriting apply
-	const grid_arrays::grid_array_pars&
-	  grid_array_pars(int N_ghost_points, int N_extend_points,
-			  fp delta_drho_dsigma)
-		const;
-
-	// compute and return reference to  struct grid::grid_pars
-	// ... result refers to internal static buffer;
-	//     the usual caveats about lifetimes/overwriting apply
-	const grid::grid_pars& grid_pars(int N_extend_points,
-					 fp delta_drho_dsigma)
-		const;
-	};
-
-//*****************************************************************************
-//*****************************************************************************
-//*****************************************************************************
-
-//
 // This class describes a +/- z patch.  It doesn't define any new
 // functions not already present in  class patch ; it "just" defines
 // non-virtual versions of all the pure virtual functions defined there.
diff --git a/src/patch/patch_edge.hh b/src/patch/patch_edge.hh
index 58191b7..6140172 100644
--- a/src/patch/patch_edge.hh
+++ b/src/patch/patch_edge.hh
@@ -7,14 +7,15 @@
 //	<stdio.h>
 //	<assert.h>
 //	<math.h>
-//	"jt/util.hh"
-//	"jt/array.hh"
-//	"jt/linear_map.hh"
-//	fp.hh
-//	coords.hh
-//	grid.hh
-//	fd_grid.hh
-//	patch.hh
+//	"stdc.h"
+//	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/linear_map.hh"
+//	"coords.hh"
+//	"grid.hh"
+//	"fd_grid.hh"
+//	"patch.hh"
 //
 
 //*****************************************************************************
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 0c7181f..c19dd42 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -4,14 +4,13 @@
 //
 // patch_interp::patch_interp
 // patch_interp::~patch_interp
-// patch_interp::compute_[min,max]_ipar
 // patch_interp::interpolate
 //
-// patch_interp::query_interpolator
 // patch_interp::verify_Jacobian_sparsity_pattern_ok
 // patch_interp::molecule_minmax_ipar_m
 // patch_interp::molecule_posn
 // patch_interp::Jacobian
+// patch_interp::query_interpolator
 //
 
 #include <stdio.h>
@@ -56,22 +55,28 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 			   const jtutil::array1d<int>& min_parindex_array_in,
 			   const jtutil::array1d<int>& max_parindex_array_in,
 			   const jtutil::array2d<fp>& interp_par_in,
+			   bool ok_to_use_min_par_corner,
+			   bool ok_to_use_max_par_corner,
 			   int interp_handle_in, int interp_par_table_handle_in)
 	: my_patch_(my_edge_in.my_patch()),
 	  my_edge_(my_edge_in),
 	  min_gfn_(my_patch().ghosted_min_gfn()),
 	  max_gfn_(my_patch().ghosted_max_gfn()),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
-	  min_ipar_(compute_min_ipar()),
-	  max_ipar_(compute_max_ipar()),
+	  min_ipar_(ok_to_use_min_par_corner
+		    ? my_edge_in.min_ipar_with_corners()
+		    : my_edge_in.min_ipar_without_corners()),
+	  max_ipar_(ok_to_use_max_par_corner
+		    ? my_edge_in.max_ipar_with_corners()
+		    : my_edge_in.max_ipar_without_corners()),
 	  min_parindex_array_(min_parindex_array_in),
 	  max_parindex_array_(max_parindex_array_in),
 	  interp_par_(interp_par_in),
 	  interp_handle_(interp_handle_in),
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
-	  // note interpolator coordinate origin = par of min_ipar() point
+	  // note interpolator coordinate origin = min_ipar_
 	  // cf. comments in "patch_interp.hh" on gridfn array subscripting
-	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())),
+	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar_)),
 	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
 	  gridfn_type_codes_      (min_gfn_, max_gfn_),
 	  gridfn_data_ptrs_       (min_gfn_, max_gfn_),
@@ -132,35 +137,6 @@ Util_TableDestroy(interp_par_table_handle_);
 //*****************************************************************************
 
 //
-// These functions compute the [min,max] range of ipar in this patch
-// from which we may use gridfn data in interpolation.  This computation
-// depends on our adjacent ghost zones' types as described in the comments
-// to our constructor, but needs only my_patch_ and my_edge_ members of
-// this object to be valid.
-//
-int patch_interp::compute_min_ipar() const
-{
-const ghost_zone& min_par_adjacent_ghost_zone
-	= my_patch()
-	  .ghost_zone_on_edge(my_edge().min_par_adjacent_edge());
-return min_par_adjacent_ghost_zone.is_symmetry()
-       ? my_edge().min_ipar_with_corners()
-       : my_edge().min_ipar_without_corners();
-}
-
-int patch_interp::compute_max_ipar() const
-{
-const ghost_zone& max_par_adjacent_ghost_zone
-	= my_patch()
-	  .ghost_zone_on_edge(my_edge().max_par_adjacent_edge());
-return max_par_adjacent_ghost_zone.is_symmetry()
-       ? my_edge().max_ipar_with_corners()
-       : my_edge().max_ipar_without_corners();
-}
-
-//*****************************************************************************
-
-//
 // This function interpolates the specified range of ghosted gridfns
 // at all the coordinates specified when we were constructed.  (It does
 // a separate interpolator call for each iperp.)  This function stores
@@ -259,89 +235,6 @@ const CCTK_INT *gridfn_type_codes_ptr
 //******************************************************************************
 
 //
-// This function queries the interpolator with whatever arguments are
-// in the parameter table.  It specifies arguments such that no actual
-// interpolation is done.
-//
-// In particular, the following interpolator arguments are set up properly
-// based on  iperp :
-//	N_dims
-//	operator_handle, param_table_handle,
-//	coord_origin, coord_delta,
-//	N_interp_points, interp_coords_type_code, interp_coords,
-//	input_array_dims	# specifies the correct par size of the
-//				# patch interpolation region at this iperp
-// The following arguments are set to specify a single input array, but
-// with a NULL data pointer so no actual data is used:
-//	N_input_arrays
-//	input_array_type_codes, input_arrays
-// The following arguments are set to specify a single output array, but
-// with a NULL data pointer so no actual data is stored:
-//	N_output_arrays
-//	output_array_type_codes, output_arrays
-//
-// Arguments:
-// iperp = (in) Specifies where in the patch interpolation region the
-//		interpolator query should be done.
-//
-// Results:
-// If the interpolator returns an error status, this function does an
-//  error_exit()  (and does not return to its caller).  Otherwise, this
-// function returns the interpolator status code.
-//
-int patch_interp::query_interpolator(const char function_name[], int iperp)
-	const
-{
-int status;
-const int N_dims = 1;
-
-const int min_parindex = min_parindex_array_(iperp);
-const int max_parindex = max_parindex_array_(iperp);
-const int N_interp_points = jtutil::how_many_in_range(min_parindex,
-						      max_parindex);
-const int interp_coords_type_code = CCTK_VARIABLE_REAL;
-const void* const interp_coords[N_dims]
-	= { static_cast<const void *>(& interp_par_(iperp, min_parindex)) };
-
-const int N_input_arrays = 1;
-const CCTK_INT input_array_dims[N_dims]
-	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
-const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
-const void* const input_arrays[N_input_arrays] = { NULL };
-
-const int N_output_arrays = 1;
-const CCTK_INT output_array_type_codes[N_output_arrays]
-	= { CCTK_VARIABLE_REAL };
-void* const output_arrays[N_output_arrays] = { NULL };
-
-status = CCTK_InterpLocalUniform(N_dims,
-				 interp_handle_, interp_par_table_handle_,
-				 &gridfn_coord_origin_, &gridfn_coord_delta_,
-				 N_interp_points,
-				    interp_coords_type_code,
-				    interp_coords,
-				 N_input_arrays,
-				    input_array_dims,
-				    input_array_type_codes,
-				    input_arrays,
-				 N_output_arrays,
-				    output_array_type_codes,
-				    output_arrays);
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::query_interpolator():\n"
-"        on behalf of patch_interp::%s()\n"
-"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
-,
-		   function_name,
-		   status, iperp, min_iperp(), max_iperp());	/*NOTREACHED*/
-
-return status;
-}
-
-//******************************************************************************
-
-//
 // This function verifies that our interpolator has a Jacobian sparsity
 // pattern which we grok: at present this means that the molecules are
 // fixed-sized hypercubes, with size/shape independent of the interpolation
@@ -505,6 +398,7 @@ int status;
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
 	const int min_parindex = min_parindex_array_(iperp);
+	const int max_parindex = max_parindex_array_(iperp);
 
 	// set up the molecule-position query in the parameter table
 	CCTK_POINTER molecule_posn_ptrs[N_dims]
@@ -523,9 +417,21 @@ int status;
 			   status);				/*NOTREACHED*/
 
 	// query the interpolator to get the molecule positions
-	// for all parindex at this iperp (the interpolator will
-	// store them directly in the  posn_buffer(iperp,*)  array)
+	// for all parindex at this iperp; the interpolator will
+	// store the  parindex-min_parindex  values (cf comments
+	// on array subscripting at the start of "patch_interp.hh")
+	// directly in the  posn_buffer(iperp,*)  array)
 	query_interpolator("molecule_positions", iperp);
+
+	// convert the molecule positions from  parindex-min_ipar
+	// to  parindex  values (again, cf comments on array subscripting
+	// at the start of "patch_interp.hh")
+		for (int parindex = min_parindex ;
+		     parindex <= max_parindex ;
+		     ++parindex)
+		{
+		posn_buffer(iperp, parindex) += min_ipar();
+		}
 	}
 
 //
@@ -658,3 +564,86 @@ if ((status < 0) || (status1 < 0) || (status2 < 0))
 ,
 		   status, status1, status2);			/*NOTREACHED*/
 }
+
+//******************************************************************************
+
+//
+// This function queries the interpolator with whatever arguments are
+// in the parameter table.  It specifies arguments such that no actual
+// interpolation is done.
+//
+// In particular, the following interpolator arguments are set up properly
+// based on  iperp :
+//	N_dims
+//	operator_handle, param_table_handle,
+//	coord_origin, coord_delta,
+//	N_interp_points, interp_coords_type_code, interp_coords,
+//	input_array_dims	# specifies the correct par size of the
+//				# patch interpolation region at this iperp
+// The following arguments are set to specify a single input array, but
+// with a NULL data pointer so no actual data is used:
+//	N_input_arrays
+//	input_array_type_codes, input_arrays
+// The following arguments are set to specify a single output array, but
+// with a NULL data pointer so no actual data is stored:
+//	N_output_arrays
+//	output_array_type_codes, output_arrays
+//
+// Arguments:
+// iperp = (in) Specifies where in the patch interpolation region the
+//		interpolator query should be done.
+//
+// Results:
+// If the interpolator returns an error status, this function does an
+//  error_exit()  (and does not return to its caller).  Otherwise, this
+// function returns the interpolator status code.
+//
+int patch_interp::query_interpolator(const char function_name[], int iperp)
+	const
+{
+int status;
+const int N_dims = 1;
+
+const int min_parindex = min_parindex_array_(iperp);
+const int max_parindex = max_parindex_array_(iperp);
+const int N_interp_points = jtutil::how_many_in_range(min_parindex,
+						      max_parindex);
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const void* const interp_coords[N_dims]
+	= { static_cast<const void *>(& interp_par_(iperp, min_parindex)) };
+
+const int N_input_arrays = 1;
+const CCTK_INT input_array_dims[N_dims]
+	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
+const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
+const void* const input_arrays[N_input_arrays] = { NULL };
+
+const int N_output_arrays = 1;
+const CCTK_INT output_array_type_codes[N_output_arrays]
+	= { CCTK_VARIABLE_REAL };
+void* const output_arrays[N_output_arrays] = { NULL };
+
+status = CCTK_InterpLocalUniform(N_dims,
+				 interp_handle_, interp_par_table_handle_,
+				 &gridfn_coord_origin_, &gridfn_coord_delta_,
+				 N_interp_points,
+				    interp_coords_type_code,
+				    interp_coords,
+				 N_input_arrays,
+				    input_array_dims,
+				    input_array_type_codes,
+				    input_arrays,
+				 N_output_arrays,
+				    output_array_type_codes,
+				    output_arrays);
+if (status < 0)
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::query_interpolator():\n"
+"        on behalf of patch_interp::%s()\n"
+"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
+,
+		   function_name,
+		   status, iperp, min_iperp(), max_iperp());	/*NOTREACHED*/
+
+return status;
+}
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index b98d39f..5693a4a 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -1,4 +1,4 @@
-// patch_interp.hh -- patch interpolation region
+// patch_interp.hh -- interpolation from a patch
 // $Id$
 
 //
@@ -9,16 +9,15 @@
 //
 //	"cctk.h"
 //
-//	"jt/util.hh"
-//	"jt/array.hh"
-//	"jt/linear_map.hh"
-//	"jt/interpolate.hh"
-//	fp.hh
-//	coords.hh
-//	grid.hh
-//	fd_grid.hh
-//	patch.hh
-//	patch_edge.hh
+//	"stdc.h"
+//	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/linear_map.hh"
+//	"coords.hh"
+//	"grid.hh"
+//	"fd_grid.hh"
+//	"patch.hh"
 //
 
 //*****************************************************************************
@@ -36,17 +35,12 @@
 // patch's coordinates into our coordinates.
 //
 
+//
 // A patch_interp defines a "patch interpolation region", the region of
 // its owning patch from which this interpolation will use gridfn data.
 //
 
 //
-// Since all patch_interp member functions are const, a patch_interp
-// object is effectively always const whether declared that way or not,
-// so there's no harm in always declaring patch_interp objects as const.
-//
-
-//
 // The way the patch coordnates are constructed, any two adjacent patches
 // share a common (perpendicular) coordinate.  Thus we only have to do
 // 1-dimensional interpolation here (in the parallel direction).  In
@@ -94,20 +88,23 @@
 // the gridfn data is a bit tricky:  The interpolator accesses the
 // gridfn data using the generic (1-D) subscripting expression
 //	data[offset + i*stride]
-// where i is the data array index.  However, we don't want to use
-// offset, because it has to be supplied in the parameter table as
-// an array subscripted by gfn, and so would require changing the
-// parameter table for each call on interpolate() (with potentially
+// where  i  is the data array index.  However, we'd rather not use
+//  offset , because it has to be supplied in the parameter table as
+// an array subscripted by  gfn , and so would require changing the
+// parameter table for each call on  interpolate()  (with potentially
 // different numbers of gridfns being interpolated).  Instead, at each
-// iperp we use i = ipar-min_ipar, so the default offset=0 makes the
-// subscripting expression zero for ipar = min_ipar.  This also makes
-// the interpolator's min_i = 0 and max_i be dims-1 (both the defaults),
-// so those also don't have to be set in the parameter table.  We set
-// the interpolator's data coordinate origin to the par coordinate for
-// min_ipar, so it correctly maps i --> par.  With this strategy we can
-// share the interpolator parameter table across all the iperp values,
-// and we don't need to modify the parameter table at all after the
-// initial setup in our constructor.
+//  iperp  we use  i = ipar-min_ipar , so the default  offset=0  makes
+// the subscripting expression zero for  ipar = min_ipar .  This also
+// makes the interpolator's  min_i = 0  and  max_i  be  dims-1  (both
+// the defaults), so those also don't have to be set in the parameter
+// table either.  We set the interpolator's data coordinate origin to
+// the  par  coordinate for  min_ipar , so it correctly maps  i --> par .
+// With this strategy we can share the interpolator parameter table
+// across all the  iperp  values, and we don't need to modify the
+// parameter table at all after the initial setup in our constructor.
+// However, we do have to adjust the molecule positions in
+//  patch_interp::molecule_posn() , since the interpolator will return
+//  i  values, while  molecule_posn()  needs  ipar  values.
 //
 
 class	patch_interp
@@ -179,15 +176,7 @@ private:
 	int min_iperp() const { return min_iperp_; }
 	int max_iperp() const { return max_iperp_; }
 
-	// compute [min,max] ipar in this patch from which we may use
-	// gridfn data in interpolation -- depends on adjacent ghost zones'
-	// types as described in comments to constructor (below), but needs
-	// only my_patch_ and my_edge_ members of this object to be valid
-	int compute_min_ipar() const;
-	int compute_max_ipar() const;
-
-	// min/max (iperp,ipar) of the gridfn data we'll use
-	// (uses cached values set up by our constructor)
+	// min/max (iperp,ipar) of the gridfn data to use for interpolation
 	int min_ipar() const { return min_ipar_; }
 	int max_ipar() const { return max_ipar_; }
 
@@ -223,6 +212,9 @@ public:
 	//	  are unused.  We keep a reference to this array, so it
 	//	  should have a lifetime at last as long as that of this
 	//	  object.
+	// ok_to_use_[min,max]_par_corner
+	//	= Boolean flags saying whether or not we should use gridfn
+	//	  data from the [min,max]_par corners in the interpolation.
 	// interp_handle_in = Cactus handle to the interpolation operator.
 	// interp_par_table_handle_in
 	//	= Cactus handle to a Cactus key/value table giving
@@ -250,6 +242,8 @@ public:
 		     const jtutil::array1d<int>& min_parindex_array_in,
 		     const jtutil::array1d<int>& max_parindex_array_in,
 		     const jtutil::array2d<fp>& interp_par_in,
+		     bool ok_to_use_min_par_corner,
+		     bool ok_to_use_max_par_corner,
 		     int interp_handle_in, int interp_par_table_handle_in);
 	~patch_interp();
 
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index fb30ded..4a3e3d1 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -50,6 +50,7 @@ using jtutil::error_exit;
 #include "patch_edge.hh"
 #include "patch_interp.hh"
 #include "ghost_zone.hh"
+#include "patch_info.hh"
 #include "patch_system.hh"
 #include "patch_system_info.hh"
 
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index e00758f..8fd9541 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -9,16 +9,19 @@
 //	<stdio.h>
 //	<assert.h>
 //	<math.h>
-//	"jt/util.hh"
-//	"jt/array.hh"
-//	"jt/linear_map.hh"
-//	fp.hh
-//	coords.hh
-//	grid.hh
-//	patch.hh
-//	patch_edge.hh
-//	ghost_zone.hh
-//	patch_interp.hh
+//	"stdc.h"
+//	"config.hh"
+//	"../jtutil/util.hh"
+//	"../jtutil/array.hh"
+//	"../jtutil/linear_map.hh"
+//	"coords.hh"
+//	"grid.hh"
+//	"fd_grid.hh"
+//	"patch.hh"
+//	"patch_edge.hh"
+//	"ghost_zone.hh"
+//	"patch_interp.hh"
+//	"patch_info.hh"
 //
 
 //******************************************************************************
diff --git a/src/patch/patch_system_info.hh b/src/patch/patch_system_info.hh
index 72b0191..dd6cf9f 100644
--- a/src/patch/patch_system_info.hh
+++ b/src/patch/patch_system_info.hh
@@ -1,6 +1,19 @@
 // patch_system_info.hh -- static info describing various types of patch systems
 // $Id$
 
+//
+// prerequisites:
+//    <stdio.h>
+//    <assert.h>
+//    <math.h>
+//    "../jtutil/util.hh"
+//    "../jtutil//array.hh"
+//    "../jtutil/linear_map.hh"
+//    "coords.hh"
+//    "grid.hh"
+//    "patch_info.hh"
+//
+
 //******************************************************************************
 
 namespace patch_system_info
diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc
index 6371305..dbde22b 100644
--- a/src/patch/test_patch_system.cc
+++ b/src/patch/test_patch_system.cc
@@ -101,6 +101,7 @@ namespace {
 void test_synchronize_Jacobians(patch_system& ps,
 				int test_gfn, int NP_test_gfn,
 				fp perturbation_amplitude,
+				bool perturb_all_y_patch_points,
 				const char Jacobian_file_name[]);
 
 void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones);
@@ -203,6 +204,7 @@ else if (STRING_EQUAL(which_test, "synchronize Jacobian"))
    then test_synchronize_Jacobians(ps,
 				   test_fn_gfn, test_fn_copy_gfn,
 				   NP_Jacobian__perturbation_amplitude,
+				 (NP_Jacobian__perturb_all_y_patch_points != 0),
 				   Jacobian_file_name);
 
 else if (STRING_EQUAL(which_test, "derivatives"))
@@ -238,32 +240,51 @@ CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system");
 //	synchronize test_gfn
 //	print ths synchronized function
 //	set up the same test function in NP_test_gfn on the nominal grid
-//		for each patch p and ghost zone pgz
+//		for each patch xp and ghost zone xgz
 //		{
 //		compute the synchronize() Jacobian for this ghost zone
-//		   (via  pgz.compute_Jacobian()  et al)
-//			for each point x in (p,pgz)
+//		   (via  xgz.compute_Jacobian()  et al)
+//			for each point x in (p,xgz)
 //			{
-//				for each point y in gz.Jacobian_y_patch()
+//				for each point y in xgz.Jacobian_y_patch()
 //				{
+//				if ( !perturb_all_y_patch_points
+//				     && (y iperp != xgz Jacobian y iperp) )
+//				   then continue;	// *** LOOP CONTROL ***
 //				const fp save_y_gridfn = NP_test_gfn at y
 //				NP_test_gfn at y += perturbation_amplitude
 //				synchronize NP_test_gfn
 //				NP_Jacobian = (NP_test_gfn at x - test_gfn at x)
 //					      / perturbation_amplitude
 //				NP_test_gfn at y = save_y_gridfn
-//				Jacobian = pgz.Jacobian(...)
-//				if (m in molecule || the Jacobians differ)
-//				   then print Jacobian, NP_Jacobian,
-//					      Jacobian error
+//				Jacobian = xgz.Jacobian(...)
+//				print Jacobians to output file
+//				if (the Jacobians differ)
+//				   then {
+//					print lots of debugging info to stdout
+//					goto done;
+//					}
 //				}
 //			}
 //		}
+//	done:
+//
+// Note that this test can be *very* slow to run!  Eg. on a 731 MHz PIII
+// laptop, a full-sphere patch system with 5 degree resolution, a ghost zone
+// iperp width of 2 points, and perturb_all_y_patch_points set to true,
+// took around 25 minutes of cpu time, and produced an 20 MB Jacobian
+// output file.  Even with perturb_all_y_patch_points set to false, it
+// took 1 minute 15 seconds cpu time, and produced a 1 MB Jacobian output
+// file.
 //
 // Arguments:
 // ps = (in out) THe patch system in/on which to do the computations.
 // test_gfn, NP_test_gfn = The gfns of two ghosted test gridfns.
 // perturbation_amplitude = The perturbation amplitude for the NP Jacobian.
+// perturb_all_y_patch_points
+//	= true  ==> When computing the NP Jacobian, try perturbing at
+//		    every point in the y patch.  This gives 
+//	  false ==> Only try perturbing at y_iper == Jacobian_y_iperp.
 // Jacobian_file_name = The name of the output file to which both Jacobians
 //			should be written.
 //
@@ -271,8 +292,11 @@ namespace {
 void test_synchronize_Jacobians(patch_system& ps,
 				int test_gfn, int NP_test_gfn,
 				fp perturbation_amplitude,
+				bool perturb_all_y_patch_points,
 				const char Jacobian_file_name[])
 {
+CCTK_VInfo(CCTK_THORNSTRING, "testing synchronize() Jacobian...");
+
 setup_sym_fn_xyz(ps, test_gfn, false);
 ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn);
 ps.print_ghosted_gridfn(test_fn_gfn, "test_fn.dat");
@@ -302,10 +326,10 @@ fprintf(fileptr, "# column 15 = Jacobian\n");
 fprintf(fileptr, "# column 16 = NP_Jacobian\n");
 fprintf(fileptr, "# column 17 = Jacobian error\n");
 
-    //*** for each patch p and ghost zone pgz
-    for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+    //*** for each patch p and ghost zone xgz
+    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
     {
-    patch& p = ps.ith_patch(pn);
+    patch& xp = ps.ith_patch(xpn);
 
     // n.b. these loops must use _int_ variables for the loop
     //      to terminate!
@@ -313,72 +337,144 @@ fprintf(fileptr, "# column 17 = Jacobian error\n");
     {
     for (int want_rho = false ; want_rho <= true ; ++want_rho)
     {
-    const patch_edge& pe = p.minmax_ang_patch_edge(want_min, want_rho);
-    ghost_zone& pgz      = p.minmax_ang_ghost_zone(want_min, want_rho);
+    const patch_edge& xe = xp.minmax_ang_patch_edge(want_min, want_rho);
+    ghost_zone& xgz      = xp.minmax_ang_ghost_zone(want_min, want_rho);
+
+    patch&            yp = xgz.Jacobian_y_patch();
+    const patch_edge& ye = xgz.Jacobian_y_edge();
 
-    pgz.compute_Jacobian(test_gfn, test_gfn,
+    CCTK_VInfo(CCTK_THORNSTRING,
+	       "   testing x patch %s, edge %s",
+	       xp.name(), xe.name());
+    CCTK_VInfo(CCTK_THORNSTRING,
+	       "           y patch %s, edge %s",
+	       yp.name(), ye.name());
+
+    xgz.compute_Jacobian(test_gfn, test_gfn,
 			 true, true, true);	// want *all* of ghost zone
-    const int Jacobian_y_min_ipar_m = pgz.Jacobian_y_min_ipar_m();
-    const int Jacobian_y_max_ipar_m = pgz.Jacobian_y_max_ipar_m();
+    const int Jacobian_y_min_ipar_m = xgz.Jacobian_y_min_ipar_m();
+    const int Jacobian_y_max_ipar_m = xgz.Jacobian_y_max_ipar_m();
 
-	//*** for each point x in (p,pgz)
-	for (int x_iperp = pgz.min_iperp() ;
-	     x_iperp <= pgz.max_iperp() ;
+	//*** for each point x in (p,xgz)
+	for (int x_iperp = xgz.min_iperp() ;
+	     x_iperp <= xgz.max_iperp() ;
 	     ++x_iperp)
 	{
-	for (int x_ipar = pgz.min_ipar(x_iperp) ;
-	     x_ipar <= pgz.max_ipar(x_iperp) ;
+	for (int x_ipar = xgz.min_ipar(x_iperp) ;
+	     x_ipar <= xgz.max_ipar(x_iperp) ;
 	     ++x_ipar)
 	{
-	const int x_irho   = pe.  irho_of_iperp_ipar(x_iperp, x_ipar);
-	const int x_isigma = pe.isigma_of_iperp_ipar(x_iperp, x_ipar);
+	const int x_irho   = xe.  irho_of_iperp_ipar(x_iperp, x_ipar);
+	const int x_isigma = xe.isigma_of_iperp_ipar(x_iperp, x_ipar);
 
-	patch&            q  = pgz.Jacobian_y_patch();
-	const patch_edge& qe = pgz.Jacobian_y_edge();
-	const int Jacobian_y_iperp = pgz.Jacobian_y_iperp(x_iperp);
+	const int Jacobian_y_iperp = xgz.Jacobian_y_iperp(x_iperp);
 	const int Jacobian_y_ipar_posn
-		= pgz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
+		= xgz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
 
 	    //*** for each point y in gz.Jacobian_patch()
-	    for (int y_irho = q.min_irho() ; y_irho <= q.max_irho() ; ++y_irho)
+	    for (int y_irho = yp.min_irho() ;
+		     y_irho <= yp.max_irho() ;
+		     ++y_irho)
 	    {
-	    for (int y_isigma = q.min_isigma() ;
-		 y_isigma <= q.max_isigma() ;
+	    for (int y_isigma = yp.min_isigma() ;
+		 y_isigma <= yp.max_isigma() ;
 		 ++y_isigma)
 	    {
+	    const int y_iperp = ye.iperp_of_irho_isigma(y_irho,y_isigma);
+	    const int y_ipar  = ye. ipar_of_irho_isigma(y_irho,y_isigma);
+
+	    if ( !perturb_all_y_patch_points && (y_iperp != Jacobian_y_iperp) )
+	       then continue;				// *** LOOP CONTROL ***
+
 	    // compute the NP Jacobian
 	    const fp save_y_gridfn
-		= q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma);
-	    q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma)
+		= yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma);
+	    yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma)
 		+= perturbation_amplitude;
 
 	    ps.synchronize_ghost_zones(NP_test_gfn, NP_test_gfn);
 	    const fp NP_Jacobian
-		= (   p.ghosted_gridfn(NP_test_gfn, x_irho,x_isigma)
-		    - p.ghosted_gridfn(   test_gfn, x_irho,x_isigma) )
+		= (   xp.ghosted_gridfn(NP_test_gfn, x_irho,x_isigma)
+		    - xp.ghosted_gridfn(   test_gfn, x_irho,x_isigma) )
 		  / perturbation_amplitude;
-	    q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) = save_y_gridfn;
+	    yp.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) = save_y_gridfn;
 
 	    // compute the query Jacobian
-	    const int y_iperp = qe.iperp_of_irho_isigma(y_irho,y_isigma);
-	    const int y_ipar  = qe. ipar_of_irho_isigma(y_irho,y_isigma);
 	    const int y_ipar_m = y_ipar - Jacobian_y_ipar_posn;
 	    const bool m_in_molecule
 		= (y_iperp == Jacobian_y_iperp)
 		  && (y_ipar_m >= Jacobian_y_min_ipar_m)
 		  && (y_ipar_m <= Jacobian_y_max_ipar_m);
-	    const fp Jacobian = pgz.Jacobian(x_iperp, x_ipar, y_ipar_m);
+	    const fp Jacobian = m_in_molecule
+				? xgz.Jacobian(x_iperp, x_ipar, y_ipar_m)
+				: 0.0;
 
 	    // print the results
-	    if (m_in_molecule || (Jacobian != NP_Jacobian))
-	       then fprintf(fileptr,
+	    const fp error = Jacobian - NP_Jacobian;
+	    fprintf(fileptr,
 "%d %d %d\t%d %d\t%d %d\t%d %d %d\t%d %d\t%d %d\t%.10g\t%.10g\t%e\n",
-			    p.patch_number(), pe.is_min(), pe.is_rho(),
-			    x_iperp, x_ipar, x_irho, x_isigma,
-			    q.patch_number(), qe.is_min(), qe.is_rho(),
-			    y_iperp, y_ipar, y_irho, y_isigma,
-			    double(Jacobian), double(NP_Jacobian),
-			    double(Jacobian-NP_Jacobian));
+		    xp.patch_number(), xe.is_min(), xe.is_rho(),
+		    x_iperp, x_ipar, x_irho, x_isigma,
+		    yp.patch_number(), ye.is_min(), ye.is_rho(),
+		    y_iperp, y_ipar, y_irho, y_isigma,
+		    double(Jacobian), double(NP_Jacobian),
+		    double(error));
+
+	    // debugging code in case the Jacobian is wrong :(
+	    if (jtutil::abs(error) > 1.0e-6)
+	       then {
+		    printf("### large Jacobian error!\n");
+
+		    printf("x: p patch %s, edge %s\n", xp.name(), xe.name());
+		    printf("y: q patch %s, edge %s\n", yp.name(), ye.name());
+
+		    const fp x_rho   = xp.rho_of_irho    (x_irho);
+		    const fp x_sigma = xp.sigma_of_isigma(x_isigma);
+		    const fp y_rho   = yp.rho_of_irho    (y_irho);
+		    const fp y_sigma = yp.sigma_of_isigma(y_isigma);
+
+		    const fp x_drho   = jtutil::degrees_of_radians(x_rho);
+		    const fp x_dsigma = jtutil::degrees_of_radians(x_sigma);
+		    const fp y_drho   = jtutil::degrees_of_radians(y_rho);
+		    const fp y_dsigma = jtutil::degrees_of_radians(y_sigma);
+
+		    printf(
+"x iperp=%d ipar=%d   irho=%d isigma=%d   drho=%g dsigma=%g\n",
+			   x_iperp, x_ipar, x_irho, x_isigma, x_drho, x_dsigma);
+		    printf(
+"y iperp=%d ipar=%d   irho=%d isigma=%d   drho=%g dsigma=%g\n",
+			   y_iperp, y_ipar, y_irho, y_isigma, y_drho, y_dsigma);
+
+		    const fp x_mu  = xp.mu_of_rho_sigma (x_rho, x_sigma);
+		    const fp x_nu  = xp.nu_of_rho_sigma (x_rho, x_sigma);
+		    const fp x_phi = xp.phi_of_rho_sigma(x_rho, x_sigma);
+		    const fp y_mu  = yp.mu_of_rho_sigma (y_rho, y_sigma);
+		    const fp y_nu  = yp.nu_of_rho_sigma (y_rho, y_sigma);
+		    const fp y_phi = yp.phi_of_rho_sigma(y_rho, y_sigma);
+
+		    const fp x_dmu  = jtutil::degrees_of_radians(x_mu);
+		    const fp x_dnu  = jtutil::degrees_of_radians(x_nu);
+		    const fp x_dphi = jtutil::degrees_of_radians(x_phi);
+		    const fp y_dmu  = jtutil::degrees_of_radians(y_mu);
+		    const fp y_dnu  = jtutil::degrees_of_radians(y_nu);
+		    const fp y_dphi = jtutil::degrees_of_radians(y_phi);
+		
+		    printf("x dmu=%g dnu=%g dphi=%g\n", x_dmu, x_dnu, x_dphi);
+		    printf("y dmu=%g dnu=%g dphi=%g\n", y_dmu, y_dnu, y_dphi);
+
+		    printf("Jacobian=%.10g\tNP_Jacobian=%.10g\terror=%e\n",
+			   double(Jacobian), double(NP_Jacobian),
+			   double(error));
+
+		    printf("Jacobian_y_[min,max]_ipar_m=[%d,%d]\n",
+			   Jacobian_y_min_ipar_m, Jacobian_y_max_ipar_m);
+		    printf("Jacobian_y_iperp=%d Jacobian_y_ipar_posn=%d\n",
+			   Jacobian_y_iperp, Jacobian_y_ipar_posn);
+		    printf("y_ipar_m=%d\n", y_ipar_m);
+
+		    printf("###\n");
+		    goto done;
+		    }
 	    }
 	    }
 	}
@@ -388,6 +484,7 @@ fprintf(fileptr, "# column 17 = Jacobian error\n");
 
     }
 
+done:
 fclose(fileptr);
 }
 	  }