diff options
Diffstat (limited to 'src/patch/coords.cc')
-rw-r--r-- | src/patch/coords.cc | 946 |
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*/ +} |