aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/patch/coords.cc946
-rw-r--r--src/patch/coords.hh294
2 files changed, 1240 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*/
+}
diff --git a/src/patch/coords.hh b/src/patch/coords.hh
new file mode 100644
index 0000000..ce455b0
--- /dev/null
+++ b/src/patch/coords.hh
@@ -0,0 +1,294 @@
+// coords.hh -- classes for mpe coordinates
+// $Id$
+//
+// wr_coord - wr(r) nonuniform gridding coordinate
+//
+// coords - misc coordinates-stuff namespace
+//
+
+//
+// prerequisites:
+// <jt/stdc.h> -- for DEGREES_TO_RADIANS and RADIANS_TO_DEGREES
+// fp.hh -- basic floating point stuff
+//
+
+//*****************************************************************************
+
+//
+// This class describes a wr nonuniform gridding coordinate,
+// as per section VII A of
+// @techreport
+// {
+// Thornburg-1999-sssf-evolution,
+// author = "Jonathan Thornburg",
+// title = "A $3+1$ Computational Scheme for Dynamic
+// Spherically Symmetric Black Hole Spacetimes
+// -- {II}: Time Evolution",
+// year = 1999, month = "June",
+// number = "UWThPh-1999-38",
+// eprint = "gr-qc/9906022",
+// note = "To be submitted to Physical Review D",
+// }
+//
+class wr_coord
+ {
+public:
+ fp r_of_wr(fp wr) const;
+ fp wr_of_r(fp r) const;
+
+ fp dwr_dr(fp r) const;
+ fp d2wr_dr2(fp r) const;
+
+ // low-level access to m210 parameters
+ fp r0() const { return r0_; }
+ fp a() const { return a_; }
+ fp b() const { return b_; }
+ fp c() const { return c_; }
+
+ // this structure is used for two purposes
+ // ... reading things in via my par_scan() function
+ // ... bundling constructor arguments together for convenience
+ struct pars
+ {
+ fp r0, a, b, c;
+
+ // element-by-element equality comparison
+ bool operator==(const struct pars &rhs);
+
+ // comparison with corresponding "stuff"
+ // (members or access functions) of a wr_coord object
+ bool operator==(const wr_coord &rhs);
+
+ // print in parameter-file fashion
+ void printme() const;
+ };
+
+ //
+ // This macro generates a comma-separated list of C structure
+ // initializers for par_scan() parameter-table entries describing
+ // members of the above structure:
+ //
+ // Sample usage:
+ // WR_COORD__PARS__MBRS2("wr_coords_pars", wr_coord_pars_str)
+ //
+ #define WR_COORD__PARS__MBRS2(name_, struct_) \
+ PAR_DSCR_VAR2(name_ ".r0", \
+ struct_.r0, \
+ PAR_NO_FLAGS, FP_SCANF_FORMAT), \
+ PAR_DSCR_VAR2(name_ ".a", \
+ struct_.a, \
+ PAR_NO_FLAGS, FP_SCANF_FORMAT), \
+ PAR_DSCR_VAR2(name_ ".b", \
+ struct_.b, \
+ PAR_NO_FLAGS, FP_SCANF_FORMAT), \
+ PAR_DSCR_VAR2(name_ ".c", \
+ struct_.c, \
+ PAR_NO_FLAGS, FP_SCANF_FORMAT) // no comma
+ // end macro
+
+ // constructor, destructor
+ wr_coord(const pars &pars_in);
+ ~wr_coord() { }
+
+private:
+ // we forbid copying and passing by value
+ // by declaring the copy constructor and assignment operator
+ // private, but never defining them
+ wr_coord(const wr_coord &rhs);
+ wr_coord& operator=(const wr_coord &rhs);
+
+private:
+ // if any of a, b, and/or c are 0.0,
+ // we treat this as an infinite value,
+ // i.e. we skip this term in the equations
+ fp r0_, a_, b_, c_;
+ };
+
+//*****************************************************************************
+
+//
+// This namespace contains coordinate conversion functions, focusing
+// on the angular coordinates. It also contains coordinate bit masks
+// and generic coordinate mnemonics.
+//
+// The angular coordinates are
+// mu = rotation about x axis = arctan(y/z)
+// nu = rotation about y axis = arctan(x/z)
+// phi = rotation about z axis = arctan(y/x)
+//
+// Terminology:
+// (a,b,c) == any one of (a,b,c) == a | b | c
+// ((a,b,c)) == any two of (a,b,c) == (a,b) | (a,c) | (b,c)
+// ((r,a,b,c)) == r and any two of (a,b,c) == (r,a,b) | (r,a,c) | (r,b,c)
+// {{a,b,c}} == any two of (a,b,c) separated by an underscore
+// == a_b | a_c | b_c
+// {{r,a,b,c}} == r, an underscore, then any two of (a,b,c) separated
+// by an underscore
+// == r_a_b | r_a_c | r_b_c
+//
+namespace coords
+ {
+ //
+ // coords:: is a *namespace*, not a class, so we need
+ // explicit inline declarations on inline functions to
+ // avoid (conflicting) duplicate defs in separate compilation
+ // units.
+ //
+
+ //
+ // FIXME:
+ // should the reference arguments in the following functions
+ // have restrict qualifiers?
+ //
+ // Bugs:
+ // The functions aren't optimally efficient -- they have
+ // lots of could-be-optimized-away divisions and trig calls.
+ // But in practice this doesn't matter, as we don't call
+ // these functions inside inner loops. (For example, we
+ // cache all the partial derivatives needed for interpatch
+ // tensor transformations, so these functions are only
+ // called to initialize the cache when setting up the
+ // patch system.)
+ //
+
+ // ((mu,nu,phi)) --> the 3rd
+ fp phi_of_mu_nu(fp mu, fp nu);
+ fp nu_of_mu_phi(fp mu, fp phi);
+ fp mu_of_nu_phi(fp nu, fp phi);
+
+ // partial derivatives of any of (mu,nu,phi)
+ // wrt any other with the 3rd held constant
+ // ... partial phi / partial (mu,nu)
+ fp partial_phi_wrt_mu_at_const_nu(fp mu, fp nu );
+ fp partial_phi_wrt_nu_at_const_mu(fp mu, fp nu );
+ // ... partial mu / partial (nu,phi)
+ fp partial_mu_wrt_nu_at_const_phi(fp nu, fp phi);
+ fp partial_mu_wrt_phi_at_const_nu(fp nu, fp phi);
+ // ... partial nu / partial (mu,phi)
+ fp partial_nu_wrt_mu_at_const_phi(fp mu, fp phi);
+ fp partial_nu_wrt_phi_at_const_mu(fp mu, fp phi);
+
+ // 2nd partial derivatives of (mu,nu,phi) among themselves
+ // ... partial^2 phi / partial (mu,nu)^2
+ fp partial2_phi_wrt_mu2_at_const_nu (fp mu, fp nu);
+ fp partial2_phi_wrt_mu_nu_at_const_nu_mu (fp mu, fp nu);
+ fp partial2_phi_wrt_nu2_at_const_mu (fp mu, fp nu);
+ // ... partial^2 mu / partial (nu,phi)^2
+ fp partial2_mu_wrt_nu2_at_const_phi (fp nu, fp phi);
+ fp partial2_mu_wrt_nu_phi_at_const_phi_nu(fp nu, fp phi);
+ fp partial2_mu_wrt_phi2_at_const_nu (fp nu, fp phi);
+ // ... partial^2 nu / partial (mu,phi)^2
+ fp partial2_nu_wrt_mu2_at_const_phi (fp mu, fp phi);
+ fp partial2_nu_wrt_mu_phi_at_const_phi_mu(fp mu, fp phi);
+ fp partial2_nu_wrt_phi2_at_const_mu (fp mu, fp phi);
+
+
+ // (r,(mu,nu,phi)) <--> (x,y,z)
+ // ... on unit sphere
+ void xyz_of_mu_nu(fp mu, fp nu, fp &x, fp &y, fp &z);
+ void xyz_of_mu_phi(fp mu, fp phi, fp &x, fp &y, fp &z);
+ void xyz_of_nu_phi(fp nu, fp phi, fp &x, fp &y, fp &z);
+ void mu_nu_of_xyz(fp x, fp y, fp z, fp &mu, fp &nu);
+ void mu_phi_of_xyz(fp x, fp y, fp z, fp &mu, fp &phi);
+ void nu_phi_of_xyz(fp x, fp y, fp z, fp &nu, fp &phi);
+ // ... generic r
+ void xyz_of_r_mu_nu(fp r, fp mu, fp nu, fp &x, fp &y, fp &z);
+ void xyz_of_r_mu_phi(fp r, fp mu, fp phi, fp &x, fp &y, fp &z);
+ void xyz_of_r_nu_phi(fp r, fp nu, fp phi, fp &x, fp &y, fp &z);
+ void r_mu_nu_of_xyz(fp x, fp y, fp z, fp &r, fp &mu, fp &nu);
+ void r_mu_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &mu, fp &phi);
+ void r_nu_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &nu, fp &phi);
+
+ // usual polar spherical (r,theta,phi) <--> (x,y,z)
+ // ... on unit sphere
+ void xyz_of_theta_phi(fp theta, fp phi, fp &x, fp &y, fp &z);
+ void theta_phi_of_xyz(fp x, fp y, fp z, fp &theta, fp &phi);
+ // ... generic r
+ void xyz_of_r_theta_phi(fp r, fp theta, fp phi, fp &x, fp &y, fp &z);
+ void r_theta_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &theta, fp &phi);
+
+ // ((mu,nu,phi)) <--> usual polar spherical (theta,phi)
+ // ... note phi is the same coordinate in both systems
+ void theta_phi_of_mu_nu(fp mu, fp nu, fp &ps_theta, fp &ps_phi);
+ void theta_phi_of_mu_phi(fp mu, fp phi, fp &ps_theta, fp &ps_phi);
+ void theta_phi_of_nu_phi(fp nu, fp phi, fp &ps_theta, fp &ps_phi);
+ void mu_nu_of_theta_phi(fp ps_theta, fp ps_phi, fp &mu, fp &nu);
+ void mu_phi_of_theta_phi(fp ps_theta, fp ps_phi, fp &mu, fp &phi);
+ void nu_phi_of_theta_phi(fp ps_theta, fp ps_phi, fp &nu, fp &phi);
+
+ //*************************************
+
+ //
+ // We want to manipulate tensor indices, and also do calculations
+ // like "which coordinate do these two patches have in common".
+ // We implement the former in the obvious way (a tensor index is
+ // just an integer), and the latter by Boolean operations on
+ // integers using the coords_set bit vectors below:
+ //
+
+ //
+ // (mu,nu,phi) coordinates
+ // ... in a non-namespace context we use the prefix mnp_
+ // to identify things relating to these coordinates
+ //
+ namespace mu_nu_phi
+ {
+ static const int N_coords = 3;
+
+ // tensor indices
+ static const int index_mu = 0;
+ static const int index_nu = 1;
+ static const int index_phi = 2;
+
+ // sets of coordinates
+ typedef int coords_set;
+
+ // singleton coordinate set for a given index
+ inline coords_set coord_set_of_index(int i) { return 0x1 << i; }
+ static const coords_set set_mu = coord_set_of_index(index_mu );
+ static const coords_set set_nu = coord_set_of_index(index_nu );
+ static const coords_set set_phi = coord_set_of_index(index_phi);
+
+ static const coords_set set_empty = 0;
+ static const coords_set set_all = set_mu | set_nu | set_phi;
+
+ // human-readable coordinate names for debugging etc
+ const char *name_of_index(int c);
+ const char *name_of_coords_set(coords_set S);
+
+ // set complement of coordinates
+ inline coords_set coords_set_not(coords_set S)
+ { return set_all - S; }
+ };
+
+ //
+ // (r,rho,sigma) coordinates
+ // ... in a non-namespace context we use the prefix rrs_
+ // to identify things relating to these coordinates
+ //
+ namespace r_rho_sigma
+ {
+ static const int N_coords = 3;
+ static const int N_angular_coords = 2;
+
+ // tensor indices
+ static const int index_r = 0;
+ static const int index_rho = 1;
+ static const int index_sigma = 2;
+
+ // human-readable coordinate names for debugging etc
+ const char *name_of_index(int c);
+ };
+
+ //*************************************
+
+ //
+ // degrees <--> radians conversions
+ //
+ // ... note coords:: is a *namespace*, not a class,
+ // so we need an explicit inline declaration to avoid
+ // conflicting duplicate defs in separate compilation units
+ inline fp dang_of_ang(fp ang) { return RADIANS_TO_DEGREES * ang; }
+ inline fp ang_of_dang(fp dang) { return DEGREES_TO_RADIANS * dang; }
+
+ };
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-13 14:54:03 +0000