initial import into this cvs repository

-- source is JT personal cvs repository ~/mpe/util/coords.{cc,hh}


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

1 files changed, 946 insertions, 0 deletions

diff --git a/src/patch/coords.cc b/src/patch/coords.cc
new file mode 100644
index 0000000..578ab3d
--- /dev/null
+++ b/src/patch/coords.cc
@@ -0,0 +1,946 @@
+// coords.cc - classes for mpe coordinates
+// $Id$
+
+//
+// wr_coord::wr_coord
+// wr_coord::pars::operator==
+// wr_coord::pars::printme
+//
+// wr_coord::r_of_wr
+//* r_of_wr_callback
+// wr_coord::wr_of_r
+//* wr_of_r_with_args - like wr_coord::wr_of_r() but with explicit args r0,a,b,c
+// wr_coord::dwr_dr
+// wr_coord::d2wr_dr
+//
+// coords::{mu,nu,phi}_of_{{mu,nu,phi}}
+// coords::partial_{mu,nu,phi}_wrt_{mu,nu,phi}_at_const_{mu,nu,phi}
+// coords::partial2_{mu,nu,phi}_wrt_{{mu,nu,phi}}_at_const_{{mu,nu,phi}}
+//
+// coords::xyz_of_{{mu,nu,phi}}
+// coords::xyz_of_r_{{mu,nu,phi}}
+// coords::{{mu,nu,phi}}_of_xyz
+// coords::r_{{mu,nu,phi}}_of_xyz
+// coords::theta_phi_of_{{mu,nu,phi}}
+//
+// coords::mu_nu_phi::name_of_index
+// coords::mu_nu_phi::name_of_coords_set
+// coords::r_rho_sigma::name_of_index
+//
+
+#include <math.h>
+#include <float.h>
+#include <assert.h>
+#include <limits.h>
+
+#include <jt/stdc.h>
+#include <jt/util.h>
+#include <jt/util++.hh>
+#include <jt/fortran.h>
+
+#include "policy.hh"
+#include "fp.hh"
+#include "misc.hh"
+#include "coords.hh"
+
+#ifdef FP_IS_FLOAT
+  #include <jt/sfmm.h>
+#endif
+#ifdef FP_IS_DOUBLE
+  #include <jt/dfmm.h>
+#endif
+
+using jtutil::arctan_xy;
+using jtutil::pow2;
+using jtutil::hypot3;
+
+// prototypes for functions local to this file
+namespace {
+fp mixed_210_wr_callback(const fp *pr);
+fp wr_of_r_with_args(fp r0, fp a, fp b, fp c, fp r);
+	  };
+
+//*****************************************************************************
+//*****************************************************************************
+// wr_coord - wr(r) nonuniform gridding coordinate
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// This function constructs a  wr_coord  object.
+//
+wr_coord::wr_coord(const pars &pars_in)
+	: r0_(pars_in.r0),
+	  a_(pars_in.a), b_(pars_in.b), c_(pars_in.c)
+{
+}
+
+//*****************************************************************************
+
+//
+// This function tests whether or not two  struct pars  objects
+// are element-for-element equal.  It returns true for equal, false for
+// not-equal.
+//
+bool wr_coord::pars::operator==(const struct pars &rhs)
+{
+return     (r0 == rhs.r0)
+	&& (a  == rhs.a)
+	&& (b  == rhs.b)
+	&& (c  == rhs.c);
+}
+
+//*************************************
+
+//
+// This function tests whether or not a  struct pars  object is
+// element-for-element equal to the corresponding "stuff" (members or
+// access functions) of a  wr_coord  object.  It returns true for equal,
+// false for not-equal.
+//
+bool wr_coord::pars::operator==(const wr_coord &rhs)
+{
+return     (r0 == rhs.r0())
+	&& (a  == rhs.a())
+	&& (b  == rhs.b())
+	&& (c  == rhs.c());
+}
+
+//*****************************************************************************
+
+//
+// This function prints a  struct_pars  object in parameter-file format.
+//
+void wr_coord::pars::printme() const
+{
+printf("\n");
+printf("wr_coord.r0\t%g\n", r0);
+printf("wr_coord.a\t%g\n", a);
+printf("wr_coord.b\t%g\n", b);
+printf("wr_coord.c\t%g\n", c);
+}
+
+//*****************************************************************************
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// This function computes the mixed-210 nonlinear grid warping function
+//  r = r(wr)  by numerically inverting the  wr = wr(r)  function.  We
+// use Brent's  zeroin()  function, as given in
+//	G. E. Forsythe, M. A. Malcolm, and C. B. Moler
+//	"Computer Methods for Mathematical Computations"
+//	Prentice-Hall, Englewood Cliffs, 1977
+// with a bounding interval computed by successive doubling or halving
+// in r, starting from the inner grid boundary r0.
+//
+// Arguments:
+// {r0,a,b,c} = (in) The mixed-210 parameters.
+// wr_star = (in) The mixed-210 wr coordinate.
+//
+// Result:
+// The function returns the r coordinate.
+//
+// Bugs:
+// - This function is somewhat inefficient.  However, in practice it's
+//   only called during initialization to set up the coefficient arrays,
+//   so the inefficiency doesn't make much difference.
+// - We implicitly assume that the Fortran  zeroin()  function uses <fp>
+//   precision; disaster is likely if this isn't true.
+//
+
+// static variables for this function
+static fp static_r0, static_a, static_b, static_c;
+static fp static_wr_star;
+
+// the function itself
+fp wr_coord::r_of_wr(fp wr_star) const
+{
+// compute bounding interval by successive doubling or halving in r
+// we assume that r=r0 corresponds to wr=0
+fp r_lo, r_hi;
+if (wr_star >= 0.0)
+   then {
+	// search outwards from r0 by successive doubling in r
+	fp r;
+		for (r = 2.0 * r0()  ;  wr_of_r(r) <= wr_star  ;  r *= 2.0)
+		{
+		}
+	r_lo = 0.5 * r;
+	r_hi = r;
+	}
+   else {
+	// search inwards from r0 by successive *halving* in r
+	fp r;
+		for (r = 0.5 * r0()  ;  wr_of_r(r) >= wr_star  ;  r *= 0.5)
+		{
+		}
+	r_lo = r;
+	r_hi = 2.0 * r;
+	}
+
+// firewall check
+if (! ((wr_of_r(r_lo) <= wr_star) && (wr_of_r(r_hi) >= wr_star)) )
+   then error_exit(PANIC_EXIT,
+"***** wr_coord::r_of_wr:\n"
+"        \"bounding interval\" [r_lo,r_hi] doesn't bound a solution!\n"
+"        (this should never happen!)\n"
+"        r0=%g   a=%g   b=%g   c=%g\n"
+"        wr_star=%g\n"
+"        r_lo=%g   wr_of_r(r0,a,b,c, r_lo)=%g\n"
+"        r_hi=%g   wr_of_r(r0,a,b,c, r_hi)=%g\n"
+,
+		   double(r0()), double(a()), double(b()), double(c()),
+		   double(wr_star),
+		   double(r_lo), double(wr_of_r(r_lo)),
+		   double(r_hi), double(wr_of_r(r_hi)));	/*NOTREACHED*/
+
+// solve  wr(r) - wr_star = 0  for r
+// using static variables to pass the parameters to the callback fn
+static_r0 = r0();
+static_a = a();
+static_b = b();
+static_c = c();
+static_wr_star = wr_star;
+fp tol = 0.0;			// solve eqn as accurately as possible
+return fmm::FORTRAN_NAME(zeroin)(&r_lo, &r_hi, &mixed_210_wr_callback, &tol);
+}
+
+//*****************************************************************************
+
+//
+// This is a callback function called by  mixed_210_r_of_wr()  (above)
+// via Brent's  zeroin()  function.  This function uses a C/Fortran-compatible
+// interface.
+//
+// Arguments (implicit, as global variables):
+// static_{r0,a,b,c} = (in) The {r0,a,b,c} mixed-210 parameters.
+// static_wr_star = (in) The wr_star value from  mixed_210_r_of_wr() .
+//
+// Arguments (explicit):
+// pr --> (in) The r coordinate.
+//
+// Results:
+// This function returns  wr(r) - wr_star .
+//
+namespace {
+fp mixed_210_wr_callback(const fp *pr)
+{
+fp wr = wr_of_r_with_args(static_r0, static_a, static_b, static_c, *pr);
+return wr - static_wr_star;
+}
+	  };
+
+//*****************************************************************************
+
+//
+// This function computes the mixed-210 grid warping function  wr(r) ,
+// as per p.363 of my sssf notes.
+//
+// Arguments:
+// r = (in) The r coordinate.
+//
+// Result:
+// The function returns the mixed-210  wr  coordinate.
+//
+fp wr_coord::wr_of_r(fp r) const
+{
+return wr_of_r_with_args(r0(), a(), b(), c(), r);
+}
+
+//*****************************************************************************
+
+//
+// This function is the same as  wr_coord::wr_of_r()  (above), except
+// that it's a non-member function with explicit arguments r0, a, b, and c.
+//
+// Arguments:
+// {r0,a,b,c} = (in) The mixed-210 parameters.
+// r = (in) The r coordinate.
+//
+// Result:
+// The function returns the mixed-210  wr  coordinate.
+//
+namespace {
+fp wr_of_r_with_args(fp r0, fp a, fp b, fp c, fp r)
+{
+fp x = r/r0;
+fp a_term = (a == 0.0)  ?  0.0  :  (r0/a) * (1.0 - 1.0/x);
+fp b_term = (b == 0.0)  ?  0.0  :  (r0/b) * log(x);
+fp c_term = (c == 0.0)  ?  0.0  :  (r0/c) * (x - 1.0);
+
+return a_term + b_term + c_term;
+}
+	  };
+
+//*****************************************************************************
+
+//
+// This function computes the derivative of the mixed-210 grid warping
+// function  wr(r) , as per p.363 of my sssf notes.
+//
+// Arguments:
+// r = (in) The r coordinate.
+//
+// Result:
+// The function returns the derivative $d wr / d r$.
+//
+fp wr_coord::dwr_dr(fp r) const
+{
+fp x = r/r0();
+fp a_term = (a() == 0.0)  ?  0.0  :  1.0 / (a() * x*x);
+fp b_term = (b() == 0.0)  ?  0.0  :  1.0 / (b() * x);
+fp c_term = (c() == 0.0)  ?  0.0  :  1.0 / c();
+
+return a_term + b_term + c_term;
+}
+
+//*****************************************************************************
+
+//
+// This function computes the second derivative of the mixed-210 grid warping
+// function  wr(r) , $d^2 wr / d r^2$.
+//
+// Arguments:
+// {r0,a,b,c} = (in) The mixed-210 parameters.
+// r = (in) The mixed-210 r coordinate.
+//
+// Result:
+// The function returns the second derivative.
+//
+fp wr_coord::d2wr_dr2(fp r) const
+{
+fp x = r/r0();
+fp a_term = (a() == 0.0)  ?  0.0  :  -2.0 / (a() * x*x*x * r0());
+fp b_term = (b() == 0.0)  ?  0.0  :  -1.0 / (b() * x*x * r0());
+
+return( a_term + b_term );
+}
+
+//*****************************************************************************
+//*****************************************************************************
+// coords - angular coordinate conversions and other misc coordinates stuff
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// these functions give any one of {mu,nu,phi} in terms of the other two
+//
+
+// ill-conditioned near z axis
+// not valid in xy plane
+fp coords::phi_of_mu_nu(fp mu, fp nu)
+{
+assert(! fuzzy<fp>::EQ(fabs(mu), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(nu), 0.5*PI));
+
+fp y_over_z = tan(mu);
+fp x_over_z = tan(nu);
+return arctan_xy(x_over_z, y_over_z);
+}
+
+// ill-conditioned near y axis
+// not valid in xz plane
+fp coords::nu_of_mu_phi(fp mu, fp phi)
+{
+fp  mu_bar = 0.5*PI - mu ;
+fp phi_bar = 0.5*PI - phi;
+assert(! fuzzy<fp>::EQ(fabs( mu_bar), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(phi_bar), 0.5*PI));
+
+fp z_over_y = tan( mu_bar);
+fp x_over_y = tan(phi_bar);
+return arctan_xy(z_over_y, x_over_y);
+}
+
+// ill-conditioned near x axis
+// not valid in yz plane
+fp coords::mu_of_nu_phi(fp nu, fp phi)
+{
+fp nu_bar = 0.5*PI - nu ;
+assert(! fuzzy<fp>::EQ(fabs(nu_bar), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(phi   ), 0.5*PI));
+
+fp z_over_x = tan(nu_bar);
+fp y_over_x = tan(phi   );
+return arctan_xy(z_over_x, y_over_x);
+}
+
+//*****************************************************************************
+
+//
+// these functions compute the partial derivative of any of (mu,nu,phi)
+// wrt any other with the 3rd held constant, as per page 13 of my mpe notes
+//
+
+fp coords::partial_phi_wrt_mu_at_const_nu(fp mu, fp nu)
+{
+fp x, y, z;
+coords::xyz_of_mu_nu(mu, nu, x, y, z);
+return (x/z) * (y*y + z*z)/(x*x + y*y);
+}
+
+fp coords::partial_phi_wrt_nu_at_const_mu(fp mu, fp nu)
+{
+fp x, y, z;
+coords::xyz_of_mu_nu(mu, nu, x, y, z);
+return (-y/z) * (x*x + z*z)/(x*x + y*y);
+}
+
+//****************************************
+
+fp coords::partial_mu_wrt_nu_at_const_phi(fp nu, fp phi)
+{
+fp x, y, z;
+coords::xyz_of_nu_phi(nu, phi, x, y, z);
+return (y/x) * (x*x + z*z) / (y*y + z*z);
+}
+
+fp coords::partial_mu_wrt_phi_at_const_nu(fp nu, fp phi)
+{
+fp x, y, z;
+coords::xyz_of_nu_phi(nu, phi, x, y, z);
+return (z/x) * (x*x + y*y) / (y*y + z*z);
+}
+
+//****************************************
+
+fp coords::partial_nu_wrt_mu_at_const_phi(fp mu, fp phi)
+{
+fp x, y, z;
+coords::xyz_of_mu_phi(mu, phi, x, y, z);
+return (x/y) * (y*y + z*z) / (x*x + z*z);
+}
+
+fp coords::partial_nu_wrt_phi_at_const_mu(fp mu, fp phi)
+{
+fp x, y, z;
+coords::xyz_of_mu_phi(mu, phi, x, y, z);
+return (-z/y) * (x*x + y*y) / (x*x + z*z);
+}
+
+//*****************************************************************************
+
+//
+// These functions compute the 2nd partial derivative of any of (mu,nu,phi)
+// wrt the other two, with the non-differentiated coordinates held constant.
+//
+// The formulas are from the Maple functions of the same name in
+//  ../maple/coords.maple , which "just" differentiate the individual
+//  {mu,nu,phi}_of_{{mu,nu,phi}}()  functions.  The Maple output from
+//  print_coord_partial2_exprs()  is include here (in comments).
+//
+
+//                                                   2    2      2    2
+//                                               y (z  + y ) x (x  - z )
+//         partial^2 phi / partial mu^2|nu = , 2 -----------------------
+//                                                     2   2    2 2
+//                                                    z  (x  + y )
+fp coords::partial2_phi_wrt_mu2_at_const_nu(fp mu, fp nu)
+{
+fp x, y, z;
+xyz_of_mu_nu(mu, nu, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return 2.0 * x*y * (xx-zz) * (yy+zz) / (zz * pow2(xx+yy));
+}
+
+
+//                                                  2    2    2    2    2    2
+//                                                (z  + y ) (z  + x ) (x  - y )
+//  partial^2 phi / partial (mu,nu)|(nu,mu) = , - -----------------------------
+//                                                         2   2    2 2
+//                                                        z  (x  + y )
+fp coords::partial2_phi_wrt_mu_nu_at_const_nu_mu(fp mu, fp nu)
+{
+fp x, y, z;
+xyz_of_mu_nu(mu, nu, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return - (xx-yy) * (xx+zz) * (yy+zz) / (zz * pow2(xx+yy));
+}
+
+//                                                   2    2       2    2
+//                                               y (z  + x ) x (-z  + y )
+//        partial^2 phi / partial nu^2|mu = , -2 ------------------------
+//                                                     2   2    2 2
+//                                                    z  (x  + y )
+fp coords::partial2_phi_wrt_nu2_at_const_mu(fp mu, fp nu)
+{
+fp x, y, z;
+xyz_of_mu_nu(mu, nu, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return -2.0 * x*y * (xx+zz) * (yy-zz) / (zz * pow2(xx+yy));
+}
+
+//****************************************
+
+//                                                     2    2    2    2
+//                                               y z (z  + x ) (x  - y )
+//         partial^2 mu / partial nu^2|phi = , 2 -----------------------
+//                                                     2   2    2 2
+//                                                    x  (z  + y )
+fp coords::partial2_mu_wrt_nu2_at_const_phi(fp nu, fp phi)
+{
+fp x, y, z;
+xyz_of_nu_phi(nu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return 2.0 * y*z * (xx-yy) * (xx+zz) / (xx * pow2(yy+zz));
+}
+
+//                                                  2    2    2    2     2    2
+//                                                (x  + y ) (z  + x ) (-z  + y )
+// partial^2 mu / partial (nu,phi)|(phi,nu) = , - ------------------------------
+//                                                         2   2    2 2
+//                                                        x  (z  + y )
+fp coords::partial2_mu_wrt_nu_phi_at_const_phi_nu(fp nu, fp phi)
+{
+fp x, y, z;
+xyz_of_nu_phi(nu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return - (xx+yy) * (xx+zz) * (yy-zz) / (xx * pow2(yy+zz));
+}
+
+//                                                    2    2      2    2
+//                                                y (x  + y ) z (x  - z )
+//         partial^2 mu / partial phi^2|nu = , -2 -----------------------
+//                                                      2   2    2 2
+//                                                     x  (z  + y )
+fp coords::partial2_mu_wrt_phi2_at_const_nu(fp nu, fp phi)
+{
+fp x, y, z;
+xyz_of_nu_phi(nu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return -2.0 * y*z * (xx+yy) * (xx-zz) / (xx * pow2(yy+zz));
+}
+
+//****************************************
+
+//                                                    2    2      2    2
+//                                                x (z  + y ) z (x  - y )
+//         partial^2 nu / partial mu^2|phi = , -2 -----------------------
+//                                                      2   2    2 2
+//                                                     y  (z  + x )
+
+fp coords::partial2_nu_wrt_mu2_at_const_phi(fp mu, fp phi)
+{
+fp x, y, z;
+xyz_of_mu_phi(mu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return -2.0 * x*z * (xx-yy) * (yy+zz) / (yy * pow2(xx+zz));
+}
+
+//                                                  2    2    2    2    2    2
+//                                                (x  + y ) (z  + y ) (x  - z )
+//   partial^2 nu / partial (mu,phi)|(phi,mu) = , -----------------------------
+//                                                         2   2    2 2
+//                                                        y  (z  + x )
+fp coords::partial2_nu_wrt_mu_phi_at_const_phi_mu(fp mu, fp phi)
+{
+fp x, y, z;
+xyz_of_mu_phi(mu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return (xx+yy) * (xx-zz) * (yy+zz) / (yy * pow2(xx+zz));
+}
+
+//                                                     2    2     2    2
+//                                               z x (x  + y ) (-z  + y )
+//        partial^2 nu / partial phi^2|mu = , -2 ------------------------
+//                                                     2   2    2 2
+//                                                    y  (z  + x )
+fp coords::partial2_nu_wrt_phi2_at_const_mu(fp mu, fp phi)
+{
+fp x, y, z;
+xyz_of_mu_phi(mu, phi, x, y, z);
+
+fp xx = x*x;
+fp yy = y*y;
+fp zz = z*z;
+
+return -2.0 * x*z * (xx+yy) * (yy-zz) / (yy * pow2(xx+zz));
+}
+
+//*****************************************************************************
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// these functions convert ((mu,nu,phi)) <--> (x,y,z) on the unit sphere
+//
+
+// valid near z axis (mu small, nu small)
+// not valid near xy plane
+void coords::xyz_of_mu_nu(fp mu, fp nu,
+			  fp &x, fp &y, fp &z)
+{
+assert(! fuzzy<fp>::EQ(fabs(mu), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(nu), 0.5*PI));
+
+fp y_over_z = tan(mu);
+fp x_over_z = tan(nu);
+fp temp = 1.0 / sqrt(1.0 + y_over_z*y_over_z + x_over_z*x_over_z);
+
+z = temp;
+x = x_over_z * temp;
+y = y_over_z * temp;
+}
+
+// valid near y axis (pi/2-mu small, pi/2-phi small)
+// not valid near xz plane
+void coords::xyz_of_mu_phi(fp mu, fp phi,
+			   fp &x, fp &y, fp &z)
+{
+fp  mu_bar = 0.5*PI - mu ;
+fp phi_bar = 0.5*PI - phi;
+assert(! fuzzy<fp>::EQ(fabs( mu_bar), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(phi_bar), 0.5*PI));
+
+fp z_over_y = tan( mu_bar);
+fp x_over_y = tan(phi_bar);
+fp temp = 1.0 / sqrt(1.0 + z_over_y*z_over_y + x_over_y*x_over_y);
+
+y = temp;
+z = z_over_y * temp;
+x = x_over_y * temp;
+}
+
+// valid near x axis (pi/2-nu small, phi small)
+// not valid near yz plane
+void coords::xyz_of_nu_phi(fp nu, fp phi,
+			   fp &x, fp &y, fp &z)
+{
+fp nu_bar = 0.5*PI - nu ;
+assert(! fuzzy<fp>::EQ(fabs(nu_bar), 0.5*PI));
+assert(! fuzzy<fp>::EQ(fabs(phi   ), 0.5*PI));
+
+fp z_over_x = tan(nu_bar);
+fp y_over_x = tan(phi   );
+fp temp = 1.0 / sqrt(1.0 + z_over_x*z_over_x + y_over_x*y_over_x);
+
+x = temp;
+z = z_over_x * temp;
+y = y_over_x * temp;
+}
+
+//****************************************
+
+// valid near z axis (mu small, nu small)
+// ill-conditioned near x or y axes
+void coords::mu_nu_of_xyz(fp x, fp y,
+			  fp z, fp &mu, fp &nu)
+{
+assert( fuzzy<fp>::NE(x,0.0) || fuzzy<fp>::NE(y,0.0) || fuzzy<fp>::NE(z,0.0) );
+
+mu = arctan_xy(z,y);
+nu = arctan_xy(z,x);
+}
+
+// valid near y axis (pi/2-mu small, pi/2-phi small)
+// ill-conditioned near x or z axes
+void coords::mu_phi_of_xyz(fp x, fp y, fp z,
+			   fp &mu, fp &phi)
+{
+assert( fuzzy<fp>::NE(x,0.0) || fuzzy<fp>::NE(y,0.0) || fuzzy<fp>::NE(z,0.0) );
+
+mu  = arctan_xy(z,y);
+phi = arctan_xy(x,y);
+}
+
+// valid near x axis (pi/2-nu small, pi/2-phi small)
+// ill-conditioned near y or z axes
+void coords::nu_phi_of_xyz(fp x, fp y, fp z,
+			   fp &nu, fp &phi)
+{
+assert( fuzzy<fp>::NE(x,0.0) || fuzzy<fp>::NE(y,0.0) || fuzzy<fp>::NE(z,0.0) );
+
+nu  = arctan_xy(z,x);
+phi = arctan_xy(x,y);
+}
+
+//*****************************************************************************
+
+//
+// these functions convert (r,(mu,nu,phi)) <--> generic (x,y,z)
+//
+
+// valid near z axis (mu small, nu small)
+// not valid near xy plane
+void coords::xyz_of_r_mu_nu(fp r, fp mu, fp nu,
+			    fp &x, fp &y, fp &z)
+{
+xyz_of_mu_nu(mu, nu, x, y, z);
+x *= r;
+y *= r;
+z *= r;
+}
+
+// valid near y axis (pi/2-mu small, pi/2-phi small)
+// not valid near xz plane
+void coords::xyz_of_r_mu_phi(fp r, fp mu, fp phi, 
+			     fp &x, fp &y, fp &z)
+{
+xyz_of_mu_phi(mu, phi, x, y, z);
+x *= r;
+y *= r;
+z *= r;
+}
+
+// valid near x axis (pi/2-nu small, phi small)
+// not valid near yz plane
+void coords::xyz_of_r_nu_phi(fp r, fp nu, fp phi,
+			     fp &x, fp &y, fp &z)
+{
+xyz_of_nu_phi(nu, phi, x, y, z);
+x *= r;
+y *= r;
+z *= r;
+}
+
+//****************************************
+
+// valid near z axis (mu small, nu small)
+// not valid near x or y axes
+void coords::r_mu_nu_of_xyz(fp x, fp y, fp z,
+			    fp &r, fp &mu, fp &nu)
+{
+r = hypot3(x, y, z);
+
+assert( fuzzy<fp>::NE(r,0.0) );
+mu_nu_of_xyz(x, y, z, mu, nu);
+}
+
+// valid near y axis (pi/2-mu small, pi/2-phi small)
+// not valid near x or z axes
+void coords::r_mu_phi_of_xyz(fp x, fp y, fp z,
+			     fp &r, fp &mu, fp &phi)
+{
+r = hypot3(x, y, z);
+
+assert( fuzzy<fp>::NE(r,0.0) );
+mu_phi_of_xyz(x, y, z, mu, phi);
+}
+
+// valid near x axis (pi/2-nu small, pi/2-phi small)
+// not valid near y or z axes
+void coords::r_nu_phi_of_xyz(fp x, fp y, fp z,
+			     fp &r, fp &nu, fp &phi)
+{
+r = hypot3(x, y, z);
+
+assert( fuzzy<fp>::NE(r,0.0) );
+nu_phi_of_xyz(x, y, z, nu, phi);
+}
+
+//*****************************************************************************
+
+//
+// these functions convert the usual polar spherical (theta,phi)
+// <--> (x,y,z) on the unit sphere
+//
+
+void coords::xyz_of_theta_phi(fp theta, fp phi,
+			      fp &x, fp &y, fp &z)
+{
+z = cos(theta);
+x = sin(theta) * cos(phi);
+y = sin(theta) * sin(phi);
+}
+
+void coords::theta_phi_of_xyz(fp x, fp y, fp z,
+			      fp &theta, fp &phi)
+{
+assert( fuzzy<fp>::NE(x,0.0) || fuzzy<fp>::NE(y,0.0) || fuzzy<fp>::NE(z,0.0) );
+
+theta = arctan_xy( z , hypot(x,y) );
+phi = arctan_xy(x, y);
+}
+
+//*****************************************************************************
+
+//
+// these functions convert the usual polar spherical (r,theta,phi) <--> (x,y,z)
+//
+
+void coords::xyz_of_r_theta_phi(fp r, fp theta, fp phi,
+				fp &x, fp &y, fp &z)
+{
+xyz_of_theta_phi(theta, phi, x, y, z);
+x *= r;
+y *= r;
+z *= r;
+}
+
+void coords::r_theta_phi_of_xyz(fp x, fp y, fp z,
+				fp &r, fp &theta, fp &phi)
+{
+r = hypot3(x, y, z);
+
+assert( fuzzy<fp>::NE(r,0.0) );
+theta_phi_of_xyz(x, y, z, theta, phi);
+}
+
+//*****************************************************************************
+
+//
+// these functions convert ((mu,nu,phi)) <--> usual polar spherical (theta,phi)
+// ... note phi is the same coordinate in both systems
+//
+
+void coords::theta_phi_of_mu_nu(fp mu, fp nu,
+				fp &ps_theta, fp &ps_phi)
+{
+fp tan_mu = tan(mu);
+fp tan_nu = tan(nu);
+ps_theta = atan(hypot(tan_mu, tan_nu));		// page 7 of my mpe notes
+ps_phi = arctan_xy(tan_nu, tan_mu);		// page 6 of my mpe notes
+}
+
+// Bugs: computes ps_phi via trig, even though it's trivially == phi
+void coords::theta_phi_of_mu_phi(fp mu, fp phi,
+				 fp &ps_theta, fp &ps_phi)
+{
+fp x, y, z;
+xyz_of_mu_phi(mu, phi, x, y, z);
+theta_phi_of_xyz(x, y, z, ps_theta, ps_phi);
+}
+
+// Bugs: computes ps_phi via trig, even though it's trivially == phi
+void coords::theta_phi_of_nu_phi(fp nu, fp phi,
+				 fp &ps_theta, fp &ps_phi)
+{
+fp x, y, z;
+xyz_of_nu_phi(nu, phi, x, y, z);
+theta_phi_of_xyz(x, y, z, ps_theta, ps_phi);
+}
+
+//****************************************
+
+void coords::mu_nu_of_theta_phi(fp ps_theta, fp ps_phi,
+				fp &mu, fp &nu)
+{
+fp x, y, z;
+xyz_of_theta_phi(ps_theta, ps_phi, x, y, z);
+mu_nu_of_xyz(x, y, z, mu, nu);
+}
+
+// Bugs: computes phi via trig, even though it's trivially == ps_phi
+void coords::mu_phi_of_theta_phi(fp ps_theta, fp ps_phi,
+				 fp &mu, fp &phi)
+{
+fp x, y, z;
+xyz_of_theta_phi(ps_theta, ps_phi, x, y, z);
+mu_phi_of_xyz(x, y, z, mu, phi);
+}
+
+// Bugs: computes phi via trig, even though it's trivially == ps_phi
+void coords::nu_phi_of_theta_phi(fp ps_theta, fp ps_phi,
+				 fp &nu, fp &phi)
+{
+fp x, y, z;
+xyz_of_theta_phi(ps_theta, ps_phi, x, y, z);
+nu_phi_of_xyz(x, y, z, nu, phi);
+}
+
+//*****************************************************************************
+
+//
+// This function computes a human-readable name from a (mu,nu,phi)
+// coordinate index
+//
+const char *coords::mu_nu_phi::name_of_index(int c)
+{
+static const char *name_table[] = { "mu", "nu", "phi" };
+
+if ((c >= 0) && (c < N_coords))
+   then return name_table[c];
+   else error_exit(PANIC_EXIT,
+"***** coords::mu_nu_phi::name_of_index:\n"
+"        c=%d isn't a valid coordinate index!\n"
+,
+		   c);						/*NOTREACHED*/
+}
+
+//*************************************
+
+//
+// This function computes a human-readable name from a (mu,nu,phi)
+// coordinates set.
+//
+const char *coords::mu_nu_phi::name_of_coords_set(coords_set S)
+{
+//
+// we have to use an if-else chain because the  coords::set_*  constants
+// aren't compile-time constants and hence aren't eligible to be switch
+// case labels
+//
+if	(S == set_empty)
+   then return "{}";
+else if (S == set_mu)
+   then return "mu";
+else if (S == set_nu)
+   then return "nu";
+else if (S == set_phi)
+   then return "phi";
+else if (S == set_mu|set_nu)
+   then return "{mu,nu}";
+else if (S == set_mu|set_phi)
+   then return "{mu,phi}";
+else if (S == set_nu|set_phi)
+   then return "{nu,phi}";
+else if (S == set_mu|set_nu|set_phi)
+   then return "{mu,nu,phi}";
+else	error_exit(PANIC_EXIT,
+"***** coords::mu_nu_phi::name_of_coords_set:\n"
+"        S=%d isn't a valid  coords_set  bit vector!\n"
+,
+		   int(S));					/*NOTREACHED*/
+}
+
+//*****************************************************************************
+
+//
+// This function computes a human-readable name from a (r,rho,sigma)
+// coordinate index
+//
+const char *coords::r_rho_sigma::name_of_index(int c)
+{
+static const char *name_table[] = { "r", "rho", "sigma" };
+
+if ((c >= 0) && (c < N_coords))
+   then return name_table[c];
+   else error_exit(PANIC_EXIT,
+"***** coords::r_rho_sigma::name_of_index:\n"
+"        c=%d isn't a valid coordinate index!\n"
+,
+		   c);						/*NOTREACHED*/
+}