// ro-4.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 4 for the Rotation class:
//
// 		Coordinate Axis Rotations
//
// * Constructing Coordinate Axis Rotations
// * Components of the matrix representation
// * Multiplication and Inversion
// * Comparisons
// * Applying to vectors

// Key lines are marked by 				    // <-------------

	// Get the Rotation class

#include "ZMutility/ZMenvironment.h"
#include "PhysicsVectors/Rotation.h"
#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/UnitVector.h"
#include "PhysicsVectors/LorentzVector.h"
ZM_USING_NAMESPACE( zmpv )
USING( std::cout )
USING( std::endl )

	// These classes themselves pull in FixedTypes, but just to
	// remind ourselves that we are using Float8, not double, for
	// 64-bit floats...

#include "ZMutility/FixedTypes.h"
ZM_USING_NAMESPACE( zmfxt )

	// This example will use iostream output, so get that --
	// in a portable way.

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

	// This example also tries to make the afore-mentioned output
	// a bit more readable:

#include "ZMutility/iomanip"
using std::setw;

	// To demonstrate the container usage, we pick an arbitrary type:

#include <vector>
using std::vector;

const Float8 PI =  3.1415926535897931;


int main() {

// ***** Constructing Coordinate Axis Rotations

	// Here is how to construct a RotationZ from the angle delta;
	//-----------------------------------------------------------

  cout << "Constructing rz1 from the angle delta \n";
  Float8 delta = PI/3;

  RotationZ rz1( delta );				    // <-------------
  cout << "rz1 = " << rz1 << endl;

	// RotationX, RotationY, and RotationZ have the same functionality
	// as a regular Rotation, but their internal implementation is more
	// compact and their constructors simpler.

// ***** Components of the matrix representation

        // Here is how to obtain individual matrix components:
        //----------------------------------------------------

  RotationX rx1( -PI/6 );
  cout << setw(15) << rx1.xx()				    // <-------------
       << setw(15) << rx1.xy()
       << setw(15) << rx1.xz() << endl;

  	// An alternative to outputting each matrix component individually
	// is to output them row by row or column by columm:

  cout << rx1.rowY() << endl << rx1.rowZ() << endl;	    // <-------------


	// Here is how to obtain a 3x3 matrix representing
 	//  a Coordinate Axis Rotation:
 	//------------------------------------------------

  RotationY ry1( .125 );
  ZMpvRep3x3 rep = ry1.rep3x3();			    // <-------------

  cout << setw(15) << rep.xx_
       << setw(15) << rep.xy_
       << setw(15) << rep.xz_ << endl;
  cout << setw(15) << rep.yx_
       << setw(15) << rep.yy_
       << setw(15) << rep.yz_ << endl;
  cout << setw(15) << rep.zx_
       << setw(15) << rep.zy_
       << setw(15) << rep.zz_ << endl;


// ***** Multiplication and Inversion

	// Here is how to multiply Axis Rotations of like form:
	//-----------------------------------------------------

  RotationX rx2 ( -5*PI/6 );
  RotationX rx3, rx4;

  rx3 = rx1*rx2;					    // <-------------
  rx4 = rx2*rx1;

  cout << "rx1 = " << rx1 << endl;
  cout << "rx2 = " << rx2 << endl;

  cout << "rx3 = " << rx3 << endl;
  cout << "rx4 = " << rx4 << endl;

  	// Multiply and assign:

  RotationZ rz2 ( 0.5 );
  rz2 *= rz1;						    // <-------------


	// Here is how to multiply Axis Rotations to give a generic Rotation:
	//-------------------------------------------------------------------

  Rotation r1 = rx1 * ry1;				    // <-------------
  cout << "rx1 = " << rx1 << endl;
  cout << "rx2 = " << rx2 << endl;
  cout << " r1 = " << r1 << endl;

	// Multiply and assign:

  r1 *= rz1;						    // <-------------
  cout << "Now, r1 =" << r1 << endl;


	// Here is how to find the inverse of an Axis Rotation:
	//-----------------------------------------------------

  rx3 = inverseOf( rx2 );				    // <-------------

  cout << "rx2 = " << rx2 << endl;
  cout << "Inverse of rx2 = " << rx3 << endl;

	// Or, you can modify a rotation directly:

  rx3.invert();						    // <-------------
  cout << "Now, rx3 = " << rx3 << endl;


// ***** Comparisons

	// Here is how to apply exact comparisons to Axis Rotations:
	//----------------------------------------------------------

  rx3 = rx2;
  if ( rx2 == rx3 ) {					    // <-------------
    cout << rx2 << " == " << rx3 << endl;
  }

	// Other comparisons > < >= <= are also provided so that
 	// sort templates can be used on lists of Rotations.

	// Here is how to determine if two Axis Rotations are "nearly equal"
	//------------------------------------------------------------------

  rx1.set( 2.0000000 );
  rx2.set( 2.0000001 );

  if ( rx1 != rx2 ) {
    cout << rx1 << " != " << rx2 << endl;
  }

  if ( rx1.isNear( rx2 ) ) {				    // <-------------
    cout << "But, " << rx1 << ".isNear() " << rx2 << endl;
  }

	// To see a measure of the "near-equality" of the two rotations:

  Float8 proximity = rx1.howNear( rx2 );		    // <-------------
  cout << "rx1.howNear(rx2) = " << proximity << endl;

	// Here is how to determine if an Axis Rotation
	// is near a generic Rotation
	//---------------------------------------------

  r1.set( UnitVector(1, 0, 0), PI);
  rx1.set( PI - .0005 );

  if ( r1.isNear( rx1 ) ) {				    // <-------------
    cout << r1.getAxis() << ", " << r1.getDelta() << " isNear() "
         << rx1 << endl;
  }


        // Here is how to control the tolerance used
	//  to determine near equality
        //------------------------------------------

  rx1.set( PI/2 );
  rx2.set( PI/2 + .01 );

  Float8 toler = 1.E-2;

  if ( rx1.isNear(rx2, toler) ) {			    // <-------------
    cout << "With a tolerance of " << toler << ",\n"
         << rx1 << " isNear() " << rx2 << endl;
  }

	// To see what the default tolerance is:

  toler = Rotation::getTolerance();			     // <------------

 if ( ! rx1.isNear(rx2) ) {
    cout << "Within the default tolerance of " << toler << ",\n"
         << rx1 << " is not near " << rx2 << endl;
  }

	// The default tolerance is 1E-6, but you can specify a
	// different tolerance as above. Or, you can change the default:

  Rotation::setTolerance(.01);				     // <------------
  if ( rx1.isNear( rx2 ) ) {
    cout << "Within the new default tolerance, "
         << Rotation::getTolerance() << ", " << endl
         << rx1 << " isNear() " << rx2 << endl;
  }

  Rotation::setTolerance( toler );
	// Leaving things as we found them.


// ***** Applying to vectors

	// Here is how to apply a RotationX to a vector to form a new vector:
	//-------------------------------------------------------------------

  SpaceVector sv1( 5, 2, -9);

  SpaceVector sv2 = rx1( sv1 );				     // <------------
  cout << "sv1 = " << sv1 << endl;
  cout << "rotate that by " << rx1 << " to get " << sv2 << endl;

	// Here is how to apply a RotationY to a
	// 4-vector to rotate the space part:
	//--------------------------------------

  LorentzVector w1( -2, 1, 13, Tcomponent(6) );

  cout << "w1 = " << w1 << endl;
  w1 = ry1(w1);						     // <------------
  cout << "rotate that by " << ry1 << " to get w1 = " << w1 << endl;

 	// Here is how to apply a RotationZ to
	// each element of a vector of vectors:
	//------------------------------------

  LorentzVector w2( sv1, 12);
  LorentzVector w3 = w1+w2;
  LorentzVector w4 = 2*w1 + w2 - 5*w3;

#ifdef DEFECT_STANDARD_CPLUSPLUS        // Prevent any problems with containers
#ifndef NT_MSVCPP

#ifdef NEVER
  cout << "Rotating a vector of LorentzVectors by " << rz1 << endl;

  typedef vector<LorentzVector> WV;
  WV  wVec;
      wVec.push_back(w1); wVec.push_back(w2);
      wVec.push_back(w3); wVec.push_back(w4);

  WV wVecRot = rz1( wVec );				     // <------------

  WV::const_iterator w_ptr;
  for ( w_ptr = wVecRot.begin(); w_ptr != wVecRot.end(); ++w_ptr ) {
    cout << *w_ptr << endl;
  }

#endif
#endif // NT_MSVCPP
#endif // DEFECT_STANDARD_CPLUSPLUS

  return 0;

} /* end of main */