try to work around some fp roundoff problems which were causing

assertion failures:
- fix a bunch of asserts to use my fuzzy-arithmetic classes
- also guard assertions on sign(z) by a test that z is fuzzily nonzero


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

1 files changed, 24 insertions, 21 deletions

diff --git a/src/patch/coords.cc b/src/patch/coords.cc
index a87469d..795929a 100644
--- a/src/patch/coords.cc
+++ b/src/patch/coords.cc
@@ -97,7 +97,8 @@ return jtutil::modulo_reduce(dang, 360.0, min_dang, max_dang);
 
 //**************************************
 
-// not valid in xy plane (cos(mu) == 0 || cos(nu) == 0)
+// valid in +/- z patch
+// not valid in xy plane (z == 0, i.e. (cos(mu) == 0 || cos(nu) == 0))
 // c.f. page 9.2 of my mpe notes
 namespace local_coords
 	  {
@@ -123,17 +124,18 @@ z = sign_z * r * temp;
 x = x_over_z * z;
 y = y_over_z * z;
 
-if (r != 0.0)
-   then {
-	assert( signum(x) == sign_x );
-	assert( signum(y) == sign_y );
-	}
+assert( jtutil::fuzzy<fp>::NE(z, 0.0) );
+if (jtutil::fuzzy<fp>::NE(x, 0.0))
+   then assert( signum(x) == sign_x );
+if (jtutil::fuzzy<fp>::NE(y, 0.0))
+   then assert( signum(y) == sign_y );
 }
 	  }
 
 //**************************************
 
-// not valid in xz plane (sin(mu) == 0 || sin(phi) == 0)
+// valid in +/- y patch
+// not valid in xz plane (y == 0, i.e. (sin(mu) == 0 || sin(phi) == 0)
 // c.f. page 9.3 of my mpe notes
 namespace local_coords
 	  {
@@ -162,17 +164,18 @@ y = sign_y * r * temp;
 z = z_over_y * y;
 x = x_over_y * y;
 
-if (r != 0.0)
-   then {
-	assert( signum(z) == sign_z );
-	assert( signum(x) == sign_x );
-	}
+assert( jtutil::fuzzy<fp>::NE(y, 0.0) );
+if (jtutil::fuzzy<fp>::NE(z, 0.0))
+   then assert( signum(z) == sign_z );
+if (jtutil::fuzzy<fp>::NE(x, 0.0))
+   then assert( signum(x) == sign_x );
 }
 	  }
 
 //**************************************
 
-// not valid in yz plane (sin(nu) == 0 || cos(phi) == 0)
+// valid in +/- x patch
+// not valid in yz plane (x == 0, i.e. (sin(nu) == 0 || cos(phi) == 0))
 // c.f. page 9.4 of my mpe notes
 namespace local_coords
 	  {
@@ -200,11 +203,11 @@ x = sign_x * r * temp;
 z = z_over_x * x;
 y = y_over_x * x;
 
-if (r != 0.0)
-   then {
-	assert( signum(z) == sign_z );
-	assert( signum(y) == sign_y );
-	}
+assert( jtutil::fuzzy<fp>::NE(x, 0.0) );
+if (jtutil::fuzzy<fp>::NE(z, 0.0))
+   then assert( signum(z) == sign_z );
+if (jtutil::fuzzy<fp>::NE(y, 0.0))
+   then assert( signum(y) == sign_y );
 }
 	  }
 
@@ -294,7 +297,7 @@ const fp tan2_nu = pow2(tan_nu);
 fp x, y, z;
 xyz_of_r_mu_nu(r,mu,nu, x,y,z);
 
-assert(r != 0.0);
+assert( jtutil::fuzzy<fp>::NE(r, 0.0) );
 const fp rinv = 1.0/r;
 partial_x_wrt_r = x*rinv;
 partial_y_wrt_r = y*rinv;
@@ -336,7 +339,7 @@ const fp tan2_phi_bar = pow2(tan_phi_bar);
 fp x, y, z;
 xyz_of_r_mu_phi(r,mu,phi, x,y,z);
 
-assert(r != 0.0);
+assert( jtutil::fuzzy<fp>::NE(r, 0.0) );
 const fp rinv = 1.0/r;
 partial_x_wrt_r = x*rinv;
 partial_y_wrt_r = y*rinv;
@@ -377,7 +380,7 @@ const fp tan2_phi    = pow2(tan_phi);
 fp x, y, z;
 xyz_of_r_nu_phi(r,nu,phi, x,y,z);
 
-assert(r != 0.0);
+assert( jtutil::fuzzy<fp>::NE(r, 0.0) );
 const fp rinv = 1.0/r;
 partial_x_wrt_r = x*rinv;
 partial_y_wrt_r = y*rinv;