// svtest.cc -- general test of all features of SpaceVector class
//----------------------------------------------------------------

#include "ZMutility/ZMenvironment.h"

#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/AxisAngle.h"
#include "PhysicsVectors/EulerAngles.h"

#include "CLHEP/Vector/ZMxpv.h"

// #include "PhysicsVectors/UnitVector.h"
	// removed - anything really needing unitvector will be moved
	// into a different test program.

ZM_USING_NAMESPACE( zmpv )

#include "ZMutility/cmath"
using CMATH_NAMESPACE::fabs;
using CMATH_NAMESPACE::log;
using CMATH_NAMESPACE::pow;
using CMATH_NAMESPACE::tan;

#include "ZMutility/iostream"
using std::cin;
using std::cout;
using std::endl;
#include "ZMutility/iomanip"
#include "ZMutility/fstream"

//#include "CLHEP/Vector/ZMxpv.h"

USING( std::cin )
USING( std::cout )
USING( std::endl )

#include "ZMutility/PI.h"
const Float8 PI =  ZMpi;

// Modifications since 2/1/98:
// ---------------------------
//
// 2/6/98 Dictionary ordering is z first:  Changed test vectors in comparisons
//	  in dictionary ordering section.
//
// 3/19/98 Added test of eta() on lhs.  Added verification of eta formula,
//	   to verify that new optimizations are correct.
//
// 3/20/98 Added tests for nearness methods deltaR, howNear, howParallel, etc.
//
// 5/27/98 Modified exact equality tests to 10**(-15) so that x86 will pass.
//
// 6/15/98 Added namespace support
//
// 9/13/00 Elimination of v.x()=x tests, to reflect merge with CLHEP
// 9/20/00 Modifications to reflect merge with CLHEP

// Since I am not including UnitVector.h, but am using X_HAT Y_HAT Z_HAT...
static const SpaceVector X_HAT (1,0,0);
static const SpaceVector Y_HAT (0,1,0);
static const SpaceVector Z_HAT (0,0,1);

bool problems( false );

SpaceVector vans;


void constructors();
void accessors();
void modifiers();
void set3();
void arithmetic();
void comparisons();
void tolerance();
void nearness();
void combining();
void intrinsic();
void input();
void direction();
void anglebetween();
void angledecomposition();
void rotations();
void rotations2();
void rotations3();


bool checkAxialRotation ( const SpaceVector &x,
			  const SpaceVector &d,
			  const Float8 &phi,
			  const SpaceVector& result );

SpaceVector goldstein107 (const SpaceVector& v,
			  Float8 phi, Float8 theta, Float8 psi);

bool near ( double a, double b, double tol ) {
  return ( fabs(a-b) < tol );
}
bool close ( double a, double b ) {
  return ( fabs(a-b) < 1.0E-10 );
}
bool eq ( double a, double b ) {
  return ( fabs(a-b) < 1.0E-14 );
}
const double teeny = 1.0E-14;


int main() {

// For purposes of this test, it is safe and necessary to ignore ZMexceptions.
// This test exercises some of them but we want to continue beyond each one.
  ZMxPhysicsVectors::setHandler( ZMexIgnoreAlways() );

  cout << "\n\n";
  cout << "   ***********\n";
  cout << "   SpaceVector\n";
  cout << "   ***********\n";

  constructors();
  accessors();
  modifiers();
  set3();

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
  // Pathological cases
#endif //  NOTYET_WAIT_FOR_EXCEPTIONS

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
// Boundaries of the range for constructors
  cout << "\n --- Ranges for constructors\n\n";
#endif // NOTYET_WAIT_FOR_EXCEPTIONS

  arithmetic();
  comparisons();
  tolerance();
  nearness();
  combining();
  intrinsic();
  input();
  direction();
  anglebetween();
  angledecomposition();
  rotations();
  rotations2();
  rotations3();

  if (problems) {
    cout << "\n ?????????????????? -- End of tests:  PROBLEMS DETECTED\n\n";
    return -1;
  } else {
    cout << "\nEnd of tests -- All automated checking passed.\n\n";
    return 0;
  }

} /* End of svtest main() */

bool checkAxialRotation ( const SpaceVector &x,
			  const SpaceVector &d,
			  const double& phi,
			  const SpaceVector& result ) {

  // Rotate x by phi around d; return true if result does not match answer
  // within .0001 tolerance; print out appropriate messages.

  SpaceVector xx (x);
  xx.rotate(d, phi);
  if ( xx.isNear(result, .0001) ) {
    cout << "OK: " << x << "Rotated " << phi << " about" << d <<
	"\n      = " << xx << " ~ " << result << endl;
  } else {
    cout << "?????" << x << "Rotated " << phi << " about" << d <<
	"\n      = " << xx << " not near " << result << endl;
    return true;
  }

  // Now form a new SpaceVector as xx rotated about the axis by -phi;
  // again check that the result is proper (matches x) very tightly

  SpaceVector yy = rotationOf(xx, d, -phi);
  if ( yy.isNear(x, .000001) ) {
    cout << "OK: " << xx << "Rotated " << -phi << " about" << d <<
	"\n      = " << yy << " ~ " << x << endl;
  } else {
    cout << "?????" << xx << "Rotated " << -phi << " about" << d <<
	"\n      = " << yy << " not near " << x << endl;
    return true;
  }
  return false;

  // How did we get those results that we have supplied?
  //
  // We rotate the vector in the following manner:
  // (1) From a set of three coordinate-like vectors with d used as the third:
  //    a - choose V1 = d cross x
  //    b - choose V2 = d cross V1
  //	c - the third is just d
  // (2) Decompose x into u V1 + v V2 + w d
  //	You know a priori that u = 0
  //	Unless d is purely in the Z direction, the two equations dealing with
  // 	the x and y directions suffice:
  //	v = (x.x*d.y - d.x*x.y) / ( V2.x*d.y - d.x*V2.y ) and similar for w.
  // (3) express x in terms of V2-hat instead of V2
  // (4) Rotate in V2-hat, V1hat space by phi
  // (5) rescale back to V1, V2
  // (6) transform back using V1, V2, and d to the x,y,z coordinates.
  // Steps 3--6 combine:
  // x' = - v sin(phi) (|V2| \ |V1|) V1 + v cos(phi) V2 + w d
  //
  // That this result matches the formula used in SpaceVector.cc is a highly
  // non-trivial test.

} /* checkAxialRotation */

SpaceVector goldstein107 (const SpaceVector& v,
			  Float8 phi, Float8 theta, Float8 psi) {

  // Perform an Euler angles rotation of v by laboriously folloing the
  // illustration on page 107 of Goldstein Classical Mechanincs.  Then
  // check into the booby hatch!

  // Note that our rotations about an axis are active, but Goldstein uses
  // (passive) frame transformations, which is the only good way to define
  // Euler angles.

  // Step 1 (diagram a):
  //
  // Rotate the X axis by phi around the Z axis to form a new axis
  // which Goldstein calls xi.  In this new frame, v is the old v rotated
  // by -phi about Z to form v'.

// #define TRACE_GOLD

  SpaceVector xi = rotationOf (X_HAT, Z_HAT, phi);
  SpaceVector eta = rotationOf (Y_HAT, Z_HAT, phi);
  SpaceVector zeta = rotationOf (Z_HAT, Z_HAT, phi);

  SpaceVector v_prime =  rotationOf (v, Z_HAT, -phi);

#ifdef TRACE_GOLD
  cout << "        goldstein107:  V is " << v << endl;
  cout << "        new frame is xi, eta, zeta =\n"
       << "        " << xi << endl
       << "        " << eta << endl
       << "        " << zeta << endl;
  cout << "        in that frame V has components\n"
       << "        " << v.dot(xi) << " , "
		     << v.dot(eta) <<" , " << v.dot(zeta) << endl;
  cout << "        and V' is " << v_prime << endl;
#endif

  // Step 2 (diagram b):
  //
  // Rotate the Z or zeta  axis by theta around that xi axis to form a new axis
  // which Goldstein calls zeta_prime.  Since in that frame a rotation about the
  // xi axis has the form we are familiar with as a RotationX, we rotate about
  // v about X_HAT to get the effect of this rotation.  Thus in the xi_prime,
  // eta_prime, zeta_prime coordinate frame, v'' is v', X-rotated by -theta.

  SpaceVector zeta_prime = rotationOf (zeta, xi, theta);
  SpaceVector eta_prime = rotationOf (eta, xi, theta);
  SpaceVector xi_prime = rotationOf (xi, xi, theta);

  SpaceVector v_double_prime = rotationOf (v_prime, X_HAT, -theta);

#ifdef TRACE_GOLD
  cout << "        Step 2" << endl;
  cout << "        new frame is xi', eta', zeta' =\n"
       << "        " << xi_prime << endl
       << "        " << eta_prime << endl
       << "        " << zeta_prime  << endl;
  cout << "        in that frame V has components\n"
       << "        " << v.dot(xi_prime) << " , "
		     << v.dot(eta_prime) <<" , " << v.dot(zeta_prime) << endl;
  cout << "        and V'' is " << v_double_prime << endl;
#endif

  // Step 3 (diagram c):
  //
  // Rotate the xi (xi-prime) axis by psi around the new zeta-prime axis.
  // Once again, the effect on V is a rotation using the simple RotateZ form
  // in those  primed axes.  This rotation produces the new coordinate frame
  // labeled x-prime, (the rotated xi-prime) y-prime (the rotated eta-prime)
  // and z-prime (same as zeta-prime or zeta).
  // We find our final v''' by rotating v'' by -psi around the Z axis.

  SpaceVector x_prime = rotationOf (xi_prime, zeta_prime, psi);
  SpaceVector y_prime = rotationOf (eta_prime, zeta_prime, psi);
  SpaceVector z_prime = rotationOf (zeta_prime, zeta_prime, psi);

  SpaceVector result = rotationOf (v_double_prime, Z_HAT, -psi );

#ifdef TRACE_GOLD
  cout << "        Step 3" << endl;
  cout << "        new frame is x', y', z' =\n"
       << "        " << x_prime << endl
       << "        " << y_prime << endl
       << "        " << z_prime << endl;
  cout << "        in that frame V has components\n"
       << "        " << v.dot(x_prime) << " , "
		     << v.dot(y_prime) << " , " << v.dot(z_prime) << endl;
  cout << "        and result is " << result << endl;
#endif

  return result;

} /* goldstein107 */


      SpaceVector v1 ( 1., 2., 3. );
const SpaceVector v1c ( 2, 3., 4. );
      SpaceVector v2 ( 10, 30, DEGREES, 60, DEGREES );
const SpaceVector v2c ( 20, 30, DEGREES, PI/3, RADIANS );
      SpaceVector v3 ( 20, 90, DEGREES, -15 );


void constructors()  {

  cout << "\n --- Constructors and Output\n\n";

	cout << "{ Cartesian Constructor }\n";

  if ( (v1.x() == 1.) && (v1.y()==2) && 3==v1.z() ) {
    cout << "OK: constructed" << v1 << endl;
  } else {
    cout << "????? constructed" << v1 << endl;
    problems=true;
  }

  if ( sizeof(v1) == 24 ) {
    cout << "OK: sizeof SpaceVector is " << sizeof (v1) << endl;
  } else {
    cout << "????? sizeof SpaceVector is " << sizeof (v1) << endl;
    problems=true;
  }

  if ( (v1c.z() == 4.) && (v1c.x()==2) && 3==v1c.y() ) {
    cout << "OK: constructed const " << v1c << endl;
  } else {
    cout << "????? constructed const " << v1c << endl;
    problems=true;
  }

	cout << "{ Spherical Constructor }\n";

  if ( v2.isNear(SpaceVector(2.5, 4.33, 8.66),.01) ) {
    cout << "OK: constructed " << v2.r() << ", " << v2.theta()*(180/PI) <<
	" degrees, " << v2.phi()*(180/PI) << " degrees = " << v2 << endl;
  } else {
    cout << "????? constructed " << v2.r() << ", " << v2.theta()*(180/PI) <<
	", " << v2.phi()*(180/PI) << " = " << v2 << endl;
    problems=true;
  }

  if ( v2c.isNear(SpaceVector(5.0, 8.66, 17.32),.01) ) {
    cout << "OK: constructed const " << v2c.r() << ", "
	<< v2c.theta()*(180/PI) << " degrees, " << v2c.phi()*(180/PI)
	<< " degrees = " << v2c << endl;
  } else {
    cout << "????? constructed const " << v2c.r() << ", "
	<< v2c.theta()*(180/PI) << ", " << v2c.phi()*(180/PI)
	<< " = " << v2c << endl;
    problems=true;
  }

  SpaceVector va2 ( 10, PI/6, RADIANS, -60, DEGREES );
  if ( va2.isNear(SpaceVector(2.5, -4.33, 8.66),.01) ) {
    cout << "OK: constructed " << va2.r() << ", " << va2.theta()*(180/PI) <<
	" degrees, " << va2.phi()*(180/PI) << " degrees = " << va2 << endl;
  } else {
    cout << "????? constructed " << va2.r() << ", " << va2.theta()*(180/PI) <<
	", " << va2.phi()*(180/PI) << " = " << va2 << endl;
    problems=true;
  }

  const SpaceVector va2c ( 20, PI/6, RADIANS, -PI/3, RADIANS );
  if ( va2c.isNear(SpaceVector(5.0, -8.66, 17.32),.01) ) {
    cout << "OK: constructed const " << va2c.r() << ", "
	<< va2c.theta()*(180/PI) << " degrees, " << va2c.phi()*(180/PI)
	<< " degrees = " << va2c << endl;
  } else {
    cout << "????? constructed const " << va2c.r() << ", "
	<< va2c.theta()*(180/PI) << ", " << va2c.phi()*(180/PI)
	<< " = " << va2c << endl;
    problems=true;
  }

  SpaceVector vb2 ( 10, 1.31696, ETA, 60, DEGREES );
  if ( vb2.isNear(SpaceVector(2.5, 4.33, 8.66),.01) ) {
    cout << "OK: constructed " << vb2.r() << ", eta = "
	<< vb2.eta()*(180/PI) << " --> " << vb2.theta()*(180/PI) <<
	" degrees, " << vb2.phi()*(180/PI) << " degrees = " << vb2 << endl;
  } else {
    cout << "????? constructed " << vb2.r() << ", " << vb2.theta()*(180/PI) <<
	", " << vb2.phi()*(180/PI) << " = " << vb2 << endl;
    problems=true;
  }

  const SpaceVector vb2c ( 20, 1.31696, ETA, PI/3, RADIANS );
  if ( vb2c.isNear(SpaceVector(5.0, 8.66, 17.32),.01) ) {
    cout << "OK: constructed const " << vb2c.r() << ", "
	<< vb2c.theta()*(180/PI) << " eta=" << vb2c.eta() <<
	", " << vb2c.phi()*(180/PI)
	<< " = " << vb2c << endl;
  } else {
    cout << "????? constructed const " << vb2c.r() << ", "
	<< vb2c.theta()*(180/PI) << " eta=" << vb2c.eta() <<
	", " << vb2c.phi()*(180/PI)
	<< " = " << vb2c << endl;
    problems=true;
  }

	cout << "{ Cylindrical Constructor }\n";

  if (v3.isNear(SpaceVector(0, 20, -15))) {
    cout << "OK: constructed  " << v3 << endl;
  } else {
    cout << "????? constructed  " << v3 << endl;
  }

  SpaceVector va3 ( 20, -PI/2, RADIANS, -15 );
  if (va3.isNear(SpaceVector(0, -20, -15))) {
    cout << "OK: constructed  " << va3 << endl;
  } else {
    cout << "????? constructed  " << va3 << endl;
    problems=true;
  }

	cout << "{ Array Constructor }\n";
  SpaceVector arrayV [100];
  cout << "OK: created array of SpaceVectors" << endl;

}  // constructors()


SpaceVector v4 (3, 4, 12);
SpaceVector v5 ( 5, -12, 0 );
SpaceVector v6 (1, 2, 1000);

double xx, yy, zz;



void accessors()  {

// Accessors

  // Straight x y z

  cout << "\n --- Accessors\n\n";

	cout << "{ get Cartesian components }\n";
  xx = v1.x(); yy = v1.y(); zz = v1.z();
  if ( (xx==1) && (yy==2) && (zz==3) ) {
    cout << "OK: v1 = " << v1 << endl;
  } else {
    cout << "????? v1 = " << v1
	 << " but xx,yy,zz = " << xx << ',' << yy << ',' << zz << endl;
    problems=true;
  }

	cout << "{ output Cartesian components }\n";
  cout << "    ... ostream<< x, y, z  ==> "
			<< v1.x() <<"," << v1.y() <<"," << v1.z() << endl;

  xx = v1c.x();
  yy = v1c.y();
  zz = v1c.z();
  if ( (xx==2) && (yy==3) && (zz==4) ) {
    cout << "OK: v1c = " << v1c << endl;
  } else {
    cout << "????? v1c = " << v1c << endl;
    problems=true;
  }
  cout << "    ... ostream<< x, y, z  ==> "
			<< v1c.x() <<"," << v1c.y() <<"," << v1c.z() << endl;

	cout << "{ use Cartesian components in expressions }\n";

  if ( 6*v1.z() == 18 ) {
    cout << "OK: 6*v1.z() = " << 6*v1.z() << endl;
  } else {
    cout << "????? 6*v1.z() = " << 6*v1.z() << endl;
    problems=true;
  }

  if ( v1.z()+4 == 7 ) {
    cout << "OK: v1.z()+4 = " << v1.z()+4 << endl;
  } else {
    cout << "????? v1.z()+4 = " << v1.z()+4 << endl;
    problems=true;
  }

  if ( 5*v1c.y() == 15 ) {
    cout << "OK: 5*v1c.z() = " << 5*v1c.z() << endl;
  } else {
    cout << "????? 5*v1c.z() = " << 5*v1c.z() << endl;
    problems=true;
  }

  if ( v1c.y()/2 == 1.5 ) {
    cout << "OK: v1c.y()/2 = " << v1c.y()/2 << endl;
  } else {
    cout << "????? v1c.y()/2 = " << v1c.y()/2 << endl;
    problems=true;
  }

  // Observation:  v.dx directly does not compile, as desired.
#ifdef ZAP1
  cout << v.dx;
#endif // ZAP1


  // Comparison of accessors

	cout << "{ compare Cartesian components }\n";
  if ( v1.x() == 1  ) {
    cout << "OK: Equality comparison on an accessor works\n";
  } else {
    cout << "????? Equality comparison on an accessor FAILS!!!!!\n";
    problems=true;
  }

  if ( v1.x() != v1.y() ) {
    cout << "OK: comp != comp works\n";
  } else {
    cout << "????? comp != comp FAILS!!!!\n";
    problems=true;
  }

  if ( v1.x() == v1.y() ) {
    cout << "????? comp == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp == comp works\n";
  }

  if ( v1.x() > v1.y() ) {
    cout << "?????comp > comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp > comp works\n";
  }


  if ( v1.y() < v1.x() ) {
    cout << "????? comp < comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp < comp works\n";
  }


  double xeq = 4.0;

  if ( xeq == v1.z() ) {
    cout << "????? Scalar == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: Scalar == comp works\n";
  }

  const double ceq = 4.0;
  if ( ceq == v1.z() ) {
    cout << "????? Scalar == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: Scalar == comp works\n";
  }

  if ( v1c.x() == 2 ) {
    cout << "OK: Equality comparison on a const vector accessor works\n";
  } else {
    cout << "?????Equality comparison on a const vector accessor FAILS!!!!!\n";
    problems=true;
  }

  if ( v1.x() != v1.y() ) {
    cout << "OK: comp != comp works\n";
  } else {
    cout << "?????comp != comp FAILS!!!!\n";
    problems=true;
  }

  if ( v1.x() == v1.y() ) {
    cout << "????? comp == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp == comp works\n";
  }

  if ( v1.x() > v1.y() ) {
    cout << "????? comp > comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp > comp works\n";
  }

  if ( v1.y() < v1.x() ) {
    cout << "????? comp < comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp < comp works\n";
  }

  xeq = 4.0;


  if ( v1c.y() < v1c.z() ) {
    cout << "OK: const comp == const comp works\n";
  } else {
    cout << "????? const comp = const comp FAILS!!!!\n";
    problems=true;
  }



  // Getting different components

	cout << "{ get Spherical components }\n";

  double p (v2.theta());
  if ( close(p,PI/6) ) {
    cout << "OK: v2.theta() = " << v2.theta()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2.theta() should be " << PI/6 <<
	" but = " << v2.theta()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2.theta() ==> " << v2.theta() << endl;

  double a = v2.phi();
  if ( close(a,PI/3) ) {
    cout << "OK: v2.phi() = " << v2.phi()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2.phi() should be " << PI/3 <<
	" but = " << v2.phi()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2.phi() ==> " << v2.phi() << endl;

  double pp = v2c.theta();
  if ( close(pp,PI/6) ) {
    cout << "OK: v2c.theta() = " << v2c.theta()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2c.theta() should be " << PI/6 <<
	" but = " << v2c.theta()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2c.theta() ==> " << v2c.theta() << endl;

  double ap (v2c.phi());
  if ( close(ap,PI/3) ) {
    cout << "OK: v2c.phi() = " << v2c.phi()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2c.phi() should be " << PI/3 <<
	" but = " << v2c.phi()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2c.phi() ==> " << v2c.phi() << endl;

  p = v2c.theta();
  if ( close(p,PI/6) ) {
    cout << "OK: v2c.theta() = " << v2c.theta()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2c.theta() should be " << PI/6 <<
	" but = " << v2c.theta()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2c.theta() ==> " << v2c.theta() << endl;

  a = v2.phi();
  if ( close(a,PI/3) ) {
    cout << "OK: v2.phi() = " << v2.phi()*180/PI << " degrees " << endl;
  } else {
    cout << "????? v2.phi() should be " << PI/3 <<
	" but = " << v2.phi()*180/PI << " degrees " << endl;
    problems=true;
  }
  cout << "    .... ostream << v2.phi() ==> " << v2.phi() << endl;

  xx =  v4.r();
  if ( xx == 13 ) {
    cout << "OK: " << v4 << ".r() ==> " << xx << endl;
  } else {
    cout << "????? " << v4 << ".r() ==> " << xx << endl;
    problems=true;
  }
  cout << "    .... ostream << v4.r() ==> " << v4.r() << endl;

  if (!v4.isNear( SpaceVector (3, 4, 12) )) {
    cout << "????? getting r changed the vector! " << v4 << endl;
    problems=true;
  }

  xx =  v4.rho();
  if ( xx == 5 ) {
    cout << "OK: " << v4 << ".rho() ==> " << xx << endl;
  } else {
    cout << "????? " << v4 << ".rho() ==> " << xx << endl;
    problems=true;
  }
  cout << "    .... ostream << v4.rho() ==> " << v4.rho() << endl;

  if (!v4.isNear( SpaceVector (3, 4, 12) )) {
    cout << "????? getting rho changed the vector! " << v4 << endl;
    problems=true;
  }


	cout << "{ Compare Spherical components }\n";



  if ( v6.theta() < 2*PI/180 ) {
    cout << "OK: "<< v6.theta() << " compares as < 2 degrees" << endl;
  } else {
    cout << "?????" << v6.theta() << " compares as not < 2 degrees!!!!" << endl;
    problems=true;
  }


	cout << "{ Compare Spherical components of const }\n";

  if ( v2c.phi() >= 59*PI/180 ) {
    cout << "OK: "<< v2c.phi() << " compares as >= 59 degrees" << endl;
  } else {
    cout << "?????" << v2c.phi() << " compares as not >= 59 degrees!" << endl;
    problems=true;
  }


	cout << "{ eta }\n";

  cout << v5 << " has theta, phi of "
	<< v5.theta() << " " << v5.phi() << endl;
  xx = v5.eta();
  if (eq(xx,0))   {
    cout << "OK: " << v5 << " has eta = " << xx << endl;
  } else {
    cout << "?????" << v5 << " has eta = " << xx << endl;
    problems=true;
  }

  SpaceVector vtest ( 5, 3, 6.751 );
  double testeta = - log( tan ( .5*vtest.theta() ) );
  if ( close (vtest.eta(), testeta ) ) {
    cout << "OK: " << vtest << " has eta = " << vtest.eta() << endl;
  } else {
    cout << "????? " << vtest << " has eta = " << vtest.eta() <<
		" != " << testeta << endl;
    problems=true;
  }

  cout << "    ...   eta for (1, 2, 1000) = " << v6.eta() << endl;
  cout << "    ... theta for (1, 2, 1000) = " << v6.theta() << endl;


  // Setting different components

	cout << "{ Set Cartesian components }\n";

  v6.setX(24);
  v6.setY(3.);
  v6.setZ(100);
  vans.set(24,3,100);
  if (v6==vans) {
    cout << "OK: set x=24, y=3, z=100 ==> " << v6 << endl;
  } else {
    cout << "OK: set x=24, y=3, z=100 ==> " << v6 << endl;
    problems=true;
  }
	cout << "{ Set Spherical components }\n";

  double p4 = v4.theta();
  double a4 = v4.phi();
  cout << "    ... Theta for " << v4 << "is   " << p4 << endl;
  cout << "    ... Phi   for " << v4 << "is   " << a4 << endl;
  v4.setTheta(0);
  if (v4.isNear(SpaceVector(0,0,13))) {
    cout << "OK: with v4.theta() = 0 v4 is " << v4 << endl;
  } else {
    cout << "????with v4.theta() = 0 v4 is " << v4 << endl;
    problems=true;
  }

  cout << "    " << v4 << "    ... Reset theta and then phi:\n     ";
  v4.setTheta(p4);
  cout << v4;
  v4.setPhi(a4);
  cout << v4 << endl;
  vans.set (3, 4, 12);
  if (v4.isNear(vans)) {
    cout << "OK: " << v4 << " = " << vans << endl;
  } else {
    cout << "????? " << v4 << " is not " << vans << endl;
    problems = true;
  }

  cout << v4 <<"    ... Now set r to 26:\n";
  v4.setR(26);
  vans.set (6, 8, 24);
  if (v4.isNear(vans)) {
    cout << "OK: " << v4 << " = " << vans << endl;
  } else {
    cout << "????? " << v4 << " is not " << vans << endl;
    problems = true;
  }

  // ADDED 980319 MF:
  SpaceVector ve ( 4, 3, 12 );
  ve.setEta(2.5);
  if ( close(ve.eta(),2.5) && close(ve.r(),13) ) {
    cout << "OK: " << ve << " has eta = " << ve.eta()
		<< " and r = " << ve.r() << endl;
  } else {
    cout << "????? " << ve << " has eta = " << ve.eta()
		<< " and r = " << ve.r() << endl;
    problems = true;
  }

#ifdef ZAP2
  v2c.setTheta(.05);
  // Properly fails to compile 
#endif

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
  // Setting components, with vector along Z axis (pos and neg)
#endif // NOTYET_WAIT_FOR_EXCEPTIONS

  // setCylTheta, setCylEta

	cout << "{ Set cylindrical theta }\n";

  v3.setPhi( -PI/6);
  cout << "    ... Vector starts as " << v3 << "change theta to 115 ==>\n";
  v3.setCylTheta(115*PI/180);
  cout << v3 << endl;
  cout << "      rho = " << v3.rho() << " theta = " << v3.theta() <<
  		" phi = " << v3.phi() << endl;
  vans.setCyl ( 20, (-PI/6), 115*PI/180 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v3 << " = " << vans << endl;
  } else {
    cout << "????? " << v3 << " is not " << vans << endl;
    problems = true;
  }

	cout << "{ Set cylindrical eta }\n";

  v3.setCylEta(-1.5);
  cout << "    ... Vector starts as " << v3 << "change CylEta to -1.5 ==>\n";
  vans.setCyl ( 20, (-PI/6), -1.5, ETA );
  cout << "     rho = " << v3.rho() << " phi = " << v3.phi() <<
		" eta = " << v3.eta() << endl;
  if (v3.isNear(vans)) {
    cout << "OK: " << v3 << " = " << vans << endl;
  } else {
    cout << "????? " << v3 << " is not " << vans << endl;
    problems = true;
  }

}  // accessors()



void modifiers()  {

  cout << "\n --- Component Modify/Assigners\n\n";

// Component modifiers ( v.x() += Scalar, etc.)

	cout << "{ Modifying cartesian components }\n";

  cout << v1 << "Add 5 to X component\n";
  v1.setX(v1.x()+5);
  vans.set(6,2,3);
  if (v1.isNear(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << "v1.y() *= v1.x\n";
  v1.setY( v1.y() * v1.x() );
  vans.set(6,12,3);
  if (v1.isNear(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << "Subtract 5 from Z component\n";
  v1.setZ(v1.z()-5);
  vans.set(6,12,-2);
  if (v1.isNear(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }


  cout << v1 << "v1.x() /= 3\n";
  v1.setX(v1.x()/3);
  vans.set(2,12,-2);
  if (v1.isNear(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

	cout << "{ Component arithmetic }\n";

  if ( (vans.x()+5) == 7 ) {
    cout << "OK: vans.x()+5 = " << vans.x()+5 << endl;
  } else {
    cout << "????? vans.x()+5 = " << vans.x()+5 << endl;
    problems = true;
  }


	cout << "{ Modifying Spherical components }\n";

  v2.setR(10);
  v2.setTheta(PI/6);
  v2.setPhi(-PI/6);

  cout << " ( " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees ";
  cout << " v2.phi() += 50 degrees\n";
  v2.setPhi(v2.phi() + 50*PI/180);
  vans.set ( 10, 30, DEGREES, 20, DEGREES);
  if (v2.isNear(vans)) {
    cout << "OK : " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees " << endl;
  } else {
    cout << "?????: " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees " << endl;
    problems = true;
  }

  cout << " ( " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees ";
  cout << " v2.theta() *= 2\n";
  v2.setTheta(v2.theta() * 2);
  vans.set ( 10, 60, DEGREES, 20, DEGREES);
  if (v2.isNear(vans)) {
    cout << "OK : " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees "<< endl;
  } else {
    cout << "?????: " << v2.r() << ", " << v2.theta()*180/PI << " degrees " <<
	v2.phi()*180/PI << " degrees "<< endl;
    problems = true;
  }

  cout << v1 << "v1.r() *= 3\n";
  v1.setR(v1.r()*3);
  vans.set ( 6, 36, -6);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

  cout << v1 << "v1.r() += sqrt(1368.)\n";
  v1.setR(v1.r()+sqrt(1368.));
  vans.set ( 12, 72, -12);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

  cout << v1 << "v1.r() /= 2\n";
  v1.setR(v1.r()/2);
  vans.set ( 6, 36, -6);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

  v1.set (2, 12, -2);

  cout << v1 << "v1.r() *= 3\n";
  v1.setR(v1.r()*3);
  vans.set ( 6, 36, -6);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

  cout << v1 << "v1.r() *= 2; v1.r() -= sqrt(1368.);\n";
  v1.setR(v1.r()*2);
  v1.setR(v1.r()-sqrt(1368.));
  vans.set ( 6, 36, -6);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
	cout << "{ Exceptionally Modifying Spherical components }\n";
  SpaceVector vzip (0,0,0);
  vzip.r() *= 0;
  cout << "vzip *= 0 works fine : vzip = " << vzip << endl;
  vzip.r() += 6;
  cout << "vzip += 6 gives vzip = " << vzip << endl;
  vzip.r() /= 0;
  cout << "vzip /= 0 gives above exception and returns vzip = " << vzip << endl;
#endif // NOTYET_WAIT_FOR_EXCEPTIONS


	cout << "{ Modifying Cylindrical rho component }\n";

  cout << v1 << "v1.rho() *= 2\n";
  v1.setRho(v1.rho()*2);
  vans.set ( 12, 72, -6);
  if (v1.isNear(vans)) {
    cout << "OK : " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }

}  // modifiers()



void set3()  {

  cout << "\n --- Setting three components at once\n\n";

	cout << "{ Cartesian components }\n";

  v1.set ( 5, 7, 9 );
  if ( (v1.x() == 5) && (v1.y() == 7) && (v1.z() == 9) ) {
    cout << "OK: " << v1 << endl;
  } else {
    cout << "?????" << v1 << endl;
    problems = true;
  }
	cout << "{ Spherical components }\n";

  v1.set ( 5, 42, DEGREES, 99, DEGREES );
  if ( close(v1.r(),5) &&
	close(v1.theta(),42*PI/180) &&
	close(v1.phi(),99*PI/180) ) {
    cout << "OK: " << v1.r() << ", " << v1.theta()*180/PI << " degrees " <<
	v1.phi()*180/PI << " degrees \n";
  } else {
    cout << "?????" << v1.r() << ", " << v1.theta()*180/PI << " degrees " <<
	v1.phi()*180/PI << " degrees \n";
  }

  v1.set ( 5, 1.2, RADIANS, -467, DEGREES );
  if ( close(v1.r(),5) &&
	close(v1.theta(),1.2) &&
	close(v1.phi(),-107*PI/180) ) {
    cout << "OK: " << v1.r() << ", " << v1.theta() << " radians " <<
	v1.phi()*180/PI << " degrees \n";
  } else {
    cout << "?????" << v1.r() << ", " << v1.theta() << " radians " <<
	v1.phi()*180/PI << " degrees \n";
    problems = true;
  }

  v1.set ( 5, 3, ETA, .5, RADIANS );
  if ( close(v1.r(),5) &&
	close(v1.eta(),3) &&
	close(v1.phi(),.5) ) {
    cout << "OK: " << v1.r() << ", eta: " << v1.eta() << " " <<
	v1.phi() << " radians \n";
  } else {
    cout << "?????" << v1.r() << ", eta: " << v1.eta() << " " <<
	v1.phi() << " radians \n";
    problems = true;
  }


	cout << "{ Cylindrical components }\n";
  v1.set ( 5, -17, DEGREES, 9 );
  if ( close(v1.rho(),5) &&
	close(v1.phi(),-17*PI/180) &&
	close(v1.z(),9) ) {
    cout << "OK: " << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees " << v1.z() << endl;
  } else {
    cout << "?????" << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees " << v1.z() << endl;
    problems = true;
  }

	cout << "{ Cylindrical with theta or eta determined }\n";
  v1.setCyl ( sqrt(2.), -PI/4, PI/4 );
  vans.set (1, -1, sqrt(2.));
  if (v1.isNear(vans)) {
    cout << "OK: " << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees, " <<
	v1.theta()*180/PI  << " degrees \n";
  } else {
    cout << "?????" << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees, " <<
	v1.theta()*180/PI  << " degrees \n";
    problems = true;
  }

  v1.setCyl ( 1, PI/2, -7, ETA );
  // Note: If we take eta too big, e.g. = -10, then exp(eta) loses enough
  // acuracy that our "close()" function (which requires for a value of 10
  // accuracy to a part in 10**(-11)) will say it is not close!
  if ( close (v1.x(),0) &&
	close(v1.y(),1) &&
	close(v1.eta(),-7) ) {
    cout << "OK: " << v1 << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees, eta = " <<
	v1.eta() << " z = " << v1.z() << endl;
  } else {
    cout << "????? " << v1 << v1.rho() << ", "<<
	v1.phi()*180/PI << " degrees, eta = " <<
	v1.eta() << " z = " << v1.z() << endl;
    problems = true;
  }

}  // set3()


const SpaceVector v3c( 3, 5, 7);
const SpaceVector v4c( 6, 2, 8);



void arithmetic()  {

// Arithmetic operations

  cout << "\n --- Arithmetic\n\n";


	cout << "{ + and - }\n";

  v1.set( 3, 5, 7);
  v2.set( 6, 2, 8);

  v3 = v1 + v2;
  vans.set ( 9, 7, 15 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v1 << "+" << v2 << "=" << v3 << endl;
  } else {
    cout << "?????" << v1 << "+" << v2 << "=" << v3 << endl;
    problems = true;
  }

  v3 = v1 - v2;
  vans.set ( -3, 3, -1 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v1 << "-" << v2 << "=" << v3 << endl;
  } else {
    cout << "?????" << v1 << "-" << v2 << "=" << v3 << endl;
    problems = true;
  }

  v3 = - v2;
  vans.set ( -6, -2, -8 );
  if (v3.isNear(vans)) {
    cout << "OK: " << "-" << v2 << "=" << v3 << endl;
  } else {
    cout << "?????" << "-" << v2 << "=" << v3 << endl;
    problems = true;
  }

  v3 = v4c + v3c;
  vans.set ( 9, 7, 15 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << "+" << v3c << "=" << v3 << endl;
  } else {
    cout << "?????" << v4c << "+" << v3c << "=" << v3 << endl;
    problems = true;
  }

  v3 = v4c - v3c;
  vans.set ( 3, -3, 1 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << "-" << v3c << "=" << v3 << endl;
  } else {
    cout << "?????" << v4c << "-" << v3c << "=" << v3 << endl;
    problems = true;
  }

  v3 = - v3c;
  vans.set ( -3, -5, -7 );
  if (v3.isNear(vans)) {
    cout << "OK: " << "-" << v3c << "=" << v3 << endl;
  } else {
    cout << "?????" << "-" << v3c << "=" << v3 << endl;
    problems = true;
  }

	cout << "{ * and / }\n";

  double pointfive = .5;
  v3 = v2 * pointfive;
  vans.set ( 3, 1, 4 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v2 << "* .5 =" << v3 << endl;
  } else {
    cout << "?????" << v2 << "* .5 =" << v3 << endl;
    problems = true;
  }

  v3 = v2 / 4;
  vans.set ( 1.5, .5, 2 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v2 << "/ 4 =" << v3 << endl;
  } else {
    cout << "?????" << v2 << "/ 4 =" << v3 << endl;
    problems = true;
  }

  v3 = 1.5 * v2;
  vans.set ( 9, 3, 12 );
  if (v3.isNear(vans)) {
    cout << "OK: 1.5 *" << v2 << "=" << v3 << endl;
  } else {
    cout << "????? 1.5 *" << v2 << "=" << v3 << endl;
    problems = true;
  }

  v3 = v4c * .5;
  vans.set ( 3, 1, 4 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << "* .5 =" << v3 << endl;
  } else {
    cout << "?????" << v4c << "* .5 =" << v3 << endl;
    problems = true;
  }

  v3 = v4c / .5;
  vans.set ( 12, 4, 16 );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << "/ .5 =" << v3 << endl;
  } else {
    cout << "?????" << v4c << "/ .5 =" << v3 << endl;
    problems = true;
  }

  v3 = 2.5 * v4c;
  vans.set (15, 5, 20 );
  if (v3.isNear(vans)) {
    cout << "OK: 2.5 *" << v4c << "=" << v3 << endl;
  } else {
    cout << "????? 2.5 *" << v4c << "=" << v3 << endl;
    problems = true;
  }

	cout << "{ Unary minus }\n";

  v3 = -v1;
  if (v3 == SpaceVector (-3, -5, -7)) {
    cout << "OK: - of " << v1 << "\n    is " << v3 << endl;
  } else {
    cout << "????? - of " << v1 << "\n    is " << v3 << endl;
    problems = true;
  }

  v3 = -v4c;
  if (v3 == SpaceVector (-6, -2, -8)) {
    cout << "OK: - of " << v4c << "\n    is " << v3 << endl;
  } else {
    cout << "????? - of " << v4c << "\n    is " << v3 << endl;
    problems = true;
  }

	cout << "{ Linear expression }\n";

  v3.set(9,3,12);
  v4 = v3 + 5 * v1 + v4c/2;
  vans.set (27, 29, 51);
  if (v4.isNear(vans)) {
    cout << "OK: " << v3 << "+ 5 *" << v1 <<
			"+" << v4c << "/ 2 =" <<  v4 << endl;
  } else {
    cout << "?????" << v3 << "+ 5 *" << v1 <<
			"+" << v4c << "/ 2 =" <<  v4 << endl;
    problems = true;
  }

	cout << "{ Modify/assign on whole vectors }\n";

  v1.set ( 4, 8, 6 );
  cout << "    ... v1 =" << v1 << endl;
  v1 *= .5;
  vans.set ( 2, 4, 3 );
  if (v1.isNear(vans)) {
    cout << "OK: v1 *= .5 ==>" << v1 << endl;
  } else {
    cout << "????? v1 *= .5 ==>" << v1 << endl;
    problems = true;
  }

  v1 /= .1;
  vans.set ( 20, 40, 30 );
  if (v1.isNear(vans)) {
    cout << "OK: v1 /= .1 ==>" << v1 << endl;
  } else {
    cout << "????? v1 /= .1 ==>" << v1 << endl;
    problems = true;
  }

  v1 += SpaceVector ( 5, -5, 10 );
  vans.set ( 25, 35, 40 );
  if (v1.isNear(vans)) {
    cout << "OK: v1 += (5, -5, 10) ==>" << v1 << endl;
  } else {
    cout << "????? v1 += (5, -5, 10) ==>" << v1 << endl;
    problems = true;
  }

  v1 -= v1/5;
  vans.set ( 20, 28, 32 );
  if (v1.isNear(vans)) {
    cout << "OK: v1 -= v1/5 ==>" << v1 << endl;
  } else {
    cout << "????? v1 -= v1/5 ==>" << v1 << endl;
    problems = true;
  }

}  // arithmetic()



void comparisons()  {

// Comparisons

  cout << "\n --- Comparisons\n\n";

  v1.set (2, 4,  6);
  v2.set (1, 5, 9);
  // Previously:
  //  const SpaceVector v3c( 3, 5, 7);
  //  const SpaceVector v4c( 6, 2, 8);
  v5 = v1;

	cout << "{ == and != }\n";

  if ( ! (v1 == v5) ) {
    cout << "????? v1 == v5 false  erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " == " << v5 << endl;

  if ( v1 == v2 ) {
    cout << "????? v1 == v2 true erroneously! v1 = "
	<< v1 << " v2 = " << v2 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " == " << v2 << endl;

  if ( v1 != v5 ) {
    cout << "????? v1 != v5 true erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " != " << v5 << endl;

  if ( ! (v1 != v2) ) {
    cout << "????? v1 != v2 false erroneously! v1 = "
	<< v1 << " v2 = " << v2 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " != " << v2 << endl;

  if ( v1 == v3c ) {
    cout << "????? v1 == v3c true erroneously! v1 = "
	<< v1 << " v3c = " << v3c << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " == " << v2 << endl;

  if ( ! (v3c != v4c) ) {
    cout << "????? v3c != v4c false erroneously! v3c = "
	<< v3c << " v4c = " << v4c << endl;
    problems = true;
  }
  cout << "OK: " << v3c << " != " << v4c << endl;


	cout << "{ Dictionary Ordering: >, <, >=, <= }\n";

// Modification 2/6/98:  Changed a,b,c to c,b,a to reflect change in ordering.

  v1.set (6, 4, 2);			// CHANGED 2/6/98 from (2,4,6)
  v2.set (9, 5, 1);			// CHANGED 2/6/98 from (1,5,9)
  // Previously:
  //  const SpaceVector v3c( 3, 5, 7);
  const SpaceVector v33c (7,5,3);	// CHANGED 2/6/98
  v5.set ( 3, 5, 2);			// CHANGED 2/6/98 from (2,5,3)
  v6 = v1;

  if ( ! (v1 <  v5) ) {
    cout << "????? v1 < v5 false erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <  " << v5 << endl;

  if ( v1 <  v6 ) {
    cout << "????? v1 < v6 true erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " <  " << v6 << endl;

  if ( v1 >  v6 ) {
    cout << "????? v1 > v6 true erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " >  " << v6 << endl;

  if ( v5 >=  v33c ) {				// CHANGED 2/6/98
    cout << "????? v5 >= v33c true erroneously! v5 = "
	<< v5 << " v5 = " << v33c << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v5 << " >=  " << v33c << endl;

   if ( ! (v1 <= v5 ) ) {
    cout << "????? v1 <= v5 false erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <= " << v5 << endl;

   if ( ! (v1 <= v6 ) ) {
    cout << "????? v1 <= v6 false erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <= " << v6 << endl;

   if ( ! (v6 <= v1 ) ) {
    cout << "????? v6 <= v1 true erroneously! v6 = "
	<< v6 << " v1 = " << v1 << endl;
    problems = true;
  }
  cout << "OK: " << v6 << " <= " << v1 << endl;

}  // comparisons()



void tolerance()  {

// Comparisons within tolerance

  cout << "\n --- IsNear(): Comparisons within tolerance\n\n";


	cout << "{ dot and mag2() as used by isNear }\n";

  v1.set( 3, 5, 1);
  v2.set( 4, 8, 7);
  xx = v1.dot(v2);
  if (xx==59) {
    cout << "OK: " << v1 << ".dot" << v2 << "= " << xx << endl;
  } else {
    cout << "?????" << v1 << ".dot" << v2 << "= " << xx << endl;
    problems = true;
  }

  xx = v1.mag2();
  if (xx==35) {
    cout << "OK: " << v1 << ".mag2() = " << xx << endl;
  } else {
    cout << "?????" << v1 << ".mag2() = " << xx << endl;
    problems = true;
  }
	cout << " { isNear() with default tolerance }\n";

  double toler = SpaceVector::setTolerance(1.);
  cout << "Default tolerance starts at (expected 2.22e-14) ==> "
							<< toler << endl;
  SpaceVector::setTolerance(toler);

  v3 = v1 + .2 * toler * v2;
  if (v3 != v1) {
    cout << "OK: " << "adding .2 * tolerance * times " << v2 << "to v1 =" << v1
			<< "gives V3 != V1 \n";
  } else {
    cout << "????? adding .2 * tolerance * times " << v2 << "to v1 =" << v1
			<< "gives V3 == V1 \n";
    problems = true;
  }

  if ( v3 .isNear(v1) ) {
    cout << "OK: " << "Added .2 * tolerance * v2 -- isNear\n";
  } else {
    cout << "????? Added .2 * tolerance -- isNear FALSE\n";
    problems = true;
  }

  v3 = v1 + .6 * toler * v2;
  if ( v1 .isNear(v3) ) {
    cout << "????? " << "Added .6 * tolerance * v2 -- isNear\n";
    problems = true;
  } else {
    cout << "OK: Added .6 * tolerance * -- isNear FALSE\n";
  }

  v3 = v1 + 1.2 * toler * v2;
  if (  v3.isNear(v1) ) {
    cout << "Added 1.2 * tolerance * v2 -- isNear TRUE\n";
  } else {
    cout << "Added 1.2 * tolerance * v2 -- isNear FALSE\n";
  }

  v3 = v1 + 1.8 * toler * v2;
  if (  v3.isNear(v1) ) {
    cout << "????? Added 1.8 * tolerance * v2 -- isNear TRUE\n";
    problems = true;
  } else {
    cout << "OK: " << "Added 1.2 * tolerance * v2 -- isNear FALSE\n";
  }

	cout << " { isNear() with specified tolerance }\n";


  v1.set ( 1, 3, 6 );
  double m = sqrt(v1.mag2());
  cout << "    .... dealing with a vector" << v1
				<< "of magnitude " << m << endl;
  double small = .01;
  v2 = v1 + small * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "with (default) tolerance " << toler << ", "
	<< v1 << "isNear" << v2 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  } else {
    cout << "OK: " << "with (default) tolerance " << toler << ",\n"
	"        " << v1 << "is not Near" << v2 << endl;
  }

  if ( v2.isNear(v1) ) {
    cout << "with tolerance " << toler << ", "  << v2 << "isNear" << v1 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  } else {
    cout << "OK: " << "with tolerance " << toler << ", "
	<< v2 << "is not Near" << v1 << endl;
  }

  if ( v1.isNear(v2 ,  .5*small/m) ) {
    cout << "with tolerance " << .5*small/m <<
	", " << v1 << "isNear" << v2 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  } else {
    cout << "OK: " << "with tolerance " << .5*small/m << ", "
	<< v1 << "is not Near" << v2 << endl;
  }

  if ( v2.isNear(v1 , small/m) ) {
    cout << "with tolerance " << small/m << ", "
	<< v2 << "isNear" << v1 << endl;
  } else {
    cout << "with tolerance " << small/m << ", "
    	<< v2 << "is not Near" << v1 << endl;
  }

  if ( v1.isNear(v2 , 1.5*small/m) ) {
    cout << "with tolerance " << 1.5*small/m << ", "
	<< v1 << "isNear" << v2 << endl;
  } else {
    cout << "with tolerance " << 1.5*small/m << ", "
    	<< v1 << "is not Near" << v2 << endl;
  }

  if ( v1.isNear(v2 , 2*small/m) ) {
    cout << "with tolerance " << 2*small/m << ", "
	<< v1 << "isNear" << v2 << endl;
  } else {
    cout << "with tolerance " << 2*small/m << ", "
    	<< v1 << "is not Near" << v2 << endl;
  }


  if ( v1.isNear(v2 , 3.5*small/m) ) {
    cout << "OK: " << "with tolerance " << 3.5*small/m << ", "
	<< v1 << "isNear" << v2 << endl;
  } else {
    cout << "with tolerance " << 3.5*small/m << ", "
    	<< v1 << "is not Near" << v2 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  }

	cout << " { setTolerance }\n";

  SpaceVector::setTolerance ( .1 );
  cout << "tolerance set to .1\n";

  v2 = v1 + .02 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "OK: " <<
		"with new default tolerance " << v1 << "isNear" << v2 << endl;
  } else {
    cout << "with new default tolerance " << v1 << "is not Near" << v2 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  }

  v2 = v1 + .05 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "OK: " <<
	"with new default tolerance " << v1 << "isNear" << v2 << endl;
  } else {
    cout << "with new default tolerance " << v1 << "is not Near" << v2 << endl;
    cout << "????? That should not be!\n";
    problems = true;
  }

  v2 = v1 + .1 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 << "isNear" << v2 << endl;
  } else {
    cout << "????? with new default tolerance " <<
			v1 << "is not Near" << v2 << endl;
    problems = true;
  }

  v2 = v1 + .2 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 << "isNear" << v2 << endl;
  } else {
    cout << "?????with new default tolerance " <<
	v1 << "is not Near" << v2 << endl;
    problems = true;
  }

  v2 = v1 + .3 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 << "isNear" << v2 << endl;
  } else {
    cout << "????? with new default tolerance " <<
			v1 << "is not Near" << v2 << endl;
    problems = true;
  }


  v2 = v1 + .45 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "?????with new default tolerance " << v1 << "isNear" << v2 << endl;
    problems = true;
  } else {
    cout << "OK: with new default tolerance " << v1 << "is not Near"
		<< v2 << endl;
  }

  v2 = v1 + .6 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "?????with new default tolerance " << v1 << "isNear" << v2 << endl;
    problems = true;
  } else {
    cout << "OK: with new default tolerance " << v1 << "is not Near" <<
		v2 << endl;
  }

  v2 = v1 + 1.2 * SpaceVector( 1, -1, 1);

  if ( v1.isNear(v2) ) {
    cout << "????? with new default tolerance " << v1 << "isNear" << v2 << endl;
    problems = true;
  } else {
    cout << "OK: " <<
	"with new default tolerance " << v1 << "is not Near" << v2 << endl;
  }

}  // tolerance()


void nearness() {

  cout << "\n --- Nearness methods \n\n";

  SpaceVector v1 ( 1, 2, 3 );
  SpaceVector v2 ( 1, 2, 5 );

  if ( close (v1.howNear(v2), 1./sqrt(5.)) ) {
    cout << "OK: " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
  } else {
    cout << "????? " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
   problems = true;
  }

  v1.setZ(-1.1);
  if ( v1.howNear(v2) == 1 ) {
    cout << "OK: " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
  } else {
    cout << "????? " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
   problems = true;
  }

  v1.set  ( 1, 2, -1 );
  v2.set  ( 1, 2, 5 );
  if ( v1.howNear(v2) == 1 ) {
    cout << "OK: " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
  } else {
    cout << "????? " << v1 << ".howNear" << v2 << "= " << v1.howNear(v2) << endl;
   problems = true;
  }

  if ( v1.howParallel(v2) == 1 ) {
    cout << "OK: " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
  } else {
    cout << "????? " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
   problems = true;
  }

  if ( v1.howOrthogonal(v2) == 0 ) {
    cout << "OK: " << v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
  } else {
    cout << "????? " << v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
   problems = true;
  }

  v1.set  ( 3, 4, -5 );
  v2.set  ( 4, 3, 2  );
  if ( v1.howParallel(v2) == 1 ) {
    cout << "OK: " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
  } else {
    cout << "????? << " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
   problems = true;
  }

  if ( close(v1.howOrthogonal(v2), 14/sqrt(1254.)) ) {
    cout << "OK: " << v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
  } else {
    cout << "????? " <<  v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
   problems = true;
  }

  v1.set  ( 6, 7, 8   );
  v2.set  ( 8, 9, 10  );
  if ( close(v1.howParallel(v2),sqrt(24.)/191) ) {
    cout << "OK: " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
  } else {
    cout << "????? << " << v1 << ".howParallel" << v2 << "= "
		<< v1.howParallel(v2) << endl;
   problems = true;
  }

  if ( v1.howOrthogonal(v2) == 1 ) {
    cout << "OK: << " << v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
  } else {
    cout << "????? << " << v1 << ".howOrthogonal" << v2 << "= "
		<< v1.howOrthogonal(v2) << endl;
   problems = true;
  }

  SpaceVector v3 ( 10, 1.2, ETA, .25, RADIANS );
  SpaceVector v4 ( 25, 0.9, ETA, -.15, RADIANS );

  if ( close(v3.deltaR(v4), .5) ) {
    cout << "OK: " << v3 << "\n    .deltaR" << v4 << "= "
		<< v3.deltaR(v4) << endl;
  } else {
    cout << "????? " << v3 << "\n    .deltaR" << v4 << "= "
		<< v3.deltaR(v4) << endl;
   problems = true;
  }

} // nearness()





void combining()  {

// Combining vectors

  cout << "\n --- Combining two vectors\n\n";

  v1.set(2, 4, 5);
  v2.set(6, 1, -4);

	cout << " { dot product }\n";

  double vdv = v1.dot(v2);
  if (vdv == -4) {
    cout << "OK: " << v1 << ".dot" << v2 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v1 << ".dot" << v2 << " ==> " << vdv << endl;
    problems = true;
  }

  vdv = v2.dot(v2);
  if (vdv == 53) {
    cout << "OK: " << v2 << ".dot" << v2 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v2 << ".dot" << v2 << " ==> " << vdv << endl;
    problems = true;
  }

	cout << " { cross product }\n";

  v3 = v1.cross(v2);
  if ( v3 == SpaceVector ( -21, 38, -22 ) ) {
    cout << "OK: " << v1 << ".cross" << v2 << " ==> " << v3 << endl;
  } else {
    cout << "????? " << v1 << ".cross" << v2 << " ==> " << v3 << endl;
    problems = true;
  }

  if ( (v3 + v2.cross(v1)) == SpaceVector (0, 0, 0) ) {
    cout << "OK: " << v1 << ".cross" << v2 << " + " <<
	v2 << ".cross" << v1 << " ==> (0, 0, 0)" << endl;
  } else {
    cout << "?????: " << v1 << ".cross" << v2 << " + " <<
	v2 << ".cross" << v1 << " ==> " << (v3 + v2.cross(v1)) << endl;
    problems = true;
  }

  if ( (v1.cross(v1)) == SpaceVector (0, 0, 0) ) {
    cout << "OK: " << v1 << ".cross" << v1 << " ==> (0, 0, 0) " << endl;
  } else {
    cout << "?????: " << v1 << ".cross" << v1 << " ==> "
	<< (v1.cross(v1)) << endl;
    problems = true;
  }


	cout << " { diff2 }\n";

  v1.set(2, 4, 5);
  v2.set(6, 1, -4);

  vdv = v1.diff2(v2);
  if (vdv == 106) {
    cout << "OK: " << v1 << ".diff2" << v2 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v1 << ".diff2" << v2 << " ==> " << vdv << endl;
    problems = true;
  }


	cout << " { isParallel }\n";

  double percent = .01;
  SpaceVector::setTolerance(percent);
  if ( SpaceVector::getTolerance() == percent ) {
    cout << "OK: new default tolerance is "
		<< SpaceVector::getTolerance() << endl;
  } else {
    cout << "????? new default tolerance is "
		<< SpaceVector::getTolerance() << endl;
    problems = true;
  }

  v1.set(2, 4, 5);
  v2.set(4, 8, 10);
  v5.set(-10, 4, 8);

  if ( v1.isParallel(v2,0) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << "exactly" << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel" << v2 << " FALSE" << endl;
    problems = true;
  }

  v4 = v2;
  v2 = v4 + .3 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + "
				<< (.3*percent*v5)  << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + .6 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + " <<
				(.6*percent*v5) << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + .9 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + " <<
			(.9*percent*v5) << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + 1.0 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + " <<
			(1.0*percent*v5) << " FALSE" << endl;
    problems = true;
  }


  // Note that we would expect to find the vectors non-parallel when
  // the multiplier used here is about 1.09 times the tolerance.

  v2 = v4 + 1.1 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + "
				<< (1.1*percent*v5) << 	endl;
    problems = true;
  } else {
    cout << "OK: NOT " << v1 << ".isParallel" << v2 << endl;
  }

  v2 = v4 + 1.2 * percent * v5;
  if ( v1.isParallel(v2) ) {
    cout << "?????: " << v1 << ".isParallel\n      " << v4 << " + "
				<< (1.2*percent*v5) << 	endl;
    problems = true;
  } else {
    cout << "OK: NOT " << v1 << ".isParallel" << v2 << endl;
  }

	cout << " { isParallel for huge vectors }\n";

  double multiplier = pow(2.0,400);
  v1 *= multiplier;

  v2 = v4 + 0.4 * percent * v5;
  v2 *= multiplier;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel" << v2 << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + 1.0 * percent * v5;
  v2 *= multiplier;
  if ( v1.isParallel(v2) ) {
    cout << "OK: " << v1 << ".isParallel" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isParallel" << v2 << " FALSE" << endl;
    problems = true;
  }

  // Note that we would expect to find the vectors non-parallel when
  // the multiplier used here is about 1.09 times the tolerance.

	// THIS TEST CURRENTLY FAILS!!!!!!!!!!

  v2 = v4 + 1.1 * percent * v5;
  v2 *= multiplier;
  if ( v1.isParallel(v2) ) {
    cout << "?????: " << v1 << ".isParallel" << v2 << endl;
    problems = true;
  } else {
    cout << "OK: NOT " << v1 << ".isParallel" << v2 << endl;
  }


	cout << " { isOrthogonal }\n";

  percent = .01;
  SpaceVector::setTolerance(percent);
  cout << "    ... now default tolerance is "
		<< SpaceVector::getTolerance() << endl;

  v1.set(6, -3, 5);
  v4.set(2, 14, 6);
  v5.set(10, -5, 10);

  // v5 is about the size of the others, and is nearly aligned with v1.
  // It turns out that when you add about 1.03 epsilon*v5 to v4, you barely
  // miss being orthogonal to v1 by our definition.

  v2 = v4;
  if ( v1.isOrthogonal(v2,0) ) {
    cout << "OK: " << v1 << ".isOrthogonal" << v2 << "exactly" << endl;
  } else {
    cout << "?????: " << v1 << ".isOrthogonal" << v2 << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + .6 * percent * v5;
  if ( v1.isOrthogonal(v2) ) {
    cout << "OK: " << v1 << ".isOrthogonal" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isOrthogonal\n      " << v4 << " + " <<
				(.6*percent*v5) << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + .9 * percent * v5;
  if ( v1.isOrthogonal(v2) ) {
    cout << "OK: " << v1 << ".isOrthogonal" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isOrthogonal\n      " << v4 << " + " <<
			(.9*percent*v5) << " FALSE" << endl;
    problems = true;
  }

  v2 = v4 + 1.0 * percent * v5;
  if ( v1.isOrthogonal(v2) ) {
    cout << "OK: " << v1 << ".isOrthogonal" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isOrthogonal\n      " << v4 << " + " <<
			(1.0*percent*v5) << " FALSE" << endl;
    problems = true;
  }


  // Note that we would expect to find the vectors non-Orthogonal when
  // the multiplier used here is about 1.09 times the tolerance.

  v2 = v4 + 1.1 * percent * v5;
  if ( v1.isOrthogonal(v2) ) {
    cout << "?????: " << v1 << ".isOrthogonal\n      " << v4 << " + "
				<< (1.1*percent*v5) << 	endl;
    problems = true;
  } else {
    cout << "OK: NOT " << v1 << ".isOrthogonal" << v2 << endl;
  }

	cout << " { isOrthogonal for huge vectors }\n";

  multiplier = pow(2.0,400);
  v1 *= multiplier;

  v2 = v4 + 1.0 * percent * v5;
  v2 *= multiplier;
  if ( v1.isOrthogonal(v2) ) {
    cout << "OK: " << v1 << ".isOrthogonal" << v2 << endl;
  } else {
    cout << "?????: " << v1 << ".isOrthogonal" << v2 << " FALSE" << endl;
    problems = true;
  }

  // Note that we would expect to find the vectors non-Orthogonal when
  // the multiplier used here is about 1.03 times the tolerance.

	// THIS TEST CURRENTLY FAILS!!!!!!!!!!

  v2 = v4 + 1.1 * percent * v5;
  v2 *= multiplier;
  if ( v1.isOrthogonal(v2) ) {
    cout << "?????: " << v1 << ".isOrthogonal" << v2 << endl;
    problems = true;
  } else {
    cout << "OK: NOT " << v1 << ".isOrthogonal" << v2 << endl;
  }

}  // combining()



void intrinsic()  {

// Simple intrinsic properties

  cout << "\n --- Simple intrinsic properties\n\n";

	cout << " { mag() }\n";
  v1.set (3, -4, -12);
  if (v1.mag() == 13) {
    cout << "OK: " << v1 << ".mag() = " << v1.mag() << endl;
  } else {
    cout << "?????" << v1 << ".mag() = " << v1.mag() << endl;
    problems = true;
  }


// Physics intrinsic properties

  cout << "\n --- Physics intrinsic properties\n\n";

	cout << " { beta() }\n";
  v1.set (15, 12, -16);
  v1 /= 100;
  if ( close(v1.beta(),0.25) ) {
    cout << "OK: " << v1 << ".beta() = " << v1.beta() << endl;
  } else {
    cout << "?????" << v1 << ".beta() = " << v1.beta() << endl;
    problems = true;
  }

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
// Test beta on unit and longer vectors.
#endif // NOTYET_WAIT_FOR_EXCEPTIONS

	cout << " { gamma() }\n";
  if ( close(v1.gamma(), sqrt(16./15.)) ) {
    cout << "OK: " << v1 << ".gamma() = " << v1.gamma() << endl;
  } else {
    cout << "?????" << v1 << ".gamma() = " << v1.gamma() << endl;
    problems = true;
  }

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
// Test gamma on unit and longer vectors.
#endif // NOTYET_WAIT_FOR_EXCEPTIONS

	cout << " { rapidity() }\n";
  cout << "    ... co-linear rapidity for " << v1 << " ==> "
		<< v1.coLinearRapidity() << endl;

  v1.set(0, 0, .5);
  if ( v1.rapidity() == v1.coLinearRapidity() ) {
    cout << "OK: " << v1 << ".rapidity() = " << v1.rapidity() << endl;
  } else {
    cout << "?????" << v1 << ".rapidity() = " << v1.rapidity() <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set(0, 0, .7);
  if ( v1.rapidity(Z_HAT) == v1.coLinearRapidity() ) {
    cout << "OK: " << v1 << ".rapidity(Z_HAT) = " << v1.rapidity(Z_HAT) << endl;
  } else {
    cout << "????" << v1 << ".rapidity(Z_HAT) = " << v1.rapidity(Z_HAT) <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set(.4, 0, 0);
  if ( v1.rapidity(X_HAT) == v1.coLinearRapidity() ) {
    cout << "OK: " << v1 << ".rapidity(X_HAT) = " << v1.rapidity(X_HAT) << endl;
  } else {
    cout << "?????" << v1 << ".rapidity(X_HAT) = " << v1.rapidity(X_HAT) <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set(0, .3, 0);
  if ( v1.rapidity(Y_HAT) == v1.coLinearRapidity() ) {
    cout << "OK: " << v1 << ".rapidity(Y_HAT) = " << v1.rapidity(Y_HAT) << endl;
  } else {
    cout << "?????" << v1 << ".rapidity(Y_HAT) = " << v1.rapidity(Y_HAT) <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set(0, -.3, 0);
  if ( close(v1.rapidity(Y_HAT), -v1.coLinearRapidity()) ) {
	// Note that == could not be used because log(.7/1.3) differs
	// in the last bit from -log(1.3/.7)!
    cout << "OK: " << v1 << ".rapidity(Y_HAT) = " << v1.rapidity(Y_HAT) << endl;
  } else {
    cout << "?????" << v1 << ".rapidity(Y_HAT) = " << v1.rapidity(Y_HAT) <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set(.3, -.3, 0);
  if ( close(v1.rapidity(X_HAT+Y_HAT), 0) ) {
    cout << "OK: " << v1 << ".rapidity((1,1,0)) = "
	<< v1.rapidity(X_HAT+Y_HAT) << endl;
  } else {
    cout << "?????" << v1 << ".rapidity((1,1,0)) = "
	<< v1.rapidity(X_HAT+Y_HAT) << endl;
    problems = true;
  }

  v1.set(.3, 0, -.3);
  if ( close(v1.rapidity(X_HAT-Z_HAT), v1.coLinearRapidity()) ) {
    cout << "OK: " << v1 << ".rapidity((1,1,0)) = "
	<< v1.rapidity(X_HAT-Z_HAT) << endl;
  } else {
    cout << "?????" << v1 << ".rapidity((1,1,0)) = "
	<< v1.rapidity(X_HAT-Z_HAT) <<
	" coLinear = " << v1.coLinearRapidity() << endl;
    problems = true;
  }

#ifdef NOTYET_WAIT_FOR_EXCEPTIONS
// Test rapidity on unit and longer vectors.
#endif // NOTYET_WAIT_FOR_EXCEPTIONS


}  // intrinsic()



void input()  {

// Input

  cout << "\n --- Input\n\n";

	cout << " { Input of entire vector }\n";

  cin >> v1;
  if (v1 == SpaceVector(99., -99., 95.)) {
    cout << "OK: input vector is " << v1 << endl;
  } else {
    cout << "????? input vector should be (99, -99, 95) but is " << v1 << endl;
    problems=true;
  }

#ifdef ZAP4
  cin >> v1c;
  // Should fail, and does.
  cout << "New value of const vector v1c is " << v1c << endl;
#endif

	cout << " { Input of cartesian components }\n";

  double xtmp, ytmp, ztmp;
  cin >> xtmp >> ytmp >> ztmp;
  v1.setX(xtmp); v1.setY(ytmp); v1.setZ(ztmp);
  if (v1 == SpaceVector(50., 60., 70.)) {
    cout << "OK: input vector is " << v1 << endl;
  } else {
    cout << "????? vector should be (50, 60, 70) but is " << v1 << endl;
    problems=true;
  }

#ifdef ZAP4
  cin >> v1c.x(); // Should fail. Does fail with:
	// svtest.cc:1152: call of overloaded method `operator
	// >>(SpaceVector::WontAlter_X_of_a_const)' is ambiguous
	// candidates are: istream::operator >>(char &)
	// ...
	// operator >>(istream &,SpaceVector::Cin_PutTo_constXPreventer)
#endif

	cout << " { Input of non-cartesian component: r }\n";

  v1.set (30, 40, 120);
  cout << "    ... will input 260 to v1.r() of " << v1 << endl;
  double rtmp;
  cin >> rtmp;
  v1.setR(rtmp);
  if (v1.isNear(SpaceVector(60, 80, 240))) {
    cout << "OK: vector v1 is " << v1 << endl;
  } else {
    cout << "????? vector should be (60, 80, 240)) but is " << v1 << endl;
    problems=true;
  }

}  // input()


SpaceVector u ( .6, .8, 0 );
SpaceVector vu ( 6, 8, 0 );



void direction()  {

// Properties relative to reference direction

  cout << "\n --- Perp and projection relative to direction\n\n";

	cout << " { Using the default ZHAT direction }\n";

  v1.set ( 3, -4, 12 );

  if ( v1.perp2() == 25 ) {
    cout << "OK: " << v1 << ".perp2() = " << v1.perp2() << endl;
  } else {
    cout << "????? " << v1 << ".perp2() = " << v1.perp2() << endl;
    problems=true;
  }

  if ( v1.perp() == 5 ) {
    cout << "OK: " << v1 << ".perp() = " << v1.perp() << endl;
  } else {
    cout << "????? " << v1 << ".perp() = " << v1.perp() << endl;
    problems=true;
  }

  if ( v1.perpPart() == SpaceVector ( 3, -4, 0 ) ) {
    cout << "OK: " << v1 << ".perpPart() = " << v1.perpPart() << endl;
  } else {
    cout << "????? " << v1 << ".perpPart() = " << v1.perpPart() << endl;
    problems=true;
  }

  if ( v1.project() == SpaceVector ( 0, 0, 12 ) ) {
    cout << "OK: " << v1 << ".project() = " << v1.project() << endl;
  } else {
    cout << "????? " << v1 << ".project() = " << v1.project() << endl;
    problems=true;
  }


	cout << " { Using a supplied reference SpaceVector }\n";

  v1.set (2,0,4);

  if ( eq(v1.perp2(vu),18.56) ) {
    cout << "OK: " << v1 << ".perp2" << vu << " = " << v1.perp2(vu) << endl;
  } else {
    cout << "????? " << v1 << ".perp2" << vu << " = " << v1.perp2(vu) << endl;
    problems=true;
  }

  if ( eq(v1.perp(vu),sqrt(18.56)) ) {
    cout << "OK: " << v1 << ".perp" << vu << " = " << v1.perp(vu) << endl;
  } else {
    cout << "????? " << v1 << ".perp" << vu << " = " << v1.perp(vu) << endl;
    problems=true;
  }

  if ( v1.perpPart(vu).isNear(SpaceVector ( 1.28, -.96, 4.0 )) ) {
    cout << "OK: " << v1 << ".perpPart" << vu << " = "
						<< v1.perpPart(vu) << endl;
  } else {
    cout << "????? " << v1 << ".perpPart" << vu << " = "
						<< v1.perpPart(vu) << endl;
    problems=true;
  }

  if ( v1.project(vu).isNear(SpaceVector ( .72, .96, 0 )) ) {
    cout << "OK: " << v1 << ".project" << vu
		<< " = " << v1.project(vu) << endl;
  } else {
    cout << "????? " << v1 << ".project" << vu
		<< " = " << v1.project(vu) << endl;
    problems=true;
  }

}  // direction()


SpaceVector dir ( 4, -3, 12 );  // NOT NORMALIZED!


void anglebetween()  {

// Angle between vector and a direction

  cout << "\n --- Angle between vector and a direction\n\n";

	cout << " { Using the default ZHAT direction }\n";

  v1.set ( 10, .35, RADIANS, -.5, RADIANS );

  double miss;
  double arads = v1.angle();
  miss = arads - .35;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".angle() = " << arads <<
		" = .25 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".angle() = " << arads <<
		" = .25 + " << miss << endl;
    problems=true;
  }

  double csx;
  v1.set ( 10, 60, DEGREES, -PI/2, RADIANS );
  csx = v1.cosTheta();
  miss = csx - .5;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cosTheta() = " << csx <<
		" = .5 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cosTheta() = " << csx <<
		" = .5 + " << miss << endl;
    problems=true;
  }

  csx = v1.cos2Theta();
  miss = csx - .25;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cos2Theta() = " << csx <<
		" = .25 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cos2Theta() = " << csx <<
		" = .25 + " << miss << endl;
    problems=true;
  }



	cout << " { Using a different reference direction }\n";

  v1.set ( 6, 16, 2);  // This is perpendicular to dir.

  arads = v1.angle(dir);
  miss = arads - PI/2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".angle" << dir << " = " << arads <<
		" = PI/2 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".angle" << dir << " = " << arads <<
		" = PI/2 + " << miss << endl;
    problems=true;
  }

  arads = v1.theta(dir);
  miss = arads - PI/2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".theta" << dir << " = " << arads <<
		" = PI/2 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".theta" << dir << " = " << arads <<
		" = PI/2 + " << miss << endl;
    problems=true;
  }

  csx = v1.cosTheta(dir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cosTheta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cosTheta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }

  csx = v1.cos2Theta(dir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cos2Theta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cos2Theta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }

  csx = v1.eta(dir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".eta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".eta" << dir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }


  v1.set ( 13, -13, 26);  	      // This is not perpendicular to dir.
  double rightc2 = 31.*31./(6.*169.); // This is the actual cos2 theta.

  arads = v1.angle(dir);

  miss = arads - acos ( sqrt(rightc2) );
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".angle" << dir << " = " << arads <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".angle" << dir << " = " << arads <<
		" = answer + " << miss << endl;
    problems=true;
  }

  arads = v1.theta(dir);
  if ( eq(arads, v1.angle(dir)) ) {
    cout << "OK: v1.angle(dir) == v1.theta(dir)\n";
  } else {
    cout << "????? v1.angle(dir) != v1.theta(dir)\n";
    problems=true;
  }

  csx = v1.cosTheta(dir);
  miss = csx - sqrt(rightc2);
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cosTheta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cosTheta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

  csx = v1.cos2Theta(dir);
  miss = csx - rightc2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cos2Theta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cos2Theta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

  csx = v1.eta(dir);
  miss = csx + log ( tan (acos(sqrt(rightc2))/2.) );
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".eta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".eta" << dir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

	cout << " { Eta with reference direction parallel }\n";

  v1 = 3 * u;
  if ( v1.eta(u) > 17 ) {
    cout << "OK: " << v1 << ".eta()" << u << " = " << v1.eta(u) << endl;
  } else {
    cout << "?????" << v1 << ".eta()" << u << " = " << v1.eta(u) << endl;
    problems=true;
  }

  v1 = -2 * u;
  if ( v1.eta(u) < -17 ) {
    cout << "OK: " << v1 << ".eta()" << u << " = " << v1.eta(u) << endl;
  } else {
    cout << "?????" << v1 << ".eta()" << u << " = " << v1.eta(u) << endl;
    problems=true;
  }



	cout << " { Angle between two SpaceVectors }\n";

  SpaceVector vdir ( 4, -3, 12 );

  v1.set ( 6, 16, 2);  // This is perpendicular to vdir.

  arads = v1.angle(vdir);
  miss = arads - PI/2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".angle" << vdir << " = " << arads <<
		" = PI/2 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".angle" << vdir << " = " << arads <<
		" = PI/2 + " << miss << endl;
    problems=true;
  }

  arads = v1.theta(vdir);
  miss = arads - PI/2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".theta" << vdir << " = " << arads <<
		" = PI/2 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".theta" << vdir << " = " << arads <<
		" = PI/2 + " << miss << endl;
    problems=true;
  }

  csx = v1.cosTheta(vdir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cosTheta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cosTheta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }

  csx = v1.cos2Theta(vdir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cos2Theta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cos2Theta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }

  csx = v1.eta(vdir);
  miss = csx;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".eta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".eta" << vdir << " = " << csx <<
		" = 0 + " << miss << endl;
    problems=true;
  }


  v1.set ( 13, -13, 26);  	      // This is not perpendicular to vdir.
  rightc2 = 31.*31./(6.*169.); 	      // This is the actual cos2 theta.

  arads = v1.angle(vdir);

  miss = arads - acos ( sqrt(rightc2) );
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".angle" << vdir << " = " << arads <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".angle" << vdir << " = " << arads <<
		" = answer + " << miss << endl;
    problems=true;
  }

  arads = v1.theta(vdir);
  if ( eq(arads, v1.angle(vdir)) ) {
    cout << "OK: v1.angle(vdir) == v1.theta(vdir)\n";
  } else {
    cout << "????? v1.angle(vdir) != v1.theta(vdir)\n";
    problems=true;
  }

  csx = v1.cosTheta(vdir);
  miss = csx - sqrt(rightc2);
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cosTheta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cosTheta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

  csx = v1.cos2Theta(vdir);
  miss = csx - rightc2;
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".cos2Theta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".cos2Theta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

  csx = v1.eta(vdir);
  miss = csx + log ( tan (acos(sqrt(rightc2))/2.) );
  if (fabs(miss) < teeny ) {
    cout << "OK: " << v1 << ".eta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
  } else {
    cout << "????? " << v1 << ".eta" << vdir << " = " << csx <<
		" = answer + " << miss << endl;
    problems=true;
  }

	cout << " { Eta with reference direction parallel to v }\n";

  v1 = 3 * u;
  v2 = 5 * u;
  if ( v1.eta(v2) > 17 ) {
    cout << "OK: " << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl;
  } else {
    cout << "?????" << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl;
    problems=true;
  }

  v1 = -2 * u;
  if ( v1.eta(v2) < -17 ) {
    cout << "OK: " << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl;
  } else {
    cout << "?????" << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl;
    problems=true;
  }

}  // anglebetween()



void angledecomposition()  {

// Decomposition of angle between vectors

  cout << "\n --- Decomposition of angle between vectors\n\n";

     cout << " { Vector and Unit, relative to the default ZHAT direction }\n";

  // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12) (not normalized) 

  v1.set ( 4, -3, -5 );
  cout << "    ... v.theta is " << v1.theta()*180/PI << " degrees; u.theta is "
		<< u.theta()*180/PI << " degrees " << endl;
//if (close (v1.polarAngle(u), -PI/4) ) {
  if (close (v1.polarAngle(u),  PI/4) ) {	     // JMM definition change?
    cout << "OK: polar angle (wrt. Z axis) from v to u is " <<
		v1.polarAngle(u)*180/PI << " degrees" << endl;
  } else {
    cout << "?????polar angle (wrt. Z axis) from v to u is " <<
		v1.polarAngle(u)*180/PI << " degrees" << endl;
    problems=true;
  }

  v1.set ( -8, 6, 5 ) ;
  cout << "    ... v.phi is " << v1.phi()*180/PI << "degrees;" <<
		" u.phi is " << u.phi()*180/PI << " degrees " << endl;
  if (close (v1.azimAngle(u), -PI/2) ) {
    cout << "OK: azim angle (wrt. Z axis) from v to u is " <<
		v1.azimAngle(u)*180/PI << " degrees" << endl;
  } else {
    cout << "?????azim angle (wrt. Z axis) from v to u is " <<
		v1.azimAngle(u)*180/PI << " degrees" << endl;
    problems=true;
  }


    cout << " { Vector and vector, relative to the default ZHAT direction }\n";

  v1.set ( 4, -3, -5 );
  cout << "    ... v.theta is " << v1.theta()*180/PI << " degrees;"
		" vu.theta is "
		<< vu.theta()*180/PI << " degrees " << endl;
  if (close (vu.polarAngle(v1), PI/4) ) {
    cout << "OK: polar angle (wrt. Z axis) from vu to v is " <<
		vu.polarAngle(v1)*180/PI << " degrees" << endl;
  } else {
    cout << "?????polar angle (wrt. Z axis) from vu to v is " <<
		vu.polarAngle(v1)*180/PI << " degrees" << endl;
    problems=true;
  }

  v1.set ( -8, 6, 5 ) ;
  cout << "    ... v.phi is " << v1.phi()*180/PI
		<< " degrees; vu.phi is " << vu.phi()*180/PI
		<<  " degrees " << endl;
  if (close (v1.azimAngle(vu), -PI/2) ) {
    cout << "OK: azim angle (wrt. Z axis) from v to vu is " <<
		v1.azimAngle(vu)*180/PI << " degrees" << endl;
  } else {
    cout << "?????azim angle (wrt. Z axis) from v to vu is " <<
		v1.azimAngle(vu)*180/PI << " degrees" << endl;
    problems=true;
  }

	cout << "\n { V U Relative to a different reference direction }\n";

  // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12)/13.

  v1.set ( .7, .8, .1 );

  cout << "  * * * polarAngle of " << v1 << "to" << u << "\n with respect to"
		<< 3.2*dir << "is " << v1.polarAngle(u, dir) *180/PI <<
		" degrees " << endl;
  cout << "  * * * azimAngle of " << v1 << "to" << u << "\n with respect to"
		<< dir << "is " << v1.azimAngle(u, dir) *180/PI <<
		" degrees " << endl;

  if ( close(v1.polarAngle(u, dir),v1.polarAngle(u, 3.2*dir)) ) {
    cout << "OK: polar angle with respect to SpaceVector reference matches\n";
  } else {
    cout << "????? polar angle with respect to SpaceVector reference matches\n";
    problems=true;
  }

  if ( close(v1.azimAngle(u, dir),v1.azimAngle(u, 3.2*dir)) ) {
    cout << "OK: azim  angle with respect to SpaceVector reference matches\n";
  } else {
    cout << "????? azim angle with respect to SpaceVector reference matches\n";
    problems=true;
  }

  v1.set ( -.6, -.8, 0 ); // -u

  if ( close(v1.polarAngle(u, dir) , 0) ) {
    cout << "OK: polarAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is " <<
		v1.polarAngle(u, dir)* 180/PI <<
		" degrees " << endl;
  } else {
    cout << "????? polarAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is " <<
		v1.polarAngle(u, dir)* 180/PI <<
		" degrees " << endl;
    problems=true;
  }

  cout <<   "    Explanation:  Angle between " << v1 << " and " << dir <<
	" is " << v1.angle(dir) <<
	  "\n                  Angle between " << u << " and " << dir <<
	" is " << u.angle(dir) << "\n";

  if ( close(fabs(v1.azimAngle(u, dir)) , PI ) ) {
    cout << "OK: azimAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is " <<
		v1.azimAngle(u, dir)* 180/PI << " degrees" << endl;
  } else {
    cout << "????? azimAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is " <<
		v1.azimAngle(u, dir)* 180/PI << " degrees" << endl;
    problems=true;
  }

  v1.set ( -.6, -.8, .1 ); // Fairly close to -u
  cout << "    * * * polarAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is "
		<< v1.polarAngle(u, dir) *180/PI << " degrees" << endl;
  cout << "    * * * azimAngle of " << v1 << "to" << u
	<< "\n with respect to" << dir << "is " <<
		v1.azimAngle(u, dir) *180/PI << " degrees" << endl;



	cout << "\n { V V Relative to a different reference direction }\n";

  // Reminder -- vu is (6, 8, 0) and dir is (4, -3, 12)/13.

  v1.set ( .7, .8, .1 );

  cout << "    * * * polarAngle of "
		<< v1 << "to" << vu <<
		"\n         with respect to"
		<< dir << "is " << v1.polarAngle(vu, dir) *180/PI
			<< " degrees" << endl;
  cout << "    * * * azimAngle of " << v1 << "to" << vu <<
		"\n         with respect to"
		<< dir << "is " << v1.azimAngle(vu, dir) *180/PI <<
		" degrees" << endl;

  v1.set ( -.6, -.8, 0 ); // Fairly close to vu

  if ( close(v1.polarAngle(vu, dir) , 0) ) {
    cout << "OK: polarAngle of " << v1 << "to" << vu
	<< "        \n with respect to" << dir << "is " <<
		v1.polarAngle(vu, dir)* 180/PI << " degrees" << endl;
  } else {
    cout << "????? polarAngle of " << v1 << "to" << vu
	<< "        \n with respect to" << dir << "is " <<
		v1.polarAngle(vu, dir)* 180/PI << " degrees" << endl;
    problems=true;
  }

  if ( close(fabs(v1.azimAngle(vu, dir)) , PI ) ) {
    cout << "OK: azimAngle of " << v1 << "to" << vu
	<< "\n with respect to" << dir << "is " <<
		v1.azimAngle(vu, dir)* 180/PI << " degrees" << endl;
  } else {
    cout << "????? azimAngle of " << v1 << "to" << vu
	<< "\n with respect to" << dir << "is " <<
		v1.azimAngle(vu, dir)* 180/PI << " degrees" << endl;
    problems=true;
  }

  v1.set ( -.6, -.8, .1 );
  cout << "    * * * polarAngle of " << v1 << "to" << vu
	<< "\n         with respect to" << dir << "is "
		<< v1.polarAngle(vu, dir) *180/PI << " degrees" << endl;
  cout << "    * * * azimAngle of " << v1 << "to" << vu
	<< "\n         with respect to" << dir << "is "
	<< v1.azimAngle(vu, dir)*180/PI  << " degrees" << endl;

}  // angledecomposition()



void rotations()  {

// ROTATIONS

  cout << "\n --- Rotations\n\n";

	cout << " { 90 degree rotations about X Y or Z axis }\n";

  v1.set ( 9, 5, -3 );

  SpaceVector::setTolerance ( .0001 );

  v2 = v1;
  v2.rotateX( PI/2 );
  v3.set ( 9, 3, 5 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateX(90 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateX(90 degrees) =" << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateX( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateX( -PI/2 )    =" << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateX( -PI/2 ) =" << v2 << endl;
    problems=true;
  }

  v2 = v1;
  v2.rotateY( PI/2 );
  v3.set ( -3, 5, -9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateY(90 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateY(90 degrees) =" << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateY( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateY( -PI/2 )    =" << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateY( -PI/2 ) =" << v2 << endl;
    problems=true;
  }

  v2 = v1;
  v2.rotateZ( PI/2 );
  v3.set ( -5, 9, -3 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateZ(90 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateZ(90 degrees) =" << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateZ( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateZ( -PI/2 )    =" << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateZ( -PI/2 ) =" << v2 << endl;
    problems=true;
  }


	cout << " { Pair of 90 degree rotations to permute x y z }\n";

  v1.set ( 9, 5, -3 );
  v3.set ( -3, 9, 5 );
  v2 = v1;
  v2.rotateX ( PI/2 );
  v2.rotateZ ( PI/2 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "becomes " << v2 << endl;
  } else {
    cout << "????? " << v1 << "becomes " << v2 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotateZ ( -PI/2 );
  v2.rotateX ( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "becomes " << v2 << endl;
  } else {
    cout << "????? " << v3 << "becomes " << v2 << endl;
    problems=true;
  }


	cout << " { Axis rotations by general angles }\n";

  v1.set ( 1, 10*sqrt(3.), 10 );
  v3.set ( 1, 0, 20 );
  v2 = v1;
  v2.rotateX ( PI/3 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateX(60 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateX(60 degrees) =" << v2 << endl;
    problems=true;
  }
  v4 = rotationXOf (v2, 2*PI/3 );
  v2.rotateX ( 2*PI/3 );
  if (v4.isNear(v2)) {
    cout << "OK: vprime = rotationXOf(v,phi) matches v.rotateX(phi)\n";
  } else {
    cout << "????? vprime = rotationXOf(v,phi) mismatches v.rotateX(phi)\n";
    problems=true;
  }
  v1.setX(v1.x()*-1);
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotate another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
				<< v2 << endl;
    problems=true;
  }

	// THIS TEST CURRENTLY FAILS!!!!!!!!!!

  v1.set ( 10, 1, 10*sqrt(3.) );
  v3.set ( 20, 1, 0 );
  v2 = v1;
  v4 =  rotationYOf (v2, PI/3 );
  v2.rotateY ( PI/3 );
  if (v4.isNear(v2)) {
    cout << "OK: vprime = rotationYOf(v,phi) matches v.rotateY(phi)\n";
  } else {
    cout << "????? vprime = rotationYOf(v,phi) mismatches v.rotateY(phi)\n";
    problems=true;
  }
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateY(60 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateY(60 degrees) =" << v2 << endl;
    problems=true;
  }
  v2.rotateY ( 2*PI/3 );
  v1.setY(v1.y()*-1);
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotate another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
			<< v2 << endl;
    problems=true;
  }

	// THIS TEST CURRENTLY FAILS!!!!!!!!!!

  v1.set ( 10*sqrt(3.), 10, 1 );
  v3.set ( 0, 20, 1 );
  v2 = v1;
  v4 =  rotationZOf (v2, PI/3 );
  v2.rotateZ ( PI/3 );
  if (v4==v2) {
    cout << "OK: vprime = rotationZOf(v,phi) matches v.rotateZ(phi)\n";
  } else {
    cout << "????? vprime = rotationZOf(v,phi) mismatches v.rotateZ(phi)\n";
    problems=true;
  }
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateZ(60 degrees) =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateZ(60 degrees) =" << v2 << endl;
    problems=true;
  }
  v2.rotateZ ( 2*PI/3 );
  v1.setZ(v1.z()*-1);
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotate another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
			<< v2 << endl;
    problems=true;
  }


	cout << " { Verifying the phi++ rule }\n";

	// THIS TEST CURRENTLY FAILS!!!!!!!!!!

  v1.set ( -13, 2, -11 );
  v2 = v1;
  v2.rotateZ ( 2.4 );
  v2.setPhi(v2.phi()-2.4);
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v1 << "rotateZ (2.4) phi -= 2.4 =" << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateZ (2.4) phi -= 2.4 =" << v2 << endl;
    problems=true;
  }

  v1.set ( 5, 21, -1 );
  v2 = v1;
  v2.rotateX ( -2.9 );
  v2.rotateZ ( -PI/2 );
  v2.rotateX ( -PI/2 );
  v2.setPhi(v2.phi()+2.9);
  v2.rotateX ( PI/2 );
  v2.rotateZ ( PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v1 <<
		"rotateX (-2.9) permute phi += 2.9 permute ="
		<< v2 << endl;
  } else {
    cout << "????? " << v1 <<
		"rotateX (-2.9) permute phi += 2.9 permute ="
		<< v2 << endl;
    problems=true;
  }

  v1.set ( 15, -2, -11 );
  v2 = v1;
  v2.rotateY ( 1.9 );
  v2.rotateX (  PI/2 );
  v2.rotateZ (  PI/2 );
  v2.setPhi(v2.phi()-1.9);
  v2.rotateZ ( -PI/2 );
  v2.rotateX ( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v1 << "rotateY (1.9) permute phi -= 1.9 permute ="
		<< v2 << endl;
  } else {
    cout << "????? " << v1 <<
		"rotateY (1.9) permute phi -= 1.9 permute ="
		<< v2 << endl;
    problems=true;
  }

} // Rotations


void rotations2() {

	cout << " { Generic rotations using XYZ axes }\n";

  v1.set (8, -5, 3);

  SpaceVector xv (4,0,0);
  SpaceVector  xu (1,0,0);

  v3 = v2 = v1;
  v2.rotate ( xu, 50*PI/180 );
  v3.rotateX ( 50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( xv, -.4 );
  v3.rotateX( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -xv, .9 );
  v3.rotateX( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  SpaceVector yv (0,5,0);
  SpaceVector  yu (0,-1,0);
  v3 = v2 = v1;
  v2.rotate ( yu, 50*PI/180 );
  v3.rotateY ( -50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( yv, -.4 );
  v3.rotateY( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -yv, .9 );
  v3.rotateY( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  SpaceVector zv (0,0,6);
  SpaceVector  zu (0,0,1);
  v3 = v2 = v1;
  v2.rotate ( zu, 50*PI/180 );
  v3.rotateZ ( 50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( zv, -.4 );
  v3.rotateZ( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -zv, .9 );
  v3.rotateZ( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }


	cout << " { Generic rotations about arbitrary axes }\n";

  // Laboriously hand-calculated :

  vu.set ( 1, 1, 0 );

  v1.set ( 3, 4, 0 );
  v3 = v2 = v1;
  v2.rotate ( vu, PI/4 );
  cout << "    " << v1 << "rotated 45 degrees about " << vu << "is" <<
		v2 << endl;
  v3.rotate ( vu, -3*PI/4 );
  cout << "    "<< v1 << "rotated -135 degrees about " << vu << "is" <<
		v3 << endl;
  double temp = v3.x(); v3.setX(v3.y()); v3.setY(temp);
  v3.setZ(v3.z()*-1);
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << " ~ " << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  // Using formula described below:

  double phi;

  v1.set(3,4,1);
  v2.set(1,1,0);
  phi = PI/4;
  v3.set( 3.6464, 3.3536, 1.2071 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(3,4,1);
  v2.set(2,2,0);
  phi = PI/2;
  v3.set( 4.2071, 2.7929, 0.7071 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(3,4,1);
  v2.set(1,1,0);
  phi = 175*PI/180;
  v3.set( 4.0597, 2.9403, -0.9346 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

#ifdef INTENTIONAL_ERROR
  // Intentionally wrong to check our checker:
  v1.set(3,4,1);
  v2.set(1,1,0);
  phi = 175*PI/180;
  v3.set( 4.0397, 2.9403, -0.9346 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;
#endif // INTENTIONAL_ERROR

  v1.set(2,-1,3);
  v2.set(4,2,-5);
  phi = 35*PI/180;
  v3.set( 1.5791, -2.7726, 1.9543 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,-1,3);
  v2.set(4,2,-5);
  phi = -10*PI/180;
  v3.set( 1.9316, -0.4214, 3.1767 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,-1,3);
  v2.set(4,2,-5);
  phi = -100*PI/180;
  v3.set( -1.4330, 2.9340, 1.8272 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(0,1,0);
  phi = 20*PI/180;
  v3.set( 3.2475, 3, 3.0747);
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(0,3,0);
  phi = PI/2;
  v3.set(4, 3, -2);
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(0,1,0);
  phi = PI;
  v3.set( -2, 3, -4 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(0,-1,0);
  phi = 20*PI/180;
  v3.set( 0.5113, 3, 4.4428 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(0,-1,0);
  phi = -20*PI/180;
  v3.set( 3.2475, 3, 3.0747);
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(1,3,4);
  v2.set(0,.01,1);
  phi = 1*PI/180;
  v3.set( 0.9482, 3.0170, 3.9998 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(1,3,4);
  v2.set(0,.01,1);
  phi = -10*PI/180;
  v3.set( 1.4988, 2.7814, 4.0022 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(1, -1, 1);
  phi = -130*PI/180;
  v3.set( 3.4531, -2.6866, -3.1397 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(1, -1, 1);
  phi = -80*PI/180;
  v3.set( 5.1537, 0.8318, -1.3220 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(1, -1, 1);
  phi = -PI/6;
  v3.set( 3.8868, 3.0415, 2.1547 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(1, -1, 1);
  phi = 20*PI/180;
  v3.set( 0.5574, 2.3638, 4.8064 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

  v1.set(2,3,4);
  v2.set(1, -1, 1);
  phi = 70*PI/180;
  v3.set( -2.4557, -0.7170, 4.7387 );
  if ( checkAxialRotation (v1, v2, phi, v3 ) ) problems = true;

}  // rotations2()



void rotations3() {
	cout << " { Rotations via AxisAngle argument }\n";

  double phi, theta, psi;

  SpaceVector::setTolerance ( .000000001);

  v1.set(2,3,5);
  v3 = v1;
  v2.set(1.5, -1, 1);
  phi = 75*PI/180;
  v1.rotate(v2, phi);
  v1 = rotationOf (v1, AxisAngle(v2, -phi) );

  if ( v1.isNear(v3) ) {
    cout << "OK: AxisAngle rotation is working" << v3 << "-->" << v1 << endl;;
  } else {
    cout << "?????AxisAngle rotation failed:  " << v3 << "-->" << v1 << endl;
    problems = true;
  }

  v1.set(21,-13,5);
  v3 = v1;
  v2.set(1.5, 2.5, -1);
  phi = 25*PI/180;
  v1.rotate(v2, phi);
  v1.rotate(AxisAngle(v2, -phi) );

  if ( v1.isNear(v3) ) {
    cout << "OK: AxisAngle rotation is working" << v3 << "-->" << v1 << endl;;
  } else {
    cout << "?????AxisAngle rotation failed:  " << v3 << "-->" << v1 << endl;
    problems = true;
  }

	cout << " { Rotations via Euler Angles }\n";

  v1.set (4, -3, -6);
  phi   = 45*PI/180;
  theta = 15*PI/180;
  psi   = -60*PI/180;
  v2 = rotationOf (v1, phi, theta, psi);
  v3 = goldstein107 (v1, phi, theta, psi);
  if ( v2.isNear(v3) ) {
    cout << "OK: Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "=" << v3 << endl;
  } else {
    cout << "????? Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "!=" << v3 << endl;
    problems = true;
  }

  v1.set (-14, 23, 6);
  phi   = -40*PI/180;
  theta = 135*PI/180;
  psi   = 20*PI/180;
  v2 = v1;
  v2.rotate(phi, theta, psi);
  v3 = goldstein107 (v1, phi, theta, psi);
  if ( v2.isNear(v3) ) {
    cout << "OK: Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "=" << v3 << endl;
  } else {
    cout << "????? Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "!=" << v3 << endl;
    problems = true;
  }

  EulerAngles xi;

  v1.set (41, 3, 20);
  phi   = 35*PI/180;
  theta = -75*PI/180;
  psi   =  50*PI/180;
  xi.set(phi, theta, psi);
  v2 = rotationOf (v1, phi, theta, psi);
  v3 = goldstein107 (v1, phi, theta, psi);
  if ( v2.isNear(v3) ) {
    cout << "OK: Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "=" << v3 << endl;
  } else {
    cout << "????? Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "!=" << v3 << endl;
    problems = true;
  }

  v1.set (-1, 2, 3);
  phi   = -90*PI/180;
  theta = 15*PI/180;
  psi   = -20*PI/180;
  xi.set(phi, theta, psi);
  v2 = v1;
  v2.rotate(phi, theta, psi);
  v3 = goldstein107 (v1, phi, theta, psi);
  if ( v2.isNear(v3) ) {
    cout << "OK: Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "=" << v3 << endl;
  } else {
    cout << "????? Euler angle rotation of" << v1 << endl
	<< "    by phi = " << phi << " theta = " << theta << " psi = " << psi
 	<< "\n     " << v2 << "!=" << v3 << endl;
    problems = true;
  }

}  // rotations3()