00001 // @(#)root/mathcore:$Id: Quaternion.cxx 22516 2008-03-07 15:14:26Z moneta $ 00002 // Authors: W. Brown, M. Fischler, L. Moneta 2005 00003 00004 /********************************************************************** 00005 * * 00006 * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team * 00007 * * 00008 * * 00009 **********************************************************************/ 00010 00011 // Implementation file for rotation in 3 dimensions, represented by quaternion 00012 // 00013 // Created by: Mark Fischler Thurs June 9 2005 00014 // 00015 // Last update: $Id: Quaternion.cxx 22516 2008-03-07 15:14:26Z moneta $ 00016 // 00017 #include "Math/GenVector/Quaternion.h" 00018 00019 #include <cmath> 00020 00021 #include "Math/GenVector/Cartesian3D.h" 00022 #include "Math/GenVector/DisplacementVector3D.h" 00023 #include "Math/GenVector/Quaternion.h" 00024 00025 #include "Math/GenVector/Rotation3Dfwd.h" 00026 #include "Math/GenVector/AxisAnglefwd.h" 00027 #include "Math/GenVector/EulerAnglesfwd.h" 00028 00029 namespace ROOT { 00030 00031 namespace Math { 00032 00033 // ========== Constructors and Assignment ===================== 00034 00035 void Quaternion::Rectify() 00036 { 00037 00038 // The vector should be a unit vector, and the first element should be 00039 // non-negative (though negative fU quaternions would work just fine, 00040 // being isomorphic to a quaternion with positive fU). 00041 00042 if ( fU < 0 ) { 00043 fU = - fU; fI = - fI; fJ = - fJ; fK = - fK; 00044 } 00045 00046 Scalar a = 1.0 / std::sqrt(fU*fU + fI*fI + fJ*fJ + fK*fK); 00047 fU *= a; 00048 fI *= a; 00049 fJ *= a; 00050 fK *= a; 00051 00052 } // Rectify() 00053 00054 00055 // ========== Operations ===================== 00056 00057 // DisplacementVector3D< Cartesian3D<double> > 00058 // Quaternion::operator() (const DisplacementVector3D< Cartesian3D<double> > & v) const 00059 // { 00060 // // apply to a 3D Vector 00061 // } 00062 00063 // Quaternion Quaternion::operator * (const Quaternion & q) const { 00064 // // combination of rotations 00065 // return Quaternion ( 00066 // fU*q.fU - fI*q.fI - fJ*q.fJ - fK*q.fK 00067 // , fU*q.fI + fI*q.fU + fJ*q.fK - fK*q.fJ 00068 // , fU*q.fJ - fI*q.fK + fJ*q.fU + fK*q.fI 00069 // , fU*q.fK + fI*q.fJ - fJ*q.fI + fK*q.fU ); 00070 //} 00071 00072 Quaternion Quaternion::operator * (const Rotation3D & r) const { 00073 // combination of rotations 00074 return operator* ( Quaternion(r) ); 00075 } 00076 00077 Quaternion Quaternion::operator * (const AxisAngle & a) const { 00078 // combination of rotations 00079 return operator* ( Quaternion(a) ); 00080 } 00081 00082 Quaternion Quaternion::operator * (const EulerAngles & e) const { 00083 // combination of rotations 00084 return operator* ( Quaternion(e) ); 00085 } 00086 00087 Quaternion Quaternion::operator * (const RotationZYX & r) const { 00088 // combination of rotations 00089 return operator* ( Quaternion(r) ); 00090 } 00091 00092 Quaternion::Scalar Quaternion::Distance(const Quaternion & q) const { 00093 // distance 00094 Scalar chordLength = std::fabs(fU*q.fU + fI*q.fI + fJ*q.fJ + fK*q.fK); 00095 if (chordLength > 1) chordLength = 1; // in case roundoff fouls us up 00096 return acos(chordLength); 00097 } 00098 00099 // ========== I/O ===================== 00100 00101 std::ostream & operator<< (std::ostream & os, const Quaternion & q) { 00102 // TODO - this will need changing for machine-readable issues 00103 // and even the human readable form may need formatiing improvements 00104 os << "\n{" << q.U() << " " << q.I() 00105 << " " << q.J() << " " << q.K() << "}\n"; 00106 return os; 00107 } 00108 00109 00110 } //namespace Math 00111 } //namespace ROOT