path: root/src/patch/coords.cc

// coords.cc - coordinate conversions etc
// $Id$

//
// local_coords::{mu,nu,phi}_of_((mu,nu,phi))
// local_coords::xyz_of_((mu,nu,phi))
// local_coords::xyz_of_r_((mu,nu,phi))
// local_coords::((mu,nu,phi))_of_xyz
// local_coords::r_((mu,nu,phi))_of_xyz
// local_coords::theta_phi_of_((mu,nu,phi))
//
// local_coords::mu_nu_phi::name_of_coords_set
//

#include <math.h>
#include <float.h>
#include <assert.h>
#include <limits.h>

#include "jt/stdc.h"
#include "jt/util.hh"

#include "fp.hh"
#include "coords.hh"

using jtutil::error_exit;
using jtutil::arctan_xy;
using jtutil::pow2;
using jtutil::hypot3;

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// 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
namespace local_coords
	  {
fp 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));

const fp y_over_z = tan(mu);
const 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
namespace local_coords
	  {
fp nu_of_mu_phi(fp mu, fp phi)
{
const fp  mu_bar = 0.5*PI - mu ;
const 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));

const fp z_over_y = tan( mu_bar);
const 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
namespace local_coords
	  {
fp mu_of_nu_phi(fp nu, fp phi)
{
const 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));

const fp z_over_x = tan(nu_bar);
const fp y_over_x = tan(phi   );
return arctan_xy(z_over_x, y_over_x);
}
	  }

//******************************************************************************

//
// these functions convert ((mu,nu,phi)) <--> (x,y,z)
//

//**************************************

// valid near z axis (mu small, nu small)
// not valid in xy plane
namespace local_coords
	  {
void xyz_of_r_mu_nu(fp r, 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));

const fp y_over_z = tan(mu);
const fp x_over_z = tan(nu);
const fp temp = 1.0 / sqrt(1.0 + pow2(y_over_z) + pow2(x_over_z));

z = r * temp;
x = r * x_over_z * temp;
y = r * y_over_z * temp;
}
	  }

//**************************************

// valid near y axis (pi/2-mu small, pi/2-phi small)
// not valid in xz plane
namespace local_coords
	  {
void xyz_of_r_mu_phi(fp r, fp mu, fp phi,   fp& x, fp& y, fp& z)
{
const fp  mu_bar = 0.5*PI - mu ;
const 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));

const fp z_over_y = tan( mu_bar);
const fp x_over_y = tan(phi_bar);
const fp temp = 1.0 / sqrt(1.0 + pow2(z_over_y) + pow2(x_over_y));

y = r * temp;
z = r * z_over_y * temp;
x = r * x_over_y * temp;
}
	  }

//**************************************

// valid near x axis (pi/2-nu small, phi small)
// not valid in yz plane
namespace local_coords
	  {
void xyz_of_r_nu_phi(fp r, fp nu, fp phi,   fp& x, fp& y, fp& z)
{
const 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));

const fp z_over_x = tan(nu_bar);
const fp y_over_x = tan(phi   );
const fp temp = 1.0 / sqrt(1.0 + pow2(z_over_x) + pow2(y_over_x));

x = r * temp;
z = r * z_over_x * temp;
y = r * y_over_x * temp;
}
	  }

//******************************************************************************

// valid near z axis (mu small, nu small)
// ill-conditioned near x or y axes
namespace local_coords
	  {
void r_mu_nu_of_xyz(fp x, fp y, fp z,   fp& r, fp& mu, fp& nu)
{
r = hypot3(x, y, z);
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
namespace local_coords
	  {
void mu_phi_of_xyz(fp x, fp y, fp z,   fp& r, fp& mu, fp& phi)
{
r = hypot3(x, y, z);
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
namespace local_coords
	  {
void nu_phi_of_xyz(fp x, fp y, fp z,   fp& r, fp& nu, fp& phi)
{
r = hypot3(x, y, z);
nu  = arctan_xy(z,x);
phi = arctan_xy(x,y);
}
	  }

//******************************************************************************

//
// these functions convert the usual polar spherical (theta,phi) <--> (x,y,z)
//

//**************************************

namespace local_coords
	  {
void xyz_of_r_theta_phi(fp r, fp theta, fp phi,   fp& x, fp& y, fp& z)
{
z = r * cos(theta);
x = r * sin(theta) * cos(phi);
y = r * sin(theta) * sin(phi);
}
	  }

//**************************************

namespace local_coords
	  {
void r_theta_phi_of_xyz(fp x, fp y, fp z,   fp& r, fp& theta, fp& phi)
{
r = hypot3(x, y, z);
theta = arctan_xy( z , hypot(x,y) );
phi = arctan_xy(x, y);
}
	  }

//******************************************************************************

//
// these functions convert ((mu,nu,phi)) <--> usual polar spherical (theta,phi)
// ... note phi is the same coordinate in both systems
//

namespace local_coords
	  {
void theta_phi_of_mu_nu(fp mu, fp nu,   fp& ps_theta, fp& ps_phi)
{
const fp tan_mu = tan(mu);
const 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
namespace local_coords
	  {
void theta_phi_of_mu_phi(fp mu, fp phi,   fp& ps_theta, fp& ps_phi)
{
fp x, y, z;
xyz_of_r_mu_phi(1.0, mu, phi,   x, y, z);
fp r;
r_theta_phi_of_xyz(x, y, z,   r, ps_theta, ps_phi);
}
	  }

//**************************************

// Bugs: computes ps_phi via trig, even though it's trivially == phi
namespace local_coords
	  {
void r_theta_phi_of_nu_phi(fp nu, fp phi,   fp& ps_theta, fp& ps_phi)
{
fp x, y, z;
xyz_of_r_nu_phi(1.0, nu, phi,   x, y, z);
fp r;
r_theta_phi_of_xyz(x, y, z,   r, ps_theta, ps_phi);
}
	  }

//******************************************************************************

namespace local_coords
	  {
void mu_nu_of_theta_phi(fp ps_theta, fp ps_phi,   fp& mu, fp& nu)
{
fp x, y, z;
xyz_of_r_theta_phi(1.0, ps_theta, ps_phi,   x, y, z);
fp r;
r_mu_nu_of_xyz(x, y, z,   r, mu, nu);
}
	  }

//**************************************

// Bugs: computes phi via trig, even though it's trivially == ps_phi
namespace local_coords
	  {
void mu_phi_of_theta_phi(fp ps_theta, fp ps_phi,   fp& mu, fp& phi)
{
fp x, y, z;
xyz_of_r_theta_phi(1.0, ps_theta, ps_phi,   x, y, z);
fp r;
r_mu_phi_of_xyz(x, y, z,   r, mu, phi);
}
	  }

//**************************************

// Bugs: computes phi via trig, even though it's trivially == ps_phi
namespace local_coords
	  {
void nu_phi_of_theta_phi(fp ps_theta, fp ps_phi,   fp& nu, fp& phi)
{
fp x, y, z;
xyz_of_r_theta_phi(1.0, ps_theta, ps_phi,   x, y, z);
fp r;
r_nu_phi_of_xyz(x, y, z,   r, nu, phi);
}
	  }

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function computes a human-readable name from a (mu,nu,phi)
// coordinates set.
//
const char *local_coords::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=0x%x isn't a valid  coords_set  bit vector!\n"
,
		   int(S));					/*NOTREACHED*/
}