// uvtest.cc -- general test of all features of UnitVector class //---------------------------------------------------------------- #include "ZMutility/ZMenvironment.h" #include "PhysicsVectors/UnitVector.h" #include "CLHEP/Vector/ZMxpv.h" ZM_USING_NAMESPACE( zmpv ) #include "ZMutility/cmath" using CMATH_NAMESPACE::sqrt; using CMATH_NAMESPACE::fabs; using CMATH_NAMESPACE::log; using CMATH_NAMESPACE::tan; using CMATH_NAMESPACE::acos; #include "ZMutility/iostream" using std::cout; using std::cin; using std::endl; #include "ZMutility/iomanip" #include "ZMutility/PI.h" const Float8 PI = ZMpi; // Modifications since 2/1/98: // --------------------------- // // 2/6/98 Broke up into smaller subroutines for compilation speed. // // 2/6/98 Removed the remnants of the tests of inappropriate methods. // // 2/6/98 Modifying theta and phi had inadvertantly been excised as a test // of private methods; these tests are restored. // // 2/6/98 Major changes in list of signatures for set() that were being tested. // A wide variety are now tested in set3(). // // 2/6/98 Arithmetic on UnitVectors is now being tested; results are SV but // the methods should not be private. // // 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/23/98 Added test of nearness methods // // 5/09/98 Added test of rectify method // // 5/27/98 Modified tests for exact equality to tests within 10**-15 to // accomodate Linux and NT math inexactness. // // 6/10/98 Modified to use teeny because transcendentals on x86 machines // were not accurate enough that our 2E-15 tests were working // // 6/10/98 Modified to check for eta wrt SPaceVector direction vvu2 // // 6/15/98 Added namespace support // 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(); bool checkAxialRotation ( const UnitVector &x, const UnitVector &d, const Float8 &phi, const UnitVector& result ); 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-15 ); } bool problems (false); 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() ); // The UnitVector tests are a block so that we can later redeclare things cout << "\n\n"; cout << " ***********\n"; cout << " UnitVector\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(); #ifdef FUTURE_SVF // Assignment from other vector type cout << "\n --- Assignments from other vector types\n\n"; cout << "\n\n"; cout << " ***********\n"; cout << " UnitVectorF\n"; cout << " ***********\n"; cout << "\n\n"; cout << " ****************************\n"; cout << " UnitVector and UnitVectorF\n"; cout << " ****************************\n"; // various cross of accessors with types of results double ax; ax = fv1.z(); cout << "... should be 3 ==> " << ax << endl; float fx; fx = v2c.x(); cout << "... should be 20 ==> " << ax << endl; #endif // FUTURE_SVF 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 uvtest main() */ bool checkAxialRotation ( const UnitVector &x, const UnitVector &d, const double& phi, const UnitVector& result ) { // Rotate x by phi around d; return true if result does not match answer // within .0001 tolerance; print out appropriate messages. UnitVector xu (x); xu.rotate(d, phi); if ( xu.isNear(result, .0002) ) { cout << "OK: " << x << " Rotated " << phi << " about" << d << "\n = " << xu << " ~ " << result << endl; } else { cout << "?????" << x << " Rotated " << phi << " about" << d << "\n = " << xu << " not near " << result << endl; return true; } // Now form a new UnitVector as xu rotated about the axis by -phi; // again check that the result is proper (matches x) very tightly UnitVector yu = rotationOf(xu, d, -phi); if ( yu.isNear(x, .000001) ) { cout << "OK: " << xu << " Rotated " << -phi << " about" << d << "\n = " << yu << " ~ " << x << endl; } else { cout << "?????" << xu << " Rotated " << -phi << " about" << d << "\n = " << yu << " 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 UnitVector.cc is a highly // non-trivial test. } /* checkAxialRotation */ UnitVector vans; UnitVector u1 ( 90, DEGREES, 90, DEGREES); UnitVector u2xx (0,1,0); UnitVector u2 ( 30, DEGREES, 60, DEGREES ); const UnitVector u2c ( 30, DEGREES, PI/3, RADIANS ); UnitVector ua2 ( PI/6, RADIANS, -60, DEGREES ); const UnitVector ua2c ( PI/6, RADIANS, -PI/3, RADIANS ); UnitVector ub2 ( 1.31696, ETA, 60, DEGREES ); const UnitVector ub2c ( 1.31696, ETA, PI/3, RADIANS ); UnitVector u1c (4,-3,2); UnitVector u3 (8,-6,4); SpaceVector v1 ( -4, 12, -3); UnitVector u4 (v1); void constructors() { // Stretching excercises: construct and output a couple of vectors. cout << "\n --- Constructors and Output\n\n"; // Note: Constructor calls were moved global when breaking up into // smaller routines. These tests verify proper working of those. cout << "{ Default Constructor }\n"; if (vans.isNear(UnitVector(0,0,1))) { cout << "OK: default UnitVector is " << vans << endl; } else { cout << "????? default UnitVector is " << vans << endl; problems=true; } cout << "\n{ Spherical Constructor }\n"; if (u1.isNear(u2xx)) { cout << "OK: Constructed u1 = " << u1 << " == u2xx " << u2xx << endl; } else { cout << "????? Constructed u1 = " << u1 << " != u2xxx " << u2xx << endl; problems=true; } if ( u2.isNear(SpaceVector(2.5, 4.33, 8.66)/10,.01) ) { cout << "OK: constructed " << u2.r() << ", " << u2.theta()*(180/PI) << "degrees, " << u2.phi()*(180/PI) << " degrees = " << u2 << endl; } else { cout << "????? constructed " << u2.r() << ", " << u2.theta()*(180/PI) << ", " << u2.phi()*(180/PI) << " = " << u2 << endl; problems=true; } if ( u2c.isNear(SpaceVector(5.0, 8.66, 17.32)/20,.01) ) { cout << "OK: constructed const " << u2c.r() << ", " << u2c.theta()*(180/PI) << "degrees, " << u2c.phi()*(180/PI) << " degrees = " << u2c << endl; } else { cout << "????? constructed const " << u2c.r() << ", " << u2c.theta()*(180/PI) << ", " << u2c.phi()*(180/PI) << " = " << u2c << endl; problems=true; } if ( ua2.isNear(SpaceVector(2.5, -4.33, 8.66)/10,.01) ) { cout << "OK: constructed " << ua2.r() << ", " << ua2.theta()*(180/PI) << "degrees, " << ua2.phi()*(180/PI) << " degrees = " << ua2 << endl; } else { cout << "????? constructed " << ua2.r() << ", " << ua2.theta()*(180/PI) << ", " << ua2.phi()*(180/PI) << " = " << ua2 << endl; problems=true; } if ( ua2c.isNear(SpaceVector(5.0, -8.66, 17.32)/20,.01) ) { cout << "OK: constructed const " << ua2c.r() << ", " << ua2c.theta()*(180/PI) << "degrees, " << ua2c.phi()*(180/PI) << " degrees = " << ua2c << endl; } else { cout << "????? constructed const " << ua2c.r() << ", " << ua2c.theta()*(180/PI) << ", " << ua2c.phi()*(180/PI) << " = " << ua2c << endl; problems=true; } if ( ub2.isNear(SpaceVector(2.5, 4.33, 8.66)/10,.01) ) { cout << "OK: constructed " << ub2.r() << ", eta = " << ub2.eta()*(180/PI) << " --> " << ub2.theta()*(180/PI) << "degrees, " << ub2.phi()*(180/PI) << " degrees = " << ub2 << endl; } else { cout << "????? constructed " << ub2.r() << ", " << ub2.theta()*(180/PI) << ", " << ub2.phi()*(180/PI) << " = " << ub2 << endl; problems=true; } if ( ub2c.isNear(SpaceVector(5.0, 8.66, 17.32)/20,.01) ) { cout << "OK: constructed const " << ub2c.r() << ", " << ub2c.theta()*(180/PI) << " eta= " << ub2c.eta() << ", " << ub2c.phi()*(180/PI) << " = " << ub2c << endl; } else { cout << "????? constructed const " << ub2c.r() << ", " << ub2c.theta()*(180/PI) << " eta= " << ub2c.eta() << ", " << ub2c.phi()*(180/PI) << " = " << ub2c << endl; problems=true; } cout << "\n{ Normalized-Cartesian Constructor }\n"; if (u1c.isNear(u3)) { cout << "OK: Constructed const u1c = " << u1c << " == u3 " << u3 << endl; } else { cout << "????? Constructed u1c = " << u1c << " != u3 " << u3 << endl; problems=true; } if ((13*u4).isNear(v1)) { cout << "OK: Constructed " << u4 << " from " << v1 << endl; } else { cout << "????? Constructed " << u4 << " from " << v1 << endl; problems=true; } } // constructors() UnitVector u5 ( 12, 3, -4); UnitVector u6 ( 1, 1, 0 ); double xx = u5.x(); double yy = u5.y(); double zz = u5.z(); void accessors() { // Accessors // Straight x y z cout << "\n --- Accessors\n\n"; cout << "{ get Cartesian components }\n"; if ( eq(13.*xx, 12) && eq(13.*yy, 3) && eq(13*zz, -4) ) { cout << "OK: u5 = " << u5 << endl; } else { cout << "????? u5 = " << u5 << endl; problems=true; } cout << "\n{ output Cartesian components }\n"; cout << " ... ostream<< x, y, z ==> " << u5.x() <<"," << u5.y() <<"," << u5.z() << endl; double rr = SpaceVector(4,-3,2).mag(); xx = u1c.x(); yy = u1c.y(); zz = u1c.z(); SpaceVector tmpv ( xx, yy, zz ); if ( (rr * tmpv).isNear(SpaceVector(4,-3,2)) ) { cout << "OK: Components are " << xx <<"," << yy <<"," << zz << endl; } else { cout << "?????Components are " << xx <<"," << yy <<"," << zz << endl; problems=true; } cout << " Output Check ... should match above components ==>\n " << u1c.x() <<"," << u1c.y() <<"," << u1c.z() << endl; cout << "\n{ use Cartesian components in expressions }\n"; if ( eq( (6 * u5.z()) , 6. * (-4./13.) ) ) { cout << "OK: 6*u5.z() ==> " << 6*u5.z() << endl; } else { cout << "????? 6*u5.z() ==> " << 6*u5.z() << endl; problems=true; } if ( eq( (u5.y()*5) , 5. * (3./13.) ) ){ cout << "OK: u5.y()*5 ==> " << u5.y()*5 << endl; } else { cout << "????? u5.y()*5 ==> " << u5.y()*5 << endl; problems=true; } if ( eq ( (6 * u1c.z()) , 6. * (2./sqrt(29.)) ) ) { cout << "OK: 6*u1c.z() ==> " << 6*u1c.z() << endl; } else { cout << "????? 6*u1c.z() ==> " << 6*u1c.z() << endl; problems=true; } if ( eq ( (u1c.x()*5) , 5. * (4./sqrt(29.)) ) ) { cout << "OK: u1c.x()*5 ==> " << u1c.x()*5 << endl; } else { cout << "????? u1c.x()*5 ==> " << u1c.x()*5 << endl; problems=true; } // Observation: u1.x_ directly does not compile, as desired. #ifdef ZAP1 cout << u1.x_; #endif // ZAP1 // Comparison of accessors cout << "\n{ compare Cartesian components }\n"; if ( u5.x() == 12./13. ) { cout << "OK: Equality comparison on an accessor works\n"; } else { cout << "????? Equality comparison on an accessor FAILS!!!!!\n"; problems=true; } if ( u5.x() != u5.y() ) { cout << "OK: comp != comp works\n"; } else { cout << "????? comp != comp FAILS!!!!\n"; problems=true; } if ( u5.x() == u5.y() ) { cout << "????? comp == comp FAILS!!!!\n"; problems=true; } else { cout << "OK: comp == comp works\n"; } // reminder: u5 is based on ( 12, 3, -4 ); // u4 is based on ( -4, 12, -3); // u1c is based on ( 4, -3, 2); if ( u5.y() > u5.x() ) { cout << "????? comp > comp FAILS!!!!\n"; problems=true; } else { cout << "OK: comp > comp works\n"; } if ( u5.x() < u5.z() ) { cout << "?????comp < comp FAILS!!!!\n"; problems=true; } else { cout << "OK: comp < comp works\n"; } double xeq = 12.0/13.; if ( xeq == u5.x() ) { cout << "OK: Scalar == comp works\n"; } else { cout << "?????Scalar == comp FAILS!!!!\n"; problems=true; } const double ceq = -4.0/13.; if ( ceq != u5.z() ) { cout << "????? Scalar != comp FAILS!!!!\n"; problems=true; } else { cout << "OK: Scalar != comp works\n"; } cout << "\n{ get Cartesian components of const }\n"; double norm432 = 1./sqrt(29.); if ( eq(u1c.x(), 4.*norm432) ) { cout << "OK: Equality comparison on a const vector accessor works\n"; } else { cout << "????? Equality comparison on a const vector accessor FAILS!!!!!\n"; cout << "\nu1c.x() = " << u1c.x() << endl; cout << "\n4*norm432 = " << 4*norm432 << endl; problems=true; } if ( u1c.x() != u1c.y() ) { cout << "OK: comp != comp works\n"; } else { cout << "?????comp != comp FAILS!!!!\n"; problems=true; } if ( u5.x() == u5.y() ) { cout << "?????comp == comp FAILS!!!!\n"; problems=true; } else { cout << "OK: comp == comp works\n"; } // u1c is based on ( 4, -3, 2); if ( u1c.y() > u1c.x() ) { cout << "????? comp > comp FAILS!!!!\n"; problems=true; } else { cout << "OK: comp > comp works\n"; } if ( u1c.z() < u1c.x() ) { cout << "OK: comp < comp works\n"; } else { cout << "????? comp < comp FAILS!!!!\n"; problems=true; } // Getting different components cout << "\n{ get Spherical components }\n"; double p = u6.theta(); if ( close(p,PI/2)) { cout << "OK: theta = " << p *180/PI<< endl; } else { cout << "????? theta = " << p*180/PI << endl; problems=true; } double a = u6.phi(); if ( close(a,PI/4)) { cout << "OK: phi = " << a *180/PI<< endl; } else { cout << "????? phi = " << a *180/PI << endl; problems=true; } const UnitVector u2cc ( 0, 1, -1 ); p = u2cc.theta(); if ( close(p,135*PI/180)) { cout << "OK: theta = " << p *180/PI<< endl; } else { cout << "????? theta = " << p *180/PI << endl; problems=true; } a = u2cc.phi(); if ( close (a,PI/2)) { cout << "OK: phi = " << a *180/PI << endl; } else { cout << "????? phi = " << a*180/PI << endl; problems=true; } xx = u4.r(); if ( xx == 1.0 ) { cout << "OK: r = " << xx << endl; } else { cout << "????? r = " << xx << endl; problems=true; } cout << "\nOutput Check ... 1.0 ==> " << u5.r() << endl; // u4 is based on ( -4, 12, -3); xx = u4.rho(); if ( close(xx,sqrt(160.)/13.) ) { cout << "OK: rho = " << xx << endl; } else { cout << "????? rho = " << xx << endl; problems=true; } cout << "\nOutput Check ... ==> " << u4.rho() << endl; cout << "\n{ eta }\n"; // u6 is based on ( 1, 1, 0 ); cout << " ..." << u6 << " has theta, phi of " << u6.theta() << u6.phi() << endl; xx = u6.eta(); if ( eq(xx,0) ) { cout << "OK: " << u6 << ".eta() = " << xx << endl; } else { cout << "????? " << u6 << ".eta() = " << xx << endl; problems=true; } cout << " ... eta for UnitVector(1, 2, -1000) = " << UnitVector(1, 2, -1000).eta() << endl; cout << " ... eta for SpaceVector(1, 2, 1000) = " << UnitVector(1, 2, 1000).eta() << endl; // Setting different components cout << "\n{ Set Spherical components }\n"; // u4 is based on ( -4, 12, -3); double p4 = u4.theta(); double a4 = u4.phi(); cout << " ... Theta for " << u4 << " is " << p4 << endl; cout << " ... Phi for " << u4 << " is " << a4 << endl; u4.setTheta( 0 ); if (u4.isNear(UnitVector(0,0,13))) { cout << "OK: with u4.theta() = 0 u4 is " << u4 << endl; } else { cout << "????with u4.theta() = 0 u4 is " << u4 << endl; problems=true; } cout << " " << u4 << " ... Reset theta and then phi:\n "; u4.setTheta( p4 ); cout << u4; u4.setPhi( a4 ); cout << u4 << endl; vans.set (-4, 12, -3); if (u4.isNear(vans)) { cout << "OK: " << u4 << " = " << vans << endl; } else { cout << "????? " << u4 << " is not " << vans << endl; problems = true; } cout << "\n Skipping test of private methods -- " "setting radial component." << endl; // ADDED 980319 MF: UnitVector ue ( 4, 3, 12 ); ue.setEta( 3.2 ); if ( close(ue.eta(),3.2) && close(ue.dot(ue), 1) ) { cout << "OK: " << ue << " has eta = " << ue.eta() << " and ue.dot(ue) = " << ue.dot(ue) << endl; } else { cout << "????? " << ue << " has eta = " << ue.eta() << " and ue.dot(ue) = " << ue.dot(ue) << endl; problems = true; } #ifdef ZAP2 u2c.theta() = .05; // Properly fails, with message // `WontAlter_Theta_of_a_const' does not define operator= #endif // ZAP2 #ifdef NOTYET_WAIT_FOR_EXCEPTIONS // Setting components, with vector along Z axis (pos and neg) #endif // NOTYET_WAIT_FOR_EXCEPTIONS // setCylTheta, setCylEta cout << "\n Skipping test of private methods -- " "setting cylindrical components." << endl; } // accessors() void modifiers() { cout << "\n --- Component Modify/Assigners\n\n"; // Component modifiers ( u.x() += Scalar, etc.) cout << "\n{ Modifying Spherical angle components }\n"; u2.setTheta( PI/6 ); u2.setPhi( -PI/6 ); cout << " u2.phi() += 50 degrees\n"; u2.setPhi( u2.phi() + 50*PI/180 ); vans.set ( 30, DEGREES, 20, DEGREES); if (u2.isNear(vans)) { cout << "OK : " << u2.theta()*180/PI << " degrees " << u2.phi()*180/PI << " degrees " << endl; } else { cout << "OK : " << u2.theta()*180/PI << " degrees " << u2.phi()*180/PI << " degrees " << endl; problems = true; } cout << " ( " << u2.r() << ", " << u2.theta()*180/PI << " degrees " << u2.phi()*180/PI << " degrees "; cout << " u2.theta() *= 2\n"; u2.setTheta( u2.theta() * 2 ); vans.set ( 60, DEGREES, 20, DEGREES); if (u2.isNear(vans)) { cout << "OK : " << u2.theta()*180/PI << " degrees " << u2.phi()*180/PI << " degrees "<< endl; } else { cout << "OK : " << u2.theta()*180/PI << " degrees " << u2.phi()*180/PI << " degrees "<< endl; problems = true; } #ifdef NOTYET_WAIT_FOR_EXCEPTIONS cout << "\n{ 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 } // modifiers() void set3() { cout << "\n --- Setting three components at once\n\n"; cout << "\n{ Cartesian components }\n"; u1.set ( 5, 7, 9 ); if ( close(u1.x(), 5./sqrt(155.)) && close(u1.y(), 7./sqrt(155.)) && close(u1.z(), 9./sqrt(155.)) ) { cout << "OK: " << u1 << endl; } else { cout << "?????" << u1 << endl; problems = true; } cout << "{ Spherical components (angles only) }\n"; u1.set ( 42, DEGREES, 99, DEGREES ); if ( close(u1.r(),1) && close(u1.theta(),42*PI/180) && close(u1.phi(),99*PI/180) ) { cout << "OK: " << u1.r() << ", " << u1.theta()*180/PI << " degrees " << u1.phi()*180/PI << " degrees\n"; } else { cout << "?????" << u1.r() << ", " << u1.theta()*180/PI << " degrees " << u1.phi()*180/PI << " degrees\n"; } u1.set ( 1.2, RADIANS, -467, DEGREES ); if ( close(u1.r(),1) && close(u1.theta(),1.2) && close(u1.phi(),-107*PI/180) ) { cout << "OK: " << u1.r() << ", " << u1.theta() << " radians " << u1.phi()*180/PI << " degrees\n"; } else { cout << "?????" << u1.r() << ", " << u1.theta() << " radians " << u1.phi()*180/PI << " degrees\n"; problems = true; } u1.set ( 3, ETA, .5, RADIANS ); if ( close(u1.r(),1) && close(u1.eta(),3) && close(u1.phi(),.5) ) { cout << "OK: " << u1.r() << ", eta: " << u1.eta() << " " << u1.phi() << " radians\n"; } else { cout << "?????" << u1.r() << ", eta: " << u1.eta() << " " << u1.phi() << " radians\n"; problems = true; } #ifdef ZAP_SET3 u1.set ( 5, 42, DEGREES, 99, DEGREES ); #endif // ZAP_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 } // set3() const UnitVector u3c( 3, 5, 7); const UnitVector u4c( 6, 2, 8); SpaceVector v3; SpaceVector v3ans; void arithmetic() { // Arithmetic operations cout << "\n --- Arithmetic\n\n"; u1.set( 3, 5, 7); u2.set( 6, 2, 8); SpaceVector u1v = u1; SpaceVector u2v = u2; cout << "\n{ + and - }\n"; // These checks rely on the SpaceVector arithmetic working properly; // those checks were done explicitly in svtest.cc. v3 = u1 + u2; v3ans = u1v + u2v; if (v3.isNear(v3ans)) { cout << "OK: " << u1 << " + " << u2 << " = " << v3 << endl; } else { cout << "?????" << u1 << " + " << u2 << " = " << v3 << endl; problems = true; } v3 = u1 - u2; v3ans = u1v - u2v; if (v3.isNear(v3ans)) { cout << "OK: " << u1 << " - " << u2 << " = " << v3 << endl; } else { cout << "?????" << u1 << " - " << u2 << " = " << v3 << endl; problems = true; } v3 = u4c + u3c; u1v = u3c; u2v = u4c; v3ans = u1v + u2v; if (v3.isNear(v3ans)) { cout << "OK: " << u4c << " + " << u3c << " = " << v3 << endl; } else { cout << "?????" << u4c << " + " << u3c << " = " << v3 << endl; problems = true; } v3 = u4c - u3c; v3ans = u2v - u1v; if (v3.isNear(v3ans)) { cout << "OK: " << u4c << " - " << u3c << " = " << v3 << endl; } else { cout << "?????" << u4c << " - " << u3c << " = " << v3 << endl; problems = true; } cout << "\n{ Unary minus }\n"; // This one gives a UnitVector result: u3 = - u2; vans.set ( -6, -2, -8 ); if (u3.isNear(vans)) { cout << "OK: " << " - " << u2 << " = " << u3 << endl; } else { cout << "?????" << " - " << u2 << " = " << u3 << endl; problems = true; } // This one gives a UnitVector result: u3 = - u3c; vans.set ( -3, -5, -7 ); if (u3.isNear(vans)) { cout << "OK: " << " - " << u3c << " = " << u3 << endl; } else { cout << "?????" << " - " << u3c << " = " << u3 << endl; problems = true; } cout << "\n{ * and / }\n"; double pointfive = .5; v3 = u2 * pointfive; v3ans.set ( 6, 2, 8 ); v3ans /= v3ans.mag(); v3ans /= 2; if (v3.isNear(v3ans)) { cout << "OK: " << u2 << " * .5 = " << v3 << endl; } else { cout << "?????" << u2 << " * .5 = " << v3 << endl; problems = true; } v3 = u2 / 4; v3ans = .25 * SpaceVector( 6, 2, 8) / sqrt(104.); if (v3.isNear(v3ans)) { cout << "OK: " << u2 << " / 4 = " << v3 << endl; } else { cout << "?????" << u2 << " / 4 = " << v3 << endl; problems = true; } v3 = 1.5 * u2; v3ans = 1.5 * SpaceVector( 6, 2, 8) / sqrt(104.); if (v3.isNear(v3ans)) { cout << "OK: 1.5 * " << u2 << " = " << v3 << endl; } else { cout << "????? 1.5 * " << u2 << " = " << v3 << endl; problems = true; } v3 = u4c * .5; v3ans = .5 * SpaceVector( 6, 2, 8) / sqrt(104.); if (v3.isNear(v3ans)) { cout << "OK: " << u4c << " * .5 = " << v3 << endl; } else { cout << "?????" << u4c << " * .5 = " << v3 << endl; problems = true; } v3 = u4c / .5; v3ans = 2.0 * SpaceVector( 6, 2, 8) / sqrt(104.); if (v3.isNear(v3ans)) { cout << "OK: " << u4c << " / .5 = " << v3 << endl; } else { cout << "?????" << u4c << " / .5 = " << v3 << endl; problems = true; } v3 = 2.5 * u4c; v3ans = 2.5 * SpaceVector( 6, 2, 8) / sqrt(104.); if (v3.isNear(v3ans)) { cout << "OK: 2.5 * " << u2 << " = " << v3 << endl; } else { cout << "????? 2.5 * " << u2 << " = " << v3 << endl; problems = true; } cout << "\n{ Linear expression }\n"; u3.set(9,3,12); v3 = u3 + 5 * u1 + u4c/2; u1v = 5 * u1; u2v = u4c/2.; v3ans = u3 + u1v + u2v; if (v3.isNear(v3ans)) { cout << "OK: " << u3 << " + 5 * " << u1 << " + " << u4c << " / 2 = " << v3 << endl; } else { cout << "?????" << u3 << " + 5 * " << u1 << " + " << u4c << " / 2 = " << v3 << endl; problems = true; } cout << "\n{ rectify } \n"; UnitVector uq ( .5, -.6, .4 ); cout << " ... Start from UnitVector " << uq << endl; SpaceVector * vuq = reinterpret_cast (&uq); (*vuq).setX( (*vuq).x() + .02 ); (*vuq).setY( (*vuq).y() + .05 ); (*vuq).setZ( (*vuq).z() + .08 ); cout << " UnitVector " << uq << " is intentionally inexact\n"; SpaceVector sq = uq; cout << " Length of that is " << sq.mag() << endl; uq.rectify(); cout << " ... rectified UnitVector is " << uq << endl; sq = uq; cout << " Length of that is " << sq.mag() << endl; if ( close (sq.mag(), 1.0) ) { cout << "OK: rectified length is 1\n"; } else { cout << "????? rectified length is not 1\n"; problems = true; } cout << "\n Skipping test of private methods -- " "modify/assign on UnitVector." << endl; } // arithmetic() void comparisons() { // Comparisons cout << "\n --- Comparisons\n\n"; u1.set (2, 4, 6); u2.set (1, 5, 9); // Previously: // const UnitVector u3c( 3, 5, 7); // const UnitVector u4c( 6, 2, 8); u5 = u1; cout << "\n{ == and != }\n"; if ( ! (u1 == u5) ) { cout << "????? u1 == u5 false erroneously! u1 = " << u1 << " u5 = " << u5 << endl; problems = true; } cout << "OK: " << u1 << " == " << u5 << endl; if ( u1 == u2 ) { cout << "????? u1 == u2 true erroneously! u1 = " << u1 << " u2 = " << u2 << endl; problems = true; } cout << "OK: " << "NOT " << u1 << " == " << u2 << endl; if ( u1 != u5 ) { cout << "????? u1 != u5 true erroneously! u1 = " << u1 << " u5 = " << u5 << endl; problems = true; } cout << "OK: " << "NOT " << u1 << " != " << u5 << endl; if ( ! (u1 != u2) ) { cout << "????? u1 != u2 false erroneously! u1 = " << u1 << " u2 = " << u2 << endl; problems = true; } cout << "OK: " << u1 << " != " << u2 << endl; if ( u1 == u3c ) { cout << "????? u1 == u3c true erroneously! u1 = " << u1 << " u3c = " << u3c << endl; problems = true; } cout << "OK: " << "NOT " << u1 << " == " << u2 << endl; if ( ! (u3c != u4c) ) { cout << "????? u3c != u4c false erroneously! u3c = " << u3c << " u4c = " << u4c << endl; problems = true; } cout << "OK: " << u3c << " != " << u4c << endl; cout << "\n{ Dictionary Ordering: >, <, >=, <= }\n"; // Modification 2/6/98: Changed a,b,c to c,b,a to reflect change in ordering. u1.set (6, 4, 2); // CHANGED 2/6/98 from (2,4,6) u2.set (9, 5, 1); // CHANGED 2/6/98 from (1,5,9) // Previously: // const UnitVector u3c( 3, 5, 7); const UnitVector u33c (7,5,3); // CHANGED 2/6/98 u5.set ( 3, 5, 2); // CHANGED 2/6/98 from (2,5,3) u6 = u1; if ( ! (u1 < u5) ) { cout << "????? u1 < u5 false erroneously! u1 = " << u1 << " u5 = " << u5 << endl; problems = true; } cout << "OK: " << u1 << " < " << u5 << endl; if ( u1 < u6 ) { cout << "????? u1 < u6 true erroneously! u1 = " << u1 << " u6 = " << u6 << endl; problems = true; } cout << "OK: " << "NOT " << u1 << " < " << u6 << endl; if ( u1 > u6 ) { cout << "????? u1 > u6 true erroneously! u1 = " << u1 << " u6 = " << u6 << endl; problems = true; } cout << "OK: " << "NOT " << u1 << " > " << u6 << endl; if ( u5 >= u3c ) { cout << "????? u5 >= u3c true erroneously! u5 = " << u5 << " u5 = " << u3c << endl; problems = true; } cout << "OK: " << "NOT " << u5 << " >= " << u3c << endl; if ( ! (u1 <= u5 ) ) { cout << "????? u1 <= u5 false erroneously! u1 = " << u1 << " u5 = " << u5 << endl; problems = true; } cout << "OK: " << u1 << " <= " << u5 << endl; if ( ! (u1 <= u6 ) ) { cout << "????? u1 <= u6 false erroneously! u1 = " << u1 << " u6 = " << u6 << endl; problems = true; } cout << "OK: " << u1 << " <= " << u6 << endl; if ( ! (u6 <= u1 ) ) { cout << "????? u6 <= u1 true erroneously! u6 = " << u6 << " u1 = " << u1 << endl; problems = true; } cout << "OK: " << u6 << " <= " << u1 << endl; } // comparisons() SpaceVector v2( 4, 8, 7 ); void tolerance() { // Comparisons within tolerance cout << "\n --- IsNear(): Comparisons within tolerance\n\n"; cout << "\n{ dot and mag2() as used by isNear }\n"; u1.set( 3, 5, 1); u2.set( 4, 8, 7); xx = u1.dot(u2); if (eq(xx , (u1.x()*u2.x()+u1.y()*u2.y()+u1.z()*u2.z()))) { cout << "OK: " << u1 << ".dot" << u2 << " = " << xx << endl; } else { cout << "?????" << u1 << ".dot" << u2 << " = " << xx << endl; problems = true; } xx = u1.mag2(); if (xx==1) { cout << "OK: " << u1 << ".mag2() = " << xx << endl; } else { cout << "?????" << u1 << ".mag2() = " << xx << endl; problems = true; } cout << "\n{ isNear() with default tolerance }\n"; v1.set( 3, 5, 1 ); double toler = SpaceVector::setTolerance(1.); cout << "Default tolerance starts at (expected 2.22e-14) ==> " << toler << endl; SpaceVector::setTolerance(toler); v3 = v1 + .2 * toler * v2; u1 = v1; u3 = v3; if (u3 != u1) { cout << "OK: " << "adding .2 * tolerance * times " << u2 << " to u1 = " << u1 << " gives U3 != U1\n"; } else { cout << "????? adding .2 * tolerance * times " << u2 << " to u1 = " << u1 << " gives U3 == U1\n"; problems = true; } if ( u3 .isNear(u1) ) { cout << "OK: " << "Added .2 * tolerance * u2 -- isNear TRUE\n"; } else { cout << "????? Added .2 * tolerance -- isNear FALSE\n"; problems = true; } v3 = v1 + .6 * toler * v2; u3 = v3; if ( u1 .isNear(u3) ) { cout << "OK: Added .6 * tolerance * -- isNear FALSE\n"; } else { cout << "????? " << "Added .6 * tolerance * u2 -- isNear\n"; problems = true; } v3 = v1 + 1.2 * toler * v2; u3 = v3; if ( u3.isNear(u1) ) { cout << "Added 1.2 * tolerance * u2 -- isNear TRUE\n"; } else { cout << "Added 1.2 * tolerance * u2 -- isNear FALSE\n"; } v3 = v1 + 1.8 * toler * v2; u3 = v3; if ( u3.isNear(u1) ) { cout << "????? Added 1.8 * tolerance * u2 -- isNear TRUE\n"; problems = true; } else { cout << "OK: " << "Added 1.8 * tolerance * u2 -- isNear FALSE\n"; } cout << "\n{ isNear() with specified tolerance }\n"; v1.set ( 1, 3, 6 ); u1 = v1; double m = sqrt(v1.mag2()); cout << " .... dealing with a vector" << u1 << "of magnitude " << m << endl; double small = .0001; v2 = v1 + small * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "with (default) tolerance " << toler << ", " << u1 << " isNear " << u2 << endl; cout << "????? That should not be!\n"; problems = true; } else { cout << "OK: " << "with (default) tolerance " << toler << ",\n" " " << u1 << " is not Near " << u2 << endl; } if ( u2.isNear(u1) ) { cout << "with tolerance " << toler << ", " << u2 << " isNear " << u1 << endl; cout << "????? That should not be!\n"; problems = true; } else { cout << "OK: " << "with tolerance " << toler << ", " << u2 << " is not Near " << u1 << endl; } if ( u1.isNear(u2 , .5*small/m) ) { cout << "with tolerance " << .5*small/m << ", " << u1 << " isNear " << u2 << endl; cout << "????? That should not be!\n"; problems = true; } else { cout << "OK: " << "with tolerance " << .5*small/m << ", " << u1 << " is not Near " << u2 << endl; } if ( u2.isNear(u1 , small/m) ) { cout << "with tolerance " << small/m << ", " << u2 << " isNear " << u1 << endl; } else { cout << "with tolerance " << small/m << ", " << u2 << " is not Near " << u1 << endl; } if ( u1.isNear(u2 , 1.5*small/m) ) { cout << "with tolerance " << 1.5*small/m << ", " << u1 << " isNear " << u2 << endl; } else { cout << "with tolerance " << 1.5*small/m << ", " << u1 << " is not Near " << u2 << endl; } if ( u1.isNear(u2 , 2*small/m) ) { cout << "with tolerance " << 2*small/m << ", " << u1 << " isNear " << u2 << endl; } else { cout << "with tolerance " << 2*small/m << ", " << u1 << " is not Near " << u2 << endl; } if ( u1.isNear(u2 , 3.5*small/m) ) { cout << "OK: " << "with tolerance " << 3.5*small/m << ", " << u1 << " isNear " << u2 << endl; } else { cout << "with tolerance " << 3.5*small/m << ", " << u1 << " is not Near " << u2 << endl; cout << "????? That should not be!\n"; problems = true; } cout << "\n{ setTolerance }\n"; SpaceVector::setTolerance ( .1 ); cout << "tolerance set to .1\n"; v2 = v1 + .02 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "OK: " << "with new default tolerance " << u1 << " isNear " << u2 << endl; } else { cout << "with new default tolerance " << u1 << " is not Near " << u2 << endl; cout << "????? That should not be!\n"; problems = true; } v2 = v1 + .05 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "OK: " << "with new default tolerance " << u1 << " isNear " << u2 << endl; } else { cout << "with new default tolerance " << u1 << " is not Near " << u2 << endl; cout << "????? That should not be!\n"; problems = true; } v2 = v1 + .1 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "OK: with new default tolerance " << u1 << " isNear " << u2 << endl; } else { cout << "????? with new default tolerance " << u1 << " is not Near " << u2 << endl; problems = true; } v2 = v1 + .2 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "OK: with new default tolerance " << u1 << " isNear " << u2 << endl; } else { cout << "?????with new default tolerance " << u1 << " is not Near " << u2 << endl; problems = true; } v2 = v1 + .3 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "OK: with new default tolerance " << u1 << " isNear " << u2 << endl; } else { cout << "????? with new default tolerance " << u1 << " is not Near " << u2 << endl; problems = true; } v2 = v1 + .45 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "?????with new default tolerance " << u1 << " isNear " << u2 << endl; problems = true; } else { cout << "OK: with new default tolerance " << u1 << " is not Near " << u2 << endl; } v2 = v1 + .6 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "?????with new default tolerance " << u1 << " isNear " << u2 << endl; problems = true; } else { cout << "OK: with new default tolerance " << u1 << " is not Near " << u2 << endl; } v2 = v1 + 1.2 * SpaceVector( 1, -1, 1); u2 = v2; if ( u1.isNear(u2) ) { cout << "????? with new default tolerance " << u1 << " isNear " << u2 << endl; problems = true; } else { cout << "OK: " << "with new default tolerance " << u1 << " is not Near " << u2 << endl; } } // tolerance() void nearness() { cout << "\n --- Nearness methods \n\n"; UnitVector un1 ( 10, 20, 30 ); UnitVector un2 ( 1, 2, 5 ); if ( close (un1.howNear(un2), sqrt(2-2*sqrt(420.)/21)) ) { cout << "OK: " << un1 << ".howNear" << un2 << "= " << un1.howNear(un2) << endl; } else { cout << "????? " << un1 << ".howNear" << un2 << "= " << un1.howNear(un2) << endl; problems = true; } un1.set(10, 20, -11); if ( close (un1.howNear(un2), sqrt(2+10/sqrt(30.*621))) ) { cout << "OK: " << un1 << ".howNear" << un2 << "= " << un1.howNear(un2) << endl; } else { cout << "????? " << un1 << ".howNear" << un2 << "= " << un1.howNear(un2) << endl; problems = true; } un1.set ( 1, 2, -1 ); SpaceVector vn2 ( 1, 2, 5 ); if ( close (un1.howNear(vn2), sqrt(31.)) ) { cout << "OK: " << un1 << ".howNear" << vn2 << "= " << un1.howNear(vn2) << endl; } else { cout << "????? " << un1 << ".howNear" << vn2 << "= " << un1.howNear(vn2) << endl; problems = true; } if ( un1.howParallel(vn2) == 1 ) { cout << "OK: " << un1 << ".howParallel" << vn2 << "= " << un1.howParallel(vn2) << endl; } else { cout << "????? " << un1 << ".howParallel" << vn2 << "= " << un1.howParallel(vn2) << endl; problems = true; } if ( eq(un1.howOrthogonal(vn2) , 0) ) { cout << "OK: " << un1 << ".howOrthogonal" << vn2 << "= " << un1.howOrthogonal(vn2) << endl; } else { cout << "????? " << un1 << ".howOrthogonal" << vn2 << "= " << un1.howOrthogonal(vn2) << endl; problems = true; } un1.set ( 3, 4, -5 ); un2.set ( 4, 3, 2 ); if ( un1.howParallel(un2) == 1 ) { cout << "OK: " << un1 << ".howParallel" << un2 << "= " << un1.howParallel(un2) << endl; } else { cout << "????? << " << un1 << ".howParallel" << un2 << "= " << un1.howParallel(un2) << endl; problems = true; } if ( close(un1.howOrthogonal(un2), 14/sqrt(1254.)) ) { cout << "OK: " << un1 << ".howOrthogonal" << un2 << "= " << un1.howOrthogonal(un2) << endl; } else { cout << "????? " << un1 << ".howOrthogonal" << un2 << "= " << un1.howOrthogonal(un2) << endl; problems = true; } un1.set ( 6, 7, 8 ); vn2.set ( 8, 9, 10 ); if ( close(un1.howParallel(vn2),sqrt(24.)/191) ) { cout << "OK: " << un1 << ".howParallel" << vn2 << "= " << un1.howParallel(vn2) << endl; } else { cout << "????? << " << un1 << ".howParallel" << vn2 << "= " << un1.howParallel(vn2) << endl; problems = true; } if ( un1.howOrthogonal(vn2) == 1 ) { cout << "OK: << " << un1 << ".howOrthogonal" << vn2 << "= " << un1.howOrthogonal(vn2) << endl; } else { cout << "????? << " << un1 << ".howOrthogonal" << vn2 << "= " << un1.howOrthogonal(vn2) << endl; problems = true; } UnitVector un3 ( 1.2, ETA, .25, RADIANS ); SpaceVector vn4 ( 25, 0.9, ETA, -.15, RADIANS ); if ( close(un3.deltaR(vn4), .5) ) { cout << "OK: " << un3 << "\n .deltaR" << vn4 << "= " << un3.deltaR(vn4) << endl; } else { cout << "????? " << un3 << "\n .deltaR" << vn4 << "= " << un3.deltaR(vn4) << endl; problems = true; } } // nearness() void combining() { // Combining vectors cout << "\n --- Combining two vectors\n\n"; v1.set(2, 4, 5); u1 = v1; v2.set(6, 1, -4); u2 = v2; cout << "\n{ dot product }\n"; double udu = u1.dot(u2); double vdv = v1.dot(v2)/(v1.mag()*v2.mag()); if ( close(udu, vdv)) { cout << "OK: " << u1 << ".dot" << u2 << " ==> " << udu << endl; } else { cout << "????? " << u1 << ".dot" << u2 << " ==> " << udu << endl; problems = true; } udu = u2.dot(u2); vdv = v2.dot(v2)/v2.mag2(); if ( close(udu, vdv )) { cout << "OK: " << u2 << ".dot" << u2 << " ==> " << udu << endl; } else { cout << "????? " << u2 << ".dot" << u2 << " ==> " << udu << endl; problems = true; } cout << "\n{ cross product }\n"; u3 = u1.cross(u2); if ( u3.isNear(UnitVector ( -21, 38, -22 )) ) { cout << "OK: " << u1 << ".cross" << u2 << " ==> " << u3 << endl; } else { cout << "????? " << u1 << ".cross" << u2 << " ==> " << u3 << endl; problems = true; } SpaceVector v3t=u3; u3 = u2.cross(u1); SpaceVector v3c=u3; SpaceVector v3s=v3t + v3c; if ( v3s.mag2() < .000000001 * u2.mag2() ) { cout << "OK: " << u1 << ".cross" << u2 << " + " << u2 << ".cross" << u1 << " ==> (0, 0, 0)" << endl; } else { cout << "?????: " << u1 << ".cross" << u2 << " + " << u2 << ".cross" << u1 << " ==> " << v3s << endl; cout << "\nu1.cross(u2) = " << v3t << endl; cout << "\nu2.cross(u1) = " << v3c << endl; cout << "\nsum = " << v3s << endl; problems = true; } double diff = (u1.cross(u1).mag()); if ( close( diff, 0) ) { cout << "OK: " << u1 << ".cross" << u1 << " ==> (0, 0, 0) " << endl; } else { cout << "?????: " << u1 << ".cross" << u1 << " ==> " << (u1.cross(u1)) << endl; problems = true; } cout << "\n{ diff2 }\n"; u1.set(2, 4, 5); v1 = u1; u2.set(6, 1, -4); v2 = u2; udu = u1.diff2(u2); vdv = v1.diff2(v2); if (eq(udu, vdv)) { cout << "OK: " << u1 << ".diff2" << u2 << " ==> " << udu << endl; } else { cout << "????? " << u1 << ".diff2" << u2 << " ==> " << udu << endl; problems = true; } cout << "\n{ 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; } u1.set(2, 4, 5); u2.set(4, 8, 10); u5.set(-10, 4, 8); //if ( u1.isParallel(u2,0) ) { // Fails with KCC, OK with gcc if ( u1.isParallel(u2,1.0E-16) ) { // JMM testing cout << "OK: " << u1 << ".isParallel" << u2 << "exactly" << endl; } else { cout << "?????: " << u1 << ".isParallel" << u2 << " FALSE" << endl; problems = true; } u4 = u2; SpaceVector v4 = u4; SpaceVector v5 = u5; v2 = v4 + .3 * percent * v5; u2 = v2; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (.3*percent*u5) << " FALSE" << endl; problems = true; } v2 = v4 + .6 * percent * v5; u2 = v2; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (.6*percent*u5) << " FALSE" << endl; problems = true; } v2 = v4 + .9 * percent * v5; u2 = v2; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (.9*percent*u5) << " FALSE" << endl; problems = true; } v2 = v4 + 1.0 * percent * v5; u2 = v2; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (1.0*percent*u5) << " 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; u2 = v2; if ( u1.isParallel(u2) ) { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (1.1*percent*u5) << endl; problems = true; } else { cout << "OK: NOT " << u1 << ".isParallel" << u2 << endl; } v2 = v4 + 1.2 * percent * v5; u2 = v2; if ( u1.isParallel(u2) ) { cout << "?????: " << u1 << ".isParallel\n " << u4 << " + " << (1.2*percent*u5) << endl; problems = true; } else { cout << "OK: NOT " << u1 << ".isParallel" << u2 << endl; } #ifdef POINTLESS cout << "\n{ isParallel for huge vectors }\n"; double multiplier = pow(2.0,400); u1 *= multiplier; u2 = u4 + 0.4 * percent * u5; u2 *= multiplier; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel" << u2 << " FALSE" << endl; problems = true; } u2 = u4 + 1.0 * percent * u5; u2 *= multiplier; if ( u1.isParallel(u2) ) { cout << "OK: " << u1 << ".isParallel" << u2 << endl; } else { cout << "?????: " << u1 << ".isParallel" << u2 << " 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. u2 = u4 + 1.1 * percent * u5; u2 *= multiplier; if ( u1.isParallel(u2) ) { cout << "?????: " << u1 << ".isParallel" << u2 << endl; problems = true; } else { cout << "OK: NOT " << u1 << ".isParallel" << u2 << endl; } #else cout << "\nSkipping huge vector test. Pointless for unit vectors." << endl; #endif // POINTLESS cout << "\n{ isOrthogonal }\n"; percent = .01; SpaceVector::setTolerance(percent); cout << " ... now default tolerance is " << SpaceVector::getTolerance() << endl; u1.set(6, -3, 5); v1 = u1; u4.set(2, 14, 6); v4 = u4; u5.set(10, -5, 10); v5 = u5; // u5 is about the size of the others, and is nearly aligned with u1. // It turns out that when you add about 1.03 epsilon*u5 to u4, you barely // miss being orthogonal to u1 by our definition. u2 = u4; if ( u1.isOrthogonal(u2,0.00000001) ) { cout << "OK: " << u1 << ".isOrthogonal" << u2 << "(almost) exactly" << endl; } else { cout << "?????: " << u1 << ".isOrthogonal" << u2 << " FALSE" << endl; problems = true; } v2 = v4 + .6 * percent * v5; u2 = v2; if ( u1.isOrthogonal(u2) ) { cout << "OK: " << u1 << ".isOrthogonal" << u2 << endl; } else { cout << "?????: " << u1 << ".isOrthogonal\n " << u4 << " + " << (.6*percent*u5) << " FALSE" << endl; problems = true; } v2 = v4 + .9 * percent * v5; u2 = v2; if ( u1.isOrthogonal(u2) ) { cout << "OK: " << u1 << ".isOrthogonal" << u2 << endl; } else { cout << "?????: " << u1 << ".isOrthogonal\n " << u4 << " + " << (.9*percent*u5) << " FALSE" << endl; problems = true; } v2 = v4 + 1.0 * percent * v5; u2 = v2; if ( u1.isOrthogonal(u2) ) { cout << "OK: " << u1 << ".isOrthogonal" << u2 << endl; } else { cout << "?????: " << u1 << ".isOrthogonal\n " << u4 << " + " << (1.0*percent*u5) << " 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; u2 = v2; if ( u1.isOrthogonal(u2) ) { cout << "?????: " << u1 << ".isOrthogonal\n " << u4 << " + " << (1.1*percent*u5) << endl; problems = true; } else { cout << "OK: NOT " << u1 << ".isOrthogonal" << u2 << endl; } #ifdef POINTLESS cout << "\n{ isOrthogonal for huge vectors }\n"; multiplier = pow(2.0,400); u1 *= multiplier; u2 = u4 + 1.0 * percent * u5; u2 *= multiplier; if ( u1.isOrthogonal(u2) ) { cout << "OK: " << u1 << ".isOrthogonal" << u2 << endl; } else { cout << "?????: " << u1 << ".isOrthogonal" << u2 << " 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. u2 = u4 + 1.1 * percent * u5; u2 *= multiplier; if ( u1.isOrthogonal(u2) ) { cout << "?????: " << u1 << ".isOrthogonal" << u2 << endl; problems = true; } else { cout << "OK: NOT " << u1 << ".isOrthogonal" << u2 << endl; } #else cout << "\nSkipping huge orthogonal vector test. Pointless for unit vectors." << endl; #endif // POINTLESS } // combining() void intrinsic() { // Simple intrinsic properties cout << "\n --- Simple intrinsic properties\n\n"; cout << " { mag() }\n"; u1.set (3, -4, -12); if (u1.mag() == 1) { cout << "OK: " << u1 << ".mag() = " << u1.mag() << endl; } else { cout << "?????" << u1 << ".mag() = " << u1.mag() << endl; problems = true; } cout << "\n{ unit() }\n"; u1.set (15, 12, -16); if (u1.unit().isNear(UnitVector(.6, .48, -.64)) ){ cout << "OK: " << u1 << ".unit() = " << u1.unit() << endl; } else { cout << "?????" << u1 << ".unit() = " << u1.unit() << endl; problems = true; } } // intrinsic() void input() { // Input cout << "\n --- Input\n\n"; cout << "\n{ Input of entire vector }\n"; cin >> u1; if (u1.isNear(UnitVector(99., -99., 95.))) { cout << "OK: input vector is " << u1 << endl; } else { cout << "????? input vector should be UnitVector(99, -99, 95) but is " << u1 << endl; problems=true; } #ifdef ZAP4 cin >> u1c; // Should fail, and does. cout << "New value of const vector u1c is " << u1c << endl; #endif // ZAP4 #ifdef ZAP4 cin >> u1c.x(); // Should fail. #endif // ZAP4 cout << "\n{ Input of individual non-cartesian component: r }\n"; cout << "\n Skipping input of individual components tests. " "Not allowed." << endl; } // input() UnitVector u ( .6, .8, 0 ); UnitVector uu ( 6, 8, 0 ); void direction() { // Properties relative to reference direction cout << "\n --- Perp and projection relative to direction\n\n"; cout << "\n{ Using the default ZHAT direction }\n"; u1.set ( 3, -4, 12 ); if ( close(u1.perp2(), 25.0/169.0) ) { cout << "OK: " << u1 << ".perp2() = " << u1.perp2() << endl; } else { cout << "????? " << u1 << ".perp2() = " << u1.perp2() << endl; problems=true; } double testres = 5./13.; double testval = u1.perp(); if ( close( testval, testres )) { cout << "OK: " << u1 << ".perp() = " << u1.perp() << endl; } else { cout << "????? " << u1 << ".perp() = " << u1.perp() << endl; problems=true; } u5 = u1.perpPart(); if ( u5.isNear(UnitVector( 3, -4, 0 ))) { cout << "OK: " << u1 << ".perpPart() = " << u1.perpPart() << endl; } else { cout << "????? " << u1 << ".perpPart() = " << u1.perpPart() << endl; problems=true; } u5 = u1.project(); if (u5.isNear(UnitVector( 0, 0, 12 ))) { cout << "OK: " << u1 << ".project() = " << u1.project() << endl; } else { cout << "????? " << u1 << ".project() = " << u1.project() << endl; problems=true; } cout << "\n{ Using a supplied reference UnitVector for direction }\n"; u1.set (2,0,4); testres = 18.56/20.; testval = u1.perp2(u); cout << "u1.perp2(u) = " << testval << " Expect " << testres << endl; if ( close(testval,testres) ) { cout << "OK: " << u1 << ".perp2" << u << " = " << u1.perp2(u) << endl; } else { cout << "????? " << u1 << ".perp2" << u << " = " << u1.perp2(u) << endl; problems=true; } testval = u1.perp(u); if ( close(testval,sqrt(testres)) ) { cout << "OK: " << u1 << ".perp" << u << " = " << u1.perp(u) << endl; } else { cout << "????? " << u1 << ".perp" << u << " = " << u1.perp(u) << endl; problems=true; } u4 = u1.perpPart(u); u5 = UnitVector ( 1.28, -.96, 4.0 ); if ( u4.isNear( u5 )) { cout << "OK: " << u1 << ".perpPart" << u << " = " << u1.perpPart(u) << endl; } else { cout << "????? " << u1 << ".perpPart" << u << " = " << u1.perpPart(u) << endl; problems=true; } u4 = u1.project(u); u5 = UnitVector ( .72, .96, 0 ); if ( u4.isNear(u5) ) { cout << "OK: " << u1 << ".project" << u << " = " << u1.project(u) << endl; } else { cout << "????? " << u1 << ".project" << u << " = " << u1.project(u) << endl; problems=true; } cout << "\n{ Using a supplied reference UnitVector }\n"; u1.set (2,0,4); testres = 18.56/20.; if ( close(u1.perp2(uu), testres) ) { cout << "OK: " << u1 << ".perp2" << uu << " = " << u1.perp2(uu) << endl; } else { cout << "????? " << u1 << ".perp2" << uu << " = " << u1.perp2(uu) << endl; } testres = sqrt(18.56/20.); testval = u1.perp(uu); if ( close( testval, testres ) ) { cout << "OK: " << u1 << ".perp" << uu << " = " << u1.perp(uu) << endl; } else { cout << "????? " << u1 << ".perp" << uu << " = " << u1.perp(uu) << endl; problems=true; } u4 = u1.perpPart(uu); u5 = UnitVector( 1.28, -.96, 4.0 ); if ( u4.isNear(u5) ) { cout << "OK: " << u1 << ".perpPart" << uu << " = " << u1.perpPart(uu) << endl; } else { cout << "????? " << u1 << ".perpPart" << uu << " = " << u1.perpPart(uu) << endl; problems=true; } u4 = u1.project(uu); u5 = UnitVector ( .72, .96, 0 ); if ( u4.isNear(u5) ) { cout << "OK: " << u1 << ".project" << uu << " = " << v1.project(uu) << endl; } else { cout << "????? " << u1 << ".project" << uu << " = " << u1.project(uu) << endl; problems=true; } } // direction() UnitVector dir ( 4, -3, 12 ); // UnitVector constructor will normalize! void anglebetween() { // Angle between vector and a direction cout << "\n --- Angle between vector and a direction\n\n"; cout << "\n{ Using the default ZHAT direction }\n"; v1.set ( 10, .35, RADIANS, -.5, RADIANS ); u1 = v1; double miss; double arads = u1.angle(); miss = arads - .35; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".angle() = " << arads << " = .25 + " << miss << endl; } else { cout << "????? " << u1 << ".angle() = " << arads << " = .25 + " << miss << endl; problems=true; } double csx; v1.set ( 10, 60, DEGREES, -PI/2, RADIANS ); u1 = v1; csx = u1.cosTheta(); miss = csx - .5; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cosTheta() = " << csx << " = .5 + " << miss << endl; } else { cout << "????? " << u1 << ".cosTheta() = " << csx << " = .5 + " << miss << endl; problems=true; } csx = u1.cos2Theta(); miss = csx - .25; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cos2Theta() = " << csx << " = .25 + " << miss << endl; } else { cout << "????? " << u1 << ".cos2Theta() = " << csx << " = .25 + " << miss << endl; problems=true; } cout << "\n{ Using a different reference direction }\n"; u1.set ( 6, 16, 2); // This is perpendicular to dir. arads = u1.angle(dir); miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".angle" << dir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << u1 << ".angle" << dir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } arads = u1.theta(dir); miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".theta" << dir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << u1 << ".theta" << dir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } csx = u1.cosTheta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cosTheta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".cosTheta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = u1.cos2Theta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cos2Theta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".cos2Theta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = u1.eta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".eta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".eta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } u1.set ( 13, -13, 26); // This is not perpendicular to dir. double rightc2 = 31.*31./(6.*169.); // This is the actual cos2 theta. arads = u1.angle(dir); miss = arads - acos ( sqrt(rightc2) ); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".angle" << dir << " = " << arads << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".angle" << dir << " = " << arads << " = answer + " << miss << endl; problems=true; } arads = u1.theta(dir); if ( eq(arads, u1.angle(dir)) ) { cout << "OK: u1.angle(dir) == u1.theta(dir)\n"; } else { cout << "????? u1.angle(dir) != u1.theta(dir)\n"; problems=true; } csx = u1.cosTheta(dir); miss = csx - sqrt(rightc2); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cosTheta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".cosTheta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = u1.cos2Theta(dir); miss = csx - rightc2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cos2Theta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".cos2Theta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = u1.eta(dir); miss = csx + log ( tan (acos(sqrt(rightc2))/2.) ); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".eta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".eta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } cout << "\n{ Eta with reference direction parallel to unit }\n"; u1 = 3 * u; if ( u1.eta(u) > 17 ) { cout << "OK: " << u1 << ".eta()" << u << " = " << u1.eta(u) << endl; } else { cout << "?????" << u1 << ".eta()" << u << " = " << u1.eta(u) << endl; problems=true; } u1 = -2 * u; if ( u1.eta(u) < -17 ) { cout << "OK: " << u1 << ".eta()" << u << " = " << u1.eta(u) << endl; } else { cout << "?????" << u1 << ".eta()" << u << " = " << u1.eta(u) << endl; problems=true; } cout << "\n{ Angle between two UnitVectors }\n"; UnitVector udir ( 4, -3, 12 ); u1.set ( 6, 16, 2); // This is perpendicular to udir. arads = u1.angle(udir); miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".angle" << udir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << u1 << ".angle" << udir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } arads = u1.theta(udir); miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".theta" << udir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << u1 << ".theta" << udir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } csx = u1.cosTheta(udir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cosTheta" << udir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".cosTheta" << udir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = u1.cos2Theta(udir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cos2Theta" << udir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".cos2Theta" << udir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = u1.eta(udir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".eta" << udir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << u1 << ".eta" << udir << " = " << csx << " = 0 + " << miss << endl; problems=true; } u1.set ( 13, -13, 26); // This is not perpendicular to udir. rightc2 = 31.*31./(6.*169.); // This is the actual cos2 theta. arads = u1.angle(udir); miss = arads - acos ( sqrt(rightc2) ); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".angle" << udir << " = " << arads << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".angle" << udir << " = " << arads << " = answer + " << miss << endl; problems=true; } arads = u1.theta(udir); if ( eq(arads, u1.angle(udir)) ) { cout << "OK: u1.angle(udir) == u1.theta(udir)\n"; } else { cout << "????? u1.angle(udir) != u1.theta(udir)\n"; problems=true; } csx = u1.cosTheta(udir); miss = csx - sqrt(rightc2); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cosTheta" << udir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".cosTheta" << udir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = u1.cos2Theta(udir); miss = csx - rightc2; if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".cos2Theta" << udir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".cos2Theta" << udir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = u1.eta(udir); miss = csx + log ( tan (acos(sqrt(rightc2))/2.) ); if (fabs(miss) < teeny ) { cout << "OK: " << u1 << ".eta" << udir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << u1 << ".eta" << udir << " = " << csx << " = answer + " << miss << endl; problems=true; } cout << "\n{ Eta with reference direction parallel to u }\n"; cout << " ... Depending on roundoffs, expect an Exception: " "PhysicsVectors-E-Infinity \n"; u1 = 3 * u; SpaceVector vvu2 = 5 * u; if ( u1.eta(vvu2) > 17 ) { // eta of 17 is effectively infinite! cout << "OK: " << u1 << ".eta()" << vvu2 << " = " << u1.eta(vvu2) << endl; } else { cout << "?????" << u1 << ".eta()" << vvu2 << " = " << u1.eta(vvu2) << endl; problems=true; } cout << " ... Depending on roundoffs, expect an Exception: " "PhysicsVectors-E-Infinity \n"; u1 = -2 * u; if ( u1.eta(vvu2) < -17 ) { cout << "OK: " << u1 << ".eta()" << vvu2 << " = " << u1.eta(vvu2) << endl; } else { cout << "?????" << u1 << ".eta()" << vvu2 << " = " << u1.eta(vvu2) << endl; problems=true; } } // anglebetween() void angledecomposition() { // Decomposition of angle between vectors cout << "\n --- Decomposition of angle between vectors\n\n"; cout << "\n{ Unit and Unit, relative to the default ZHAT direction }\n"; // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12)/13. u1.set ( 4, -3, -5 ); cout << " ... u.theta is " << u1.theta()*180/PI << " degrees; u.theta is " << u.theta()*180/PI << " degrees " << endl; //if (close (u1.polarAngle(u), -PI/4) ) { if (close (u1.polarAngle(u), PI/4) ) { // JMM definition changed? cout << "OK: polar angle (wrt. Z axis) from u1 to u is " << u1.polarAngle(u)*180/PI << " degrees" << endl; } else { cout << "?????polar angle (wrt. Z axis) from u1 to u is " << u1.polarAngle(u)*180/PI << " degrees" << endl; problems=true; } u1.set ( -8, 6, 5 ) ; cout << " ... u.phi is " << u1.phi()*180/PI << "degrees;" << " u.phi is " << u.phi()*180/PI << " degrees " << endl; if (close (u1.azimAngle(u), -PI/2) ) { cout << "OK: azim angle (wrt. Z axis) from u1 to u is " << u1.azimAngle(u)*180/PI << " degrees" << endl; } else { cout << "?????azim angle (wrt. Z axis) from u1 to u is " << u1.azimAngle(u)*180/PI << " degrees" << endl; problems=true; } cout << " { Unit and vector, relative to the default ZHAT direction }\n"; u1.set ( 4, -3, -5 ); cout << " ... u.theta is " << u1.theta()*180/PI << " degrees;" " uu.theta is " << uu.theta()*180/PI << " degrees " << endl; if (close (uu.polarAngle(u1), PI/4) ) { cout << "OK: polar angle (wrt. Z axis) from uu to u is " << uu.polarAngle(u1)*180/PI << " degrees" << endl; } else { cout << "?????polar angle (wrt. Z axis) from uu to u is " << uu.polarAngle(u1)*180/PI << " degrees" << endl; problems=true; } u1.set ( -8, 6, 5 ) ; cout << " ... u.phi is " << u1.phi()*180/PI << " degrees; uu.phi is " << uu.phi()*180/PI << " degrees " << endl; if (close (u1.azimAngle(uu), -PI/2) ) { cout << "OK: azim angle (wrt. Z axis) from u to uu is " << u1.azimAngle(uu)*180/PI << " degrees" << endl; } else { cout << "?????azim angle (wrt. Z axis) from u to uu is " << u1.azimAngle(uu)*180/PI << " degrees" << endl; problems=true; } cout << "\n{ U U Relative to a different reference direction }\n"; // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12)/13. u1.set ( .7, .8, .1 ); cout << " * * * polarAngle of " << u1 << " to " << u << "\n with respect to" << 3.2*dir << " is " << u1.polarAngle(u, dir) *180/PI << " degrees " << endl; cout << " * * * azimAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.azimAngle(u, dir) *180/PI << " degrees " << endl; if ( close(u1.polarAngle(u, dir),u1.polarAngle(u, 3.2*dir)) ) { cout << "OK: polar angle with respect to UnitVector reference matches\n"; } else { cout << "????? polar angle with respect to UnitVector reference matches\n"; problems=true; } if ( close(u1.azimAngle(u, dir),u1.azimAngle(u, 3.2*dir)) ) { cout << "OK: azim angle with respect to UnitVector reference matches\n"; } else { cout << "????? azim angle with respect to UnitVector reference matches\n"; problems=true; } u1.set ( -.6, -.8, 0 ); // -u if ( close(u1.polarAngle(u, dir) , 0) ) { cout << "OK: polarAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.polarAngle(u, dir)* 180/PI << " degrees " << endl; } else { cout << "????? polarAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.polarAngle(u, dir)* 180/PI << " degrees " << endl; problems=true; } cout << " Explanation: Angle between " << u1 << " and " << dir << " is " << u1.angle(dir) << "\n Angle between " << u << " and " << dir << " is " << u.angle(dir) << "\n"; if ( close(fabs(u1.azimAngle(u, dir)) , PI ) ) { cout << "OK: azimAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.azimAngle(u, dir)* 180/PI << " degrees" << endl; } else { cout << "????? azimAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.azimAngle(u, dir)* 180/PI << " degrees" << endl; problems=true; } u1.set ( -.6, -.8, .1 ); // Fairly close to -u cout << " * * * polarAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.polarAngle(u, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << u1 << " to " << u << "\n with respect to" << dir << " is " << u1.azimAngle(u, dir) *180/PI << " degrees" << endl; cout << "\n{ U V Relative to a different reference direction }\n"; // Reminder -- uu is (6, 8, 0) and dir is (4, -3, 12)/13. u1.set ( .7, .8, .1 ); cout << " * * * polarAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.polarAngle(uu, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.azimAngle(uu, dir) *180/PI << " degrees" << endl; u1.set ( -.6, -.8, 0 ); // Fairly close to uu if ( close(u1.polarAngle(uu, dir) , 0) ) { cout << "OK: polarAngle of " << u1 << " to " << uu << " \n with respect to" << dir << " is " << u1.polarAngle(uu, dir)* 180/PI << " degrees" << endl; } else { cout << "????? polarAngle of " << u1 << " to " << uu << " \n with respect to" << dir << " is " << u1.polarAngle(uu, dir)* 180/PI << " degrees" << endl; problems=true; } if ( close(fabs(u1.azimAngle(uu, dir)) , PI ) ) { cout << "OK: azimAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.azimAngle(uu, dir)* 180/PI << " degrees" << endl; } else { cout << "????? azimAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.azimAngle(uu, dir)* 180/PI << " degrees" << endl; problems=true; } u1.set ( -.6, -.8, .1 ); cout << " * * * polarAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.polarAngle(uu, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << u1 << " to " << uu << "\n with respect to" << dir << " is " << u1.azimAngle(uu, dir)*180/PI << " degrees" << endl; } // angledecomposition() void rotations() { // ROTATIONS cout << "\n --- Rotations\n\n"; cout << "\n{ 90 degree rotations about X Y or Z axis }\n"; u1.set ( 9, 5, -3 ); SpaceVector::setTolerance ( .0001 ); u2 = u1; u2.rotateX( PI/2 ); u3.set ( 9, 3, 5 ); if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateX(90 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateX(90 degrees) = " << u2 << endl; problems=true; } u2 = u3; u2.rotateX( -PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u3 << "rotateX( -PI/2 ) = " << u2 << endl; } else { cout << "????? " << u3 << "rotateX( -PI/2 ) = " << u2 << endl; problems=true; } u2 = u1; u2.rotateY( PI/2 ); u3.set ( -3, 5, -9 ); if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateY(90 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateY(90 degrees) = " << u2 << endl; problems=true; } u2 = u3; u2.rotateY( -PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u3 << "rotateY( -PI/2 ) = " << u2 << endl; } else { cout << "????? " << u3 << "rotateY( -PI/2 ) = " << u2 << endl; problems=true; } u2 = u1; u2.rotateZ( PI/2 ); u3.set ( -5, 9, -3 ); if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateZ(90 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateZ(90 degrees) = " << u2 << endl; problems=true; } u2 = u3; u2.rotateZ( -PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u3 << "rotateZ( -PI/2 ) = " << u2 << endl; } else { cout << "????? " << u3 << "rotateZ( -PI/2 ) = " << u2 << endl; problems=true; } cout << "\n{ Pair of 90 degree rotations to permute x y z }\n"; u1.set ( 9, 5, -3 ); u3.set ( -3, 9, 5 ); u2 = u1; u2.rotateX ( PI/2 ); u2.rotateZ ( PI/2 ); if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "becomes " << u2 << endl; } else { cout << "????? " << u1 << "becomes " << u2 << endl; problems=true; } u3 = u2; u2.rotateZ ( -PI/2 ); u2.rotateX ( -PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u3 << "becomes " << u2 << endl; } else { cout << "????? " << u3 << "becomes " << u2 << endl; problems=true; } cout << "\n{ Axis rotations by general angles }\n"; u1.set ( 1, 10*sqrt(3.), 10 ); u3.set ( 1, 0, 20 ); u2 = u1; u2.rotateX ( PI/3 ); if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateX(60 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateX(60 degrees) = " << u2 << endl; problems=true; } u4 = rotationXOf (u2, 2*PI/3 ); u2.rotateX ( 2*PI/3 ); if (u4==u2) { cout << "OK: uprime = rotationXOf(u,phi) matches u.rotateX(phi)\n"; } else { cout << "????? uprime = rotationXOf(u,phi) mismatches u.rotateX(phi)\n"; problems=true; } v1 = u1; v1.setX( v1.x() * -1 ); u1 = v1; if ( u2.isNear(-u1) ) { cout << "OK: rotate another 120 degrees gets back to " << u2 << endl; } else { cout << "????? rotate another 120 degrees gets back to " << u2 << endl; problems=true; } u1.set ( 10, 1, 10*sqrt(3.) ); u3.set ( 20, 1, 0 ); u2 = u1; u4 = rotationYOf (u2, PI/3 ); u2.rotateY ( PI/3 ); if (u4==u2) { cout << "OK: uprime = rotationYOf(u,phi) matches u.rotateY(phi)\n"; } else { cout << "????? uprime = rotationYOf(u,phi) mismatches u.rotateY(phi)\n"; problems=true; } if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateY(60 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateY(60 degrees) = " << u2 << endl; problems=true; } u2.rotateY ( 2*PI/3 ); v1 = u1; v1.setY( v1.y() * -1 ); u1 = v1; if ( u2.isNear(-u1) ) { cout << "OK: rotate another 120 degrees gets back to " << u2 << endl; } else { cout << "????? rotate another 120 degrees gets back to " << u2 << endl; problems=true; } u1.set ( 10*sqrt(3.), 10, 1 ); u3.set ( 0, 20, 1 ); u2 = u1; u4 = rotationZOf (u2, PI/3 ); u2.rotateZ ( PI/3 ); if (u4==u2) { cout << "OK: uprime = rotationZOf(u,phi) matches u.rotateZ(phi)\n"; } else { cout << "????? uprime = rotationZOf(u,phi) mismatches u.rotateZ(phi)\n"; problems=true; } if ( u2.isNear(u3) ) { cout << "OK: " << u1 << "rotateZ(60 degrees) = " << u2 << endl; } else { cout << "????? " << u1 << "rotateZ(60 degrees) = " << u2 << endl; problems=true; } u2.rotateZ ( 2*PI/3 ); v1 = u1; v1.setZ( v1.z() * -1 ); u1 = v1; if ( u2.isNear(-u1) ) { cout << "OK: rotate another 120 degrees gets back to " << u2 << endl; } else { cout << "????? rotate another 120 degrees gets back to " << u2 << endl; problems=true; } cout << "\n{ Verifying the phi++ rule }\n"; u1.set ( -13, 2, -11 ); u2 = u1; u2.rotateZ ( 2.4 ); u2.setPhi( u2.phi() - 2.4 ); if ( u2.isNear(u1) ) { cout << "OK: " << u1 << "rotateZ (2.4) phi -= 2.4 = " << u2 << endl; } else { cout << "????? " << u1 << "rotateZ (2.4) phi -= 2.4 = " << u2 << endl; problems=true; } u1.set ( 5, 21, -1 ); u2 = u1; u2.rotateX ( -2.9 ); u2.rotateZ ( -PI/2 ); u2.rotateX ( -PI/2 ); u2.setPhi( u2.phi() + 2.9 ); u2.rotateX ( PI/2 ); u2.rotateZ ( PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u1 << "rotateX (-2.9) permute phi += 2.9 permute = " << u2 << endl; } else { cout << "????? " << u1 << "rotateX (-2.9) permute phi += 2.9 permute = " << u2 << endl; problems=true; } u1.set ( 15, -2, -11 ); u2 = u1; u2.rotateY ( 1.9 ); u2.rotateX ( PI/2 ); u2.rotateZ ( PI/2 ); u2.setPhi( u2.phi() - 1.9 ); u2.rotateZ ( -PI/2 ); u2.rotateX ( -PI/2 ); if ( u2.isNear(u1) ) { cout << "OK: " << u1 << "rotateY (1.9) permute phi -= 1.9 permute = " << u2 << endl; } else { cout << "????? " << u1 << "rotateY (1.9) permute phi -= 1.9 permute = " << u2 << endl; problems=true; } cout << "\n{ Generic rotations using XYZ axes }\n"; u1.set (8, -5, 3); SpaceVector xv (4,0,0); UnitVector xu (1,0,0); u3 = u2 = u1; u2.rotate ( xu, 50*PI/180 ); u3.rotateX ( 50*PI/180 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( xv, -.4 ); u3.rotateX( -.4 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( -xv, .9 ); u3.rotateX( -.9 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } SpaceVector yv (0,5,0); UnitVector yu (0,-1,0); u3 = u2 = u1; u2.rotate ( yu, 50*PI/180 ); u3.rotateY ( -50*PI/180 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( yv, -.4 ); u3.rotateY( -.4 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( -yv, .9 ); u3.rotateY( -.9 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } SpaceVector zv (0,0,6); UnitVector zu (0,0,1); u3 = u2 = u1; u2.rotate ( zu, 50*PI/180 ); u3.rotateZ ( 50*PI/180 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( zv, -.4 ); u3.rotateZ( -.4 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } u3 = u2; u2.rotate ( -zv, .9 ); u3.rotateZ( -.9 ); if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } cout << "\n{ Generic rotations about arbitrary axes }\n"; // Laboriously hand-calculated : uu.set ( 1, 1, 0 ); u1.set ( 3, 4, 0 ); u3 = u2 = u1; u2.rotate ( uu, PI/4 ); cout << " " << u1 << "rotated 45 degrees about " << uu << " is" << u2 << endl; u3.rotate ( uu, -3*PI/4 ); cout << " "<< u1 << "rotated -135 degrees about " << uu << " is" << u3 << endl; v3 = u3; double temp = v3.x(); v3.setX( v3.y() ); v3.setY( temp ); v3.setZ( v3.z() * -1 ); u3 = v3; if ( u2.isNear(u3) ) { cout << "OK: " << u2 << " ~ " << u3 << endl; } else { cout << "????? " << u2 << " not near " << u3 << endl; problems=true; } // Using formula described below: double phi; u1.set(3,4,1); u2.set(1,1,0); phi = PI/4; u3.set( 3.6464, 3.3536, 1.2071 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(3,4,1); u2.set(2,2,0); phi = PI/2; u3.set( 4.2071, 2.7929, 0.7071 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(3,4,1); u2.set(1,1,0); phi = 175*PI/180; u3.set( 4.0597, 2.9403, -0.9346 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; #ifdef INTENTIONAL_ERROR // Intentionally wrong to check our checker: u1.set(3,4,1); u2.set(1,1,0); phi = 175*PI/180; u3.set( 4.0397, 2.9403, -0.9346 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; #endif // INTENTIONAL_ERROR u1.set(2,-1,3); u2.set(4,2,-5); phi = 35*PI/180; u3.set( 1.5791, -2.7726, 1.9543 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,-1,3); u2.set(4,2,-5); phi = -10*PI/180; u3.set( 1.9316, -0.4214, 3.1767 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,-1,3); u2.set(4,2,-5); phi = -100*PI/180; u3.set( -1.4330, 2.9340, 1.8272 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(0,1,0); phi = 20*PI/180; u3.set( 3.2475, 3, 3.0747); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(0,3,0); phi = PI/2; u3.set(4, 3, -2); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(0,1,0); phi = PI; u3.set( -2, 3, -4 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(0,-1,0); phi = 20*PI/180; u3.set( 0.5113, 3, 4.4428 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(0,-1,0); phi = -20*PI/180; u3.set( 3.2475, 3, 3.0747); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(1,3,4); u2.set(0,.01,1); phi = 1*PI/180; u3.set( 0.9482, 3.0170, 3.9998 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(1,3,4); u2.set(0,.01,1); phi = -10*PI/180; u3.set( 1.4988, 2.7814, 4.0022 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(1, -1, 1); phi = -130*PI/180; u3.set( 3.4531, -2.6866, -3.1397 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(1, -1, 1); phi = -80*PI/180; u3.set( 5.1537, 0.8318, -1.3220 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(1, -1, 1); phi = -PI/6; u3.set( 3.8868, 3.0415, 2.1547 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(1, -1, 1); phi = 20*PI/180; u3.set( 0.5574, 2.3638, 4.8064 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; u1.set(2,3,4); u2.set(1, -1, 1); phi = 70*PI/180; u3.set( -2.4557, -0.7170, 4.7387 ); if ( checkAxialRotation (u1, u2, phi, u3 ) ) problems = true; } // rotations()