// coords.hh -- classes for mpe coordinates
// $Id$
//
// coords - misc coordinates-stuff namespace
//

//
// prerequisites:
//    <jt/stdc.h> -- for DEGREES_TO_RADIANS and RADIANS_TO_DEGREES
//    fp.hh -- basic floating point stuff
//

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

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

	};