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

//
// local_coords::fuzzy_EQ_{ang,dang}
// local_coords::modulo_reduce_{ang,dang}
//
// (r,(mu,nu,phi)) <--> (x,y,z)
// ((mu,nu,phi)) --> the 3rd
// usual polar spherical (r,theta,phi) <--> (x,y,z)
// ((mu,nu,phi)) <--> usual polar spherical (theta,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"
using jtutil::fuzzy;

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

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

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

//
// these functions test if two angles are fuzzily equal mod 2*pi radians
// (360 degrees)
//

namespace local_coords
	  {

bool fuzzy_EQ_ang(fp ang1, fp ang2)
{
return fuzzy<fp>::is_integer( (ang2 - ang1)/(2.0*PI) );
}

bool fuzzy_EQ_dang(fp dang1, fp dang2)
{
return fuzzy<fp>::is_integer( (dang2 - dang1)/360.0 );
}

	  }

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

//
// modulo-reduce  {ang,dang}  to be (fuzzily) within the range
// [min,max]_{ang,dang}, or error_exit() if no such value exists
//

namespace local_coords
	  {

fp modulo_reduce_ang(fp ang, fp min_ang, fp max_ang)
{
return jtutil::modulo_reduce(ang, 2.0*PI, min_ang, max_ang);
}

fp modulo_reduce_dang(fp dang, fp min_dang, fp max_dang)
{
return jtutil::modulo_reduce(dang, 360.0, min_dang, max_dang);
}

	  }

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

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

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

// not valid in xy plane (cos(mu) == 0 || cos(nu) == 0)
// c.f. page 9.2 of my mpe notes
namespace local_coords
	  {
void xyz_of_r_mu_nu(fp r, fp mu, fp nu,   fp& x, fp& y, fp& z)
{
const fp sign_y        = signum(sin(mu));
const fp sign_z_via_mu = signum(cos(mu));
assert( fuzzy<fp>::NE(cos(mu), 0.0) );
const fp y_over_z = tan(mu);

const fp sign_x        = signum(sin(nu));
const fp sign_z_via_nu = signum(cos(nu));
assert( fuzzy<fp>::NE(cos(nu), 0.0) );
const fp x_over_z = tan(nu);

// failure of next assert() ==> inconsistent input (mu,nu)
assert(sign_z_via_mu == sign_z_via_nu);
const fp sign_z = sign_z_via_mu;

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

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

assert( signum(x) == sign_x );
assert( signum(y) == sign_y );
}
	  }

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

// not valid in xz plane (sin(mu) == 0 || sin(phi) == 0)
// c.f. page 9.3 of my mpe notes
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;

const fp sign_z            = signum(sin(mu_bar));
const fp sign_y_via_mu_bar = signum(cos(mu_bar));
assert( fuzzy<fp>::NE(cos(mu_bar), 0.0) );
const fp z_over_y = tan(mu_bar);

const fp sign_x             = signum(sin(phi_bar));
const fp sign_y_via_phi_bar = signum(cos(phi_bar));
assert( fuzzy<fp>::NE(cos(phi_bar), 0.0) );
const fp x_over_y = tan(phi_bar);

// failure of next assert() ==> inconsistent input (mu,phi)
assert(sign_y_via_mu_bar == sign_y_via_phi_bar);
const fp sign_y = sign_y_via_mu_bar;

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

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

assert( signum(z) == sign_z );
assert( signum(x) == sign_x );
}
	  }

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

// not valid in yz plane (sin(nu) == 0 || cos(phi) == 0)
// c.f. page 9.4 of my mpe notes
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 ;

const fp sign_z            = signum(sin(nu_bar));
const fp sign_x_via_nu_bar = signum(cos(nu_bar));
assert( fuzzy<fp>::NE(cos(nu_bar), 0.0) );
const fp z_over_x = tan(nu_bar);

const fp sign_y         = signum(sin(phi));
const fp sign_x_via_phi = signum(cos(phi));
assert( fuzzy<fp>::NE(cos(phi), 0.0) );
const fp y_over_x = tan(phi   );

// failure of next assert() ==> inconsistent input (nu,phi)
assert(sign_x_via_nu_bar == sign_x_via_phi);
const fp sign_x = sign_x_via_nu_bar;

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

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

assert( signum(z) == sign_z );
assert( signum(y) == sign_y );
}
	  }

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

//
// 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 (cos(mu) == 0 || cos(nu) == 0)
namespace local_coords
	  {
fp phi_of_mu_nu(fp mu, fp nu)
{
fp x, y, z;
xyz_of_r_mu_nu(1.0,mu,nu, x,y,z);
return phi_of_xy(x, y);
}
	  }

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

// ill-conditioned near y axis
// not valid in xz plane (sin(mu) == 0 || sin(phi == 0)
namespace local_coords
	  {
fp nu_of_mu_phi(fp mu, fp phi)
{
fp x, y, z;
xyz_of_r_mu_phi(1.0,mu,phi, x,y,z);
return nu_of_xz(x, z);
}
	  }

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

// ill-conditioned near x axis
// not valid in yz plane (sin(nu) == 0 || cos(phi) == 0)
namespace local_coords
	  {
fp mu_of_nu_phi(fp nu, fp phi)
{
fp x, y, z;
xyz_of_r_nu_phi(1.0,nu,phi, x,y,z);
return mu_of_yz(y, z);
}
	  }


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

namespace local_coords
	  {
fp r_of_xyz(fp x, fp y, fp z)
{
return hypot3(x, y, z);
}

fp mu_of_yz(fp y, fp z)
{
return arctan_xy(z,y);
}

fp nu_of_xz(fp x, fp z)
{
return arctan_xy(z,x);
}

fp phi_of_xy(fp x, fp y)
{
return 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
	  {
fp theta_of_xyz(fp x, fp y, fp z)
{
return arctan_xy(z , hypot(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)
{
fp x, y, z;
xyz_of_r_mu_nu(1.0,mu,nu, x,y,z);

ps_theta = theta_of_xyz(x, y, z);
ps_phi = phi_of_xy(x, y);
}
	  }

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

// 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);

ps_theta = theta_of_xyz(x, y, z);
ps_phi = phi_of_xy(x, y);
assert( fuzzy_EQ_ang(ps_phi, phi) );
}
	  }

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

// Bugs: computes ps_phi via trig, even though it's trivially == phi
namespace local_coords
	  {
void 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);

ps_theta = theta_of_xyz(x, y, z);
ps_phi = phi_of_xy(x, y);
assert( fuzzy_EQ_ang(ps_phi, 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);

mu = mu_of_yz(y, z);
nu = nu_of_xz(x, z);
}
	  }

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

// 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);

mu  =  mu_of_yz(y, z);
phi = phi_of_xy(x, y);
assert( fuzzy_EQ_ang(phi, ps_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);

nu  =  nu_of_xz(x, z);
phi = phi_of_xy(x, y);
assert( fuzzy_EQ_ang(phi, ps_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  local_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,
"***** local_coords::mu_nu_phi::name_of_coords_set:\n"
"        S=0x%x isn't a valid  coords_set  bit vector!\n"
,
		   int(S));					/*NOTREACHED*/
}