// test_coords.cc -- test driver #2 for coordinate systems/conversions
// $Id$

#include <stdio.h>
#include <assert.h>
#include <math.h>

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

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

using jtutil::error_exit;
using jtutil::radians_of_degrees;
using jtutil::degrees_of_radians;
using namespace local_coords;

// prototypes
namespace {
void test_r_mu_nu (fp x, fp y, fp z);
void test_r_mu_phi(fp x, fp y, fp z);
void test_r_nu_phi(fp x, fp y, fp z);
void test_r_theta_phi(fp x, fp y, fp z);
	  };

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

//
// This program is a test driver for the local_coords:: coordinate-conversion
// functions.  It tries a large number of coordinate conversions, and checks
// them all for mutual consistency.
//
int main(int argc, const char *argv[])
{
bool verbose_flag = (argc == 2) && STRING_EQUAL(argv[1], "--verbose");
const fp xyz_minmax = 1.0;
const fp xyz_delta = 0.1;

	// these loops go *down* so we start in (+,+,+) quadrant
	// ==> get some successful tests before hard stuff
	// ==> help test the tests themselves
	for (fp x = xyz_minmax ; fuzzy<fp>::GE(x,xyz_minmax) ; x -= xyz_delta)
	{
	for (fp y = xyz_minmax ; fuzzy<fp>::GE(y,xyz_minmax) ; y -= xyz_delta)
	{
	for (fp z = xyz_minmax ; fuzzy<fp>::GE(z,xyz_minmax) ; z -= xyz_delta)
	{
	// avoid places where angular coords might not be defined
	if (    fuzzy<fp>::EQ(x,0.0)
	     || fuzzy<fp>::EQ(y,0.0)
	     || fuzzy<fp>::EQ(z,0.0)    )
	   then continue;

	if (verbose_flag)
	   printf("testing x=%g y=%g z=%g\n", x,y,z);

	test_r_mu_nu (x, y, z);
	test_r_mu_phi(x, y, z);
	test_r_nu_phi(x, y, z);

	test_r_theta_phi(x, y, z);
	}
	}
	}

printf("all ok!\n");
return 0;
}

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

namespace {
void test_r_mu_nu(fp x, fp y, fp z)
{
fp r, mu, nu;
r_mu_nu_of_xyz(x,y,z, r,mu,nu);

const fp phi = phi_of_mu_nu(mu,nu);
assert( fuzzy<fp>::EQ(mu, mu_of_nu_phi(nu,phi)) );
assert( fuzzy<fp>::EQ(nu, nu_of_mu_phi(mu,phi)) );

fp xx, yy, zz;
xyz_of_r_mu_nu(r,mu,nu, xx,yy,zz);

assert( fuzzy<fp>::EQ(x, xx) );
assert( fuzzy<fp>::EQ(y, yy) );
assert( fuzzy<fp>::EQ(z, zz) );
}
	  }

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

namespace {
void test_r_mu_phi(fp x, fp y, fp z)
{
fp r, mu, phi;
r_mu_phi_of_xyz(x,y,z, r,mu,phi);

const fp nu = nu_of_mu_phi(mu,phi);
assert( fuzzy<fp>::EQ(mu, mu_of_nu_phi(nu,phi)) );
assert( fuzzy<fp>::EQ(phi, phi_of_mu_nu(mu,nu)) );

fp xx, yy, zz;
xyz_of_r_mu_phi(r,mu,phi, xx,yy,zz);

assert( fuzzy<fp>::EQ(x, xx) );
assert( fuzzy<fp>::EQ(y, yy) );
assert( fuzzy<fp>::EQ(z, zz) );
}
	  }

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

namespace {
void test_r_nu_phi(fp x, fp y, fp z)
{
fp r, nu, phi;
r_nu_phi_of_xyz(x,y,z, r,nu,phi);

const fp mu = mu_of_nu_phi(nu,phi);
assert( fuzzy<fp>::EQ(nu, nu_of_mu_phi(mu,phi)) );
assert( fuzzy<fp>::EQ(phi, phi_of_mu_nu(mu,nu)) );

fp xx, yy, zz;
xyz_of_r_nu_phi(r,nu,phi, xx,yy,zz);

assert( fuzzy<fp>::EQ(x, xx) );
assert( fuzzy<fp>::EQ(y, yy) );
assert( fuzzy<fp>::EQ(z, zz) );
}
	  }

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

namespace {
void test_r_theta_phi(fp x, fp y, fp z)
{
fp r, theta, phi;
r_theta_phi_of_xyz(x,y,z, r,theta,phi);

fp mu, nu;
mu_nu_of_theta_phi(theta,phi, mu,nu);
fp theta2, phi2;
theta_phi_of_mu_nu(mu,nu, theta2,phi2);
assert( fuzzy<fp>::EQ(theta, theta2) );
assert( fuzzy<fp>::EQ(phi, phi2) );

fp xx, yy, zz;
xyz_of_r_theta_phi(r,theta,phi, xx,yy,zz);

assert( fuzzy<fp>::EQ(x, xx) );
assert( fuzzy<fp>::EQ(y, yy) );
assert( fuzzy<fp>::EQ(z, zz) );
}
	  }