#include <angles.h>
#include <mathlib.h>


doublec angles::radtodeg = 180.0/PI;
doublec angles::degtorad = 1.0/radtodeg;

doublec angles::circlemapto( doublec x, doublec y )
{
  point2<double> X(x,y);
  return circlemapto(X);
}

doublec angles::circlemapto( point2<double> const & X )
{
  doublec d = X.distance();
  if (zero<double>::test(d))
    return 0.0;

  point2<double> X2(X);
  X2 /= d;

  return circleunit(X2);
}

doublec angles::circleunit( point2<double> const & X )
{
  assert(zero<double>::test(X.x*X.x+X.y*X.y-1.0));

  double theta = abs(asin(X.y));
  if (X.x<0.0)
  {
    if (X.y<0.0)
    {
      // Quadrant 3
    
      return theta + PI;
    }
   
    // Quadrant 2
    return PI - theta;
  }

  // Quadrant 4
  if (X.y<0.0)
    return -theta+PI*2.0;
  
  // Quadrant 1
  return theta;
}

/*
// KEEP  Chelton 2007-04-09
void angles::rotationsxy
(
  double & rotx, 
  double & roty, 
  point3<double> const & X 
)
{
  roty = angles::circlemapto(X.z,X.x);
  rotx = - angles::circlemapto(X.z,X.y);
}
*/

void angles::rotateaboutaxis
(
  double & theta,
  point3<double> & axis,
  point3<double> const & nml
)
{
  static point3<double> const zaxis(0.0,0.0,1.0);

  point3<double> normal(nml);
  normal.normalize();

  // If nml coincies with the z axis no change necessary.
  if (zero<double>::test((zaxis-normal).dot()))
  {
    theta=0.0;
    return;
  }

  crossproduct::evalxyz(axis,zaxis,normal);
  theta = acos(normal.z);
}