add some new fns to determine what angular patch contains

the ray from the origin out to a given (x,y,z)
--> these will be needed to suppot checking whether an arbitrary (x,y,z)
is inside/outside the horizon, and also to support moving the origin point
to follow a moving BH


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

2 files changed, 100 insertions, 3 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index b80f865..b526c32 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -18,8 +18,11 @@
 // patch_system::N_patches_of_type
 // patch_system::name_of_type
 // patch_system::type_of_name
+// patch_system::plus_or_minus_xyz_patch
 // patch_system::patch_number_of_name
 //
+// patch_system::patch_containing_local_xyz
+//
 // patch_system::patch_irho_isigma_of_gpn
 // patch_system::ghosted_patch_irho_isigma_of_gpn
 //
@@ -1026,6 +1029,34 @@ else	error_exit(PANIC_EXIT,
 //******************************************************************************
 
 //
+// This function finds a (the) patch with a specified sign and xyz ctype.
+// If no such patch exists, it does an  error_exit()  (and doesn't return
+// to the caller).
+//
+// Bugs:
+// - This function could be implemented to be very fast, but right now
+//   it just does a sequential search through all the patches. :(
+//
+const patch& patch_system::plus_or_minus_xyz_patch(bool is_plus, char ctype)
+	const
+{
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	const patch& p = ith_patch(pn);
+	if ((p.is_plus() == is_plus) && (p.ctype() == ctype))
+	   then return p;
+	}
+
+error_exit(ERROR_EXIT,
+"***** patch_system::plus_or_minus_xyz_patch():\n"
+"        can't find any %c%c patch!"
+,
+	   (is_plus ? '+' : '-'), ctype);			/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
 // This function finds a patch from its human-readable name, and returns
 // the patch number, or does an error_exit() if no patch is found with
 // the specified name.
@@ -1049,6 +1080,57 @@ error_exit(ERROR_EXIT,
 //******************************************************************************
 
 //
+// This function finds what patch contains (the ray from the origin to)
+// a given local (x,y,z) position.
+//
+// If there are multiple patches containing the position, we return the
+// one which would still contain it if patches didn't overlap; if multiple
+// patches satisfy this criterion then it's arbitrary which one we return.
+//
+// If no patch contains the position (this can only happen for a
+// non--full-sphere patch system), then we do an  error_exit()  (and
+// don't return to the caller).
+//
+// Arguments:
+// (x,y,z) = The local coordinates to be converted.
+//
+// Results:
+// This function returns a reference to the containing patch.
+//
+const patch& patch_system::patch_containing_local_xyz(fp x, fp y, fp z)
+	const
+{
+if ((x == 0.0) && (y == 0.0) && (z == 0.0))
+   then error_exit(ERROR_EXIT,
+"***** patch_system::patch_containing_local_xyz(): position is at the origin!");
+								/*NOTREACHED*/
+
+// to which axis is (x,y,z) closest?
+// ... or equivalently, which of |x|, |y|, and |z| is largest?
+const fp abs_x = jtutil::abs(x);
+const fp abs_y = jtutil::abs(y);
+const fp abs_z = jtutil::abs(z);
+
+if	((abs_z >= abs_x) && (abs_z >= abs_y))
+   then return plus_or_minus_xyz_patch(z > 0.0, 'z');	// +/- z patch
+else if ((abs_x >= abs_y) && (abs_x >= abs_z))
+   then return plus_or_minus_xyz_patch(x > 0.0, 'x');	// +/- x patch
+else if ((abs_y >= abs_x) && (abs_y >= abs_z))
+   then return plus_or_minus_xyz_patch(y > 0.0, 'y');	// +/- y patch
+else	error_exit(ERROR_EXIT,
+"***** patch_system::patch_containing_local_xyz():\n"
+"        unknown (wierd!) ordering of |x|, |y|, and |z|!\n"
+"        (this should never happen!)\n"
+"        (x,y,z)=(%g,%g,%g)\n"
+,
+		   double(x), double(y), double(z));
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function decodes a 0-origin grid point number into a
 // (patch,irho,isigma) triple.
 //
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 6eb70d6..08c15d8 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -87,11 +87,22 @@ public:
 	fp global_z_of_local_z(fp local_z) const
 		{ return global_coords_.global_z_of_local_z(local_z); }
 
-	// ... get global (x,y,z) coordinates of local origin point
+	// get global (x,y,z) coordinates of local origin point
 	fp origin_x() const { return global_coords_.origin_x(); }
 	fp origin_y() const { return global_coords_.origin_y(); }
 	fp origin_z() const { return global_coords_.origin_z(); }
 
+	// find patch containing (ray from origin to) given local (x,y,z)
+	// ... if there are multiple patches containing the position,
+	//     we return the one which would still contain it if patches
+	//     didn't overlap; if multiple patches satisfy this criterion
+	//     then it's arbitrary which one we return
+	// ... if no patch contains the position (for a non--full-sphere
+	//     patch system), or the position is at the origin, then
+	//     we do an  error_exit()  and don't return to the caller
+	const patch& patch_containing_local_xyz(fp x, fp y, fp z)
+		const;
+
  
 	//
 	// ***** meta-info about the entire patch system *****
@@ -104,11 +115,15 @@ public:
 	int N_patches() const { return N_patches_; }
 
 	// get patches by patch number
-	const patch &ith_patch(int pn) const
+	const patch& ith_patch(int pn) const
 		{ return * all_patches_[pn]; }
-	patch &ith_patch(int pn)
+	patch& ith_patch(int pn)
 		{ return * all_patches_[pn]; }
 
+	// find a patch by +/- xyz "ctype"
+	const patch& plus_or_minus_xyz_patch(bool is_plus, char ctype)
+		const;
+
 	// find a patch by name, return patch number; error_exit() if not found
 	int patch_number_of_name(const char* name) const;