fix subscripting of gridfn arrays by the interpolator

--> now set offset as well as stride and min/max subscript
... also add comments clarifying this


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

1 files changed, 20 insertions, 0 deletions

diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 0708b9f..a213b43 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -85,6 +85,22 @@
 // call.  [Different iperp values involve different sets of (1-D)
 // gridfn data, and so inherently require distinct interpolator calls.]
 //
+// Setting up the array subscripting for the interpolator to access
+// 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.  It's convenient to take i=ipar
+// (this way the mapping from i to par is just given by the usual
+// origin and delta values for the par coordinate).  The stride is
+// fixed by the gridfn storage mapping (we get this from the patch).
+// We then calculate the offset such that  offset + min_ipar*stride == 0 ,
+// and then pass the data pointer for each interpolator call as a pointer
+// to the min_ipar grid point at that iperp.  (We also need to set the
+// min and max i in the interpolator parameter table.)  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 after the
+// initial setup.
+//
 
 class	patch_interp
 	{
@@ -183,6 +199,10 @@ private:
 	const patch& my_patch_;
 	const patch_edge& my_edge_;
 
+	// range of gfn we can handle
+	// (any given interpolate() call may specify a subrange)
+	const int min_gfn_, max_gfn_;
+
 	// patch interpolation region,
 	// i.e. range of (iperp,ipar) in this patch from which
 	// we will use gridfn data in interpolation