// lt-3.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 3 for the LorentzTransformation
// 		  and LorentzBoost class:
//
// 		Arithmetic and Comparison
//
// * Linear Arithmetic
// * Inversion
// * Comparing LorentzTransformations

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

	// Get the LorentzTransformation class

#include "ZMutility/ZMenvironment.h"

#include "PhysicsVectors/LorentzTransformation.h"
#include "PhysicsVectors/LorentzVector.h"
#include "PhysicsVectors/Rotation.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 it the afore-mentioned output
	// look a bit more readable:

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

const Float8 PI =  3.1415926535897931;


int main() {

	// (See tutorial example lt-1 for how to construct and access
	//  LorentzBoosts and example lt-2 for how to construct and access
	//  LorentzTransformations.)

  LorentzBoost b1( UnitVector(.6, -.5, .1), .7 );
  LorentzBoost b2( SpaceVector( .08, -.85, .50 ));
  LorentzBoost b3( .8, .5, .1 );
  LorentzBoost b4, b5, b6;

  Rotation r1( PI/3, PI, -PI/2);

  LorentzTransformation lt1( b1, r1 );
  LorentzTransformation lt2( r1, b2 );
  LorentzTransformation lt3( b3, r1 );
  LorentzTransformation lt4, lt5, lt6;

// ***** Linear Arithmetic

	// Here is how to mutliply LorentzTransformations:
	//------------------------------------------------

  cout << "lt1 = " << lt1 << endl;
  cout << "lt2 = " << lt2 << endl;
  cout << "lt3 = " << lt3 << endl << endl;
  cout << "lt1*lt2 = " << lt1*lt2 << endl;		    // <-------------

  lt4 = lt2*lt3;					    // <-------------
  cout << "lt4 = lt2*lt3 = " << lt4 << endl;

	// You can also use the modify/assign operators:

  lt3 *= lt1;						    // <-------------
  cout << "lt3 = " << lt3 << endl;


	// Here is how to multiply LorentzBoosts:
	//---------------------------------------

  cout << "b1 = " << b1 << endl;
  cout << "b2 = " << b2 << endl;
  cout << "b3 = " << b3 << endl << endl;
  cout << "b1*b2 = " << b1*b2 << endl; 			    // <-------------

	// You can NOT assign the product of two LorentzBoosts to
	// another LorentzBoost. The reason is that the product could
	// essentially be a rotation, which can not be represented in
	// a LorentzBoost's symmetric matrix. So, the simplest way to find
	// the LorentzBoost product of two LorentzBoosts is to decompose
	// the resulting LorentzTransformation into a boost.

  (b2*b3).decompose( b4, r1 );				    // <-------------
  cout << "b4 = b2*b3 = " << b4 << endl;
  cout << " and a rotation by " << r1 << endl;

        // You can NOT, for the same reason as above, use modify/assign ops:

        // WRONG:  b3 *= b1;


// ***** Inversion

	// Here is how to find the inverse of a LorentzTransformation:
	//------------------------------------------------------------

  lt5 = inverseOf( lt1 );				    // <-------------

  cout << "lt1 = " << lt1 << endl;
  cout << "The inverse of lt1 = " << lt5 << endl;

	// Here is how to invert a LorentzTransformation in place:
	//--------------------------------------------------------

  cout << "lt4 = " << lt4 << endl;
  lt4.invert();						    // <-------------
  cout << "After inverting, lt4 = " << lt4 << endl;

	// Here is how to find the inverse of a LorentzBoost:
	//---------------------------------------------------

  cout << "b2 = " << b2 << endl;
  cout << "But the inverse of b2 = " << inverseOf( b2 )     // <-------------
       << endl;


	// Here is how to invert a LorentzBoost in place:
	//-----------------------------------------------

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


// ***** Comparing LorentzTransformations

	// Here is how to apply exact comparisons to LorentzTransformations:
	//------------------------------------------------------------------

  lt6 = lt1;

  if ( lt1 == lt6 ) {					    // <-------------
    cout << lt1 << " == " << lt6 << endl;
  }

  lt6.invert();

  if ( lt1 != lt6 ) {					    // <-------------
    cout << lt1 << " != " << lt6 << endl;
  }

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


	// Here is how to determine if two LorentzTransformations
	// are "nearly equal":
	//-------------------------------------------------------

  lt5.set( b1, Rotation( PI, PI/2, -PI/5) );
  lt6.set( b1, Rotation( PI, PI/2, -PI/5 + .0000001) );

  if ( lt5 != lt6 ) {
    cout << lt5 << " != " << lt6 << endl;
  }

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

	// Here is how to apply exact comparisons to LorentzBoosts:
	//---------------------------------------------------------

  b6 = b1;

  if ( b1 == b6 ) {					    // <-------------
    cout << b1 << " ==\n" << b6 << endl;
  }

  b6.invert();

  if ( b1 != b6 ) {					    // <-------------
    cout << b1 << " !=\n" << b6 << endl;
  }

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


        // Here is how to determine if two LorentzBoosts are "nearly equal":
        //------------------------------------------------------------------

  b5.set( .5985, .0001, .7394 );
  b6.set( .5985, .0001001, .7394 );

  if ( b5 != b6 ) {
    cout << b5 << " !=\n" << b6 << endl;
  }

  if ( b5.isNear(b6) ) {                                    // <-------------
    cout << "But, " << b5 << " isNear\n" << b6 << endl;
  }

	// Here is how to control the tolerance:
	//--------------------------------------

  lt5.set( b2, Rotation( PI/2, -PI, PI/6) );
  lt6.set( b2, Rotation( PI/2, -PI, PI/6 + .00001) );
  b5.set( .45, .81, .25 );
  b6.set( .45, .81, .250001 );

  Float8 toler = 1.E-2;

  if ( lt5.isNear(lt6, toler) ) {		   	    // <-------------
    cout << "Within a tolerance of " << toler << "\n"
         << lt5 << " isNear " << lt6 << endl;
  }

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

        // To see what the default tolerance is:

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

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

  if ( ! b5.isNear(b6) ) {
    cout << "Within the same default tolerance, \n"
         << b5 << " is not near\n" << b6 << endl;
  }

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

  LorentzBoost::setTolerance(.01);                          // <-------------

	// Notice that above the method was LorentzTransformation::getTolerance
	// Because the "get" and "set" Tolerance methods are a part of the
	// base class, calling them from any derived class works just the same
	// So, the tolerance for a LorentzTransformation is the same as the
	// tolerance for a LorentzBoost or a LorentzBoostX, etc.

  if ( lt5.isNear(lt6) ) {
    cout << "Within the new default tolerance, "
         << LorentzTransformation::getTolerance() << ", " << endl
         << lt5 << " isNear() " << lt6 << endl;
  }

  if ( b5.isNear(b6) ) {
    cout << "Within the new default tolerance, "
         << LorentzTransformation::getTolerance() << ", " << endl
         << b5 << " isNear()\n" << b6 << endl;
  }

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


return 0;

} /* end of main */