// svtest.cc -- test of features of SpaceVector class that depend on UnitVector //----------------------------------------------------------------------------- #include "ZMutility/ZMenvironment.h" #include "PhysicsVectors/SpaceVector.h" #include "PhysicsVectors/AxisAngle.h" #include "PhysicsVectors/EulerAngles.h" #include "PhysicsVectors/UnitVector.h" #include "CLHEP/Vector/ZMxpv.h" ZM_USING_NAMESPACE( zmpv ) #include "ZMutility/cmath" using CMATH_NAMESPACE::fabs; using CMATH_NAMESPACE::log; using CMATH_NAMESPACE::pow; using CMATH_NAMESPACE::tan; #include "ZMutility/iostream" #include "ZMutility/iomanip" #include "ZMutility/fstream" //#include "CLHEP/Vector/ZMxpv.h" USING( std::cin ) USING( std::cout ) USING( std::endl ) #include "ZMutility/PI.h" const Float8 PI = ZMpi; bool problems( false ); SpaceVector vans; void direction(); void anglebetween(); void angledecomposition(); void rotations2(); bool checkAxialRotation ( const SpaceVector &x, const SpaceVector &d, const Float8 &phi, const SpaceVector& result ); SpaceVector goldstein107 (const SpaceVector& v, Float8 phi, Float8 theta, Float8 psi); bool close ( double a, double b ) { return ( fabs(a-b) < 1.0E-10 ); } bool eq ( double a, double b ) { return ( fabs(a-b) < 1.0E-14 ); } const double teeny = 1.0E-14; int main() { // For purposes of this test, it is safe and necessary to ignore ZMexceptions. // This test exercises some of them but we want to continue beyond each one. ZMxPhysicsVectors::setHandler( ZMexIgnoreAlways() ); cout << "\n\n"; cout << " ***********\n"; cout << " SpaceVector\n"; cout << " ***********\n"; direction(); anglebetween(); angledecomposition(); rotations2(); 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 svtestu main() */ bool checkAxialRotation ( const SpaceVector &x, const SpaceVector &d, const double& phi, const SpaceVector& result ) { // Rotate x by phi around d; return true if result does not match answer // within .0001 tolerance; print out appropriate messages. SpaceVector xx (x); xx.rotate(d, phi); if ( xx.isNear(result, .0001) ) { cout << "OK: " << x << "Rotated " << phi << " about" << d << "\n = " << xx << " ~ " << result << endl; } else { cout << "?????" << x << "Rotated " << phi << " about" << d << "\n = " << xx << " not near " << result << endl; return true; } // Now form a new SpaceVector as xx rotated about the axis by -phi; // again check that the result is proper (matches x) very tightly SpaceVector yy = rotationOf(xx, d, -phi); if ( yy.isNear(x, .000001) ) { cout << "OK: " << xx << "Rotated " << -phi << " about" << d << "\n = " << yy << " ~ " << x << endl; } else { cout << "?????" << xx << "Rotated " << -phi << " about" << d << "\n = " << yy << " not near " << x << endl; return true; } return false; // How did we get those results that we have supplied? // // We rotate the vector in the following manner: // (1) From a set of three coordinate-like vectors with d used as the third: // a - choose V1 = d cross x // b - choose V2 = d cross V1 // c - the third is just d // (2) Decompose x into u V1 + v V2 + w d // You know a priori that u = 0 // Unless d is purely in the Z direction, the two equations dealing with // the x and y directions suffice: // v = (x.x*d.y - d.x*x.y) / ( V2.x*d.y - d.x*V2.y ) and similar for w. // (3) express x in terms of V2-hat instead of V2 // (4) Rotate in V2-hat, V1hat space by phi // (5) rescale back to V1, V2 // (6) transform back using V1, V2, and d to the x,y,z coordinates. // Steps 3--6 combine: // x' = - v sin(phi) (|V2| \ |V1|) V1 + v cos(phi) V2 + w d // // That this result matches the formula used in SpaceVector.cc is a highly // non-trivial test. } /* checkAxialRotation */ SpaceVector goldstein107 (const SpaceVector& v, Float8 phi, Float8 theta, Float8 psi) { // Perform an Euler angles rotation of v by laboriously folloing the // illustration on page 107 of Goldstein Classical Mechanincs. Then // check into the booby hatch! // Note that our rotations about an axis are active, but Goldstein uses // (passive) frame transformations, which is the only good way to define // Euler angles. // Step 1 (diagram a): // // Rotate the X axis by phi around the Z axis to form a new axis // which Goldstein calls xi. In this new frame, v is the old v rotated // by -phi about Z to form v'. // #define TRACE_GOLD UnitVector xi = rotationOf (X_HAT, Z_HAT, phi); UnitVector eta = rotationOf (Y_HAT, Z_HAT, phi); UnitVector zeta = rotationOf (Z_HAT, Z_HAT, phi); SpaceVector v_prime = rotationOf (v, Z_HAT, -phi); #ifdef TRACE_GOLD cout << " goldstein107: V is " << v << endl; cout << " new frame is xi, eta, zeta =\n" << " " << xi << endl << " " << eta << endl << " " << zeta << endl; cout << " in that frame V has components\n" << " " << v.dot(xi) << " , " << v.dot(eta) <<" , " << v.dot(zeta) << endl; cout << " and V' is " << v_prime << endl; #endif // Step 2 (diagram b): // // Rotate the Z or zeta axis by theta around that xi axis to form a new axis // which Goldstein calls zeta_prime. Since in that frame a rotation about the // xi axis has the form we are familiar with as a RotationX, we rotate about // v about X_HAT to get the effect of this rotation. Thus in the xi_prime, // eta_prime, zeta_prime coordinate frame, v'' is v', X-rotated by -theta. UnitVector zeta_prime = rotationOf (zeta, xi, theta); UnitVector eta_prime = rotationOf (eta, xi, theta); UnitVector xi_prime = rotationOf (xi, xi, theta); SpaceVector v_double_prime = rotationOf (v_prime, X_HAT, -theta); #ifdef TRACE_GOLD cout << " Step 2" << endl; cout << " new frame is xi', eta', zeta' =\n" << " " << xi_prime << endl << " " << eta_prime << endl << " " << zeta_prime << endl; cout << " in that frame V has components\n" << " " << v.dot(xi_prime) << " , " << v.dot(eta_prime) <<" , " << v.dot(zeta_prime) << endl; cout << " and V'' is " << v_double_prime << endl; #endif // Step 3 (diagram c): // // Rotate the xi (xi-prime) axis by psi around the new zeta-prime axis. // Once again, the effect on V is a rotation using the simple RotateZ form // in those primed axes. This rotation produces the new coordinate frame // labeled x-prime, (the rotated xi-prime) y-prime (the rotated eta-prime) // and z-prime (same as zeta-prime or zeta). // We find our final v''' by rotating v'' by -psi around the Z axis. UnitVector x_prime = rotationOf (xi_prime, zeta_prime, psi); UnitVector y_prime = rotationOf (eta_prime, zeta_prime, psi); UnitVector z_prime = rotationOf (zeta_prime, zeta_prime, psi); SpaceVector result = rotationOf (v_double_prime, Z_HAT, -psi ); #ifdef TRACE_GOLD cout << " Step 3" << endl; cout << " new frame is x', y', z' =\n" << " " << x_prime << endl << " " << y_prime << endl << " " << z_prime << endl; cout << " in that frame V has components\n" << " " << v.dot(x_prime) << " , " << v.dot(y_prime) << " , " << v.dot(z_prime) << endl; cout << " and result is " << result << endl; #endif return result; } /* goldstein107 */ SpaceVector v1 ( 1., 2., 3. ); const SpaceVector v1c ( 2, 3., 4. ); SpaceVector v2 ( 10, 30, DEGREES, 60, DEGREES ); const SpaceVector v2c ( 20, 30, DEGREES, PI/3, RADIANS ); SpaceVector v3 ( 20, 90, DEGREES, -15 ); UnitVector u ( .6, .8, 0 ); SpaceVector vu ( 6, 8, 0 ); void direction() { // Properties relative to reference direction cout << "\n --- Perp and projection relative to direction\n\n"; cout << " { Using a supplied reference UnitVector for direction }\n"; v1.set (2,0,4); if ( eq(v1.perp2(u),18.56) ) { cout << "OK: " << v1 << ".perp2" << u << " = " << v1.perp2(u) << endl; } else { cout << "????? " << v1 << ".perp2" << u << " = " << v1.perp2(u) << endl; problems=true; } if ( eq(v1.perp(u),sqrt(18.56)) ) { cout << "OK: " << v1 << ".perp" << u << " = " << v1.perp(u) << endl; } else { cout << "????? " << v1 << ".perp" << u << " = " << v1.perp(u) << endl; problems=true; } if ( v1.perpPart(u).isNear(SpaceVector ( 1.28, -.96, 4.0 )) ) { cout << "OK: " << v1 << ".perpPart" << u << " = " << v1.perpPart(u) << endl; } else { cout << "????? " << v1 << ".perpPart" << u << " = " << v1.perpPart(u) << endl; problems=true; } if ( v1.project(u).isNear(SpaceVector ( .72, .96, 0 )) ) { cout << "OK: " << v1 << ".project" << u << " = " << v1.project(u) << endl; } else { cout << "????? " << v1 << ".project" << u << " = " << v1.project(u) << 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 << " { Using a different reference direction }\n"; v1.set ( 6, 16, 2); // This is perpendicular to dir. double arads = v1.angle(dir); double miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".angle" << dir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << v1 << ".angle" << dir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } arads = v1.theta(dir); miss = arads - PI/2; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".theta" << dir << " = " << arads << " = PI/2 + " << miss << endl; } else { cout << "????? " << v1 << ".theta" << dir << " = " << arads << " = PI/2 + " << miss << endl; problems=true; } double csx = v1.cosTheta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".cosTheta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << v1 << ".cosTheta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = v1.cos2Theta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".cos2Theta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << v1 << ".cos2Theta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } csx = v1.eta(dir); miss = csx; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".eta" << dir << " = " << csx << " = 0 + " << miss << endl; } else { cout << "????? " << v1 << ".eta" << dir << " = " << csx << " = 0 + " << miss << endl; problems=true; } v1.set ( 13, -13, 26); // This is not perpendicular to dir. double rightc2 = 31.*31./(6.*169.); // This is the actual cos2 theta. arads = v1.angle(dir); miss = arads - acos ( sqrt(rightc2) ); if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".angle" << dir << " = " << arads << " = answer + " << miss << endl; } else { cout << "????? " << v1 << ".angle" << dir << " = " << arads << " = answer + " << miss << endl; problems=true; } arads = v1.theta(dir); if ( eq(arads, v1.angle(dir)) ) { cout << "OK: v1.angle(dir) == v1.theta(dir)\n"; } else { cout << "????? v1.angle(dir) != v1.theta(dir)\n"; problems=true; } csx = v1.cosTheta(dir); miss = csx - sqrt(rightc2); if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".cosTheta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << v1 << ".cosTheta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = v1.cos2Theta(dir); miss = csx - rightc2; if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".cos2Theta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << v1 << ".cos2Theta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } csx = v1.eta(dir); miss = csx + log ( tan (acos(sqrt(rightc2))/2.) ); if (fabs(miss) < teeny ) { cout << "OK: " << v1 << ".eta" << dir << " = " << csx << " = answer + " << miss << endl; } else { cout << "????? " << v1 << ".eta" << dir << " = " << csx << " = answer + " << miss << endl; problems=true; } cout << " { Eta with reference direction parallel to unit }\n"; v1 = 3 * u; if ( v1.eta(u) > 17 ) { cout << "OK: " << v1 << ".eta()" << u << " = " << v1.eta(u) << endl; } else { cout << "?????" << v1 << ".eta()" << u << " = " << v1.eta(u) << endl; problems=true; } v1 = -2 * u; if ( v1.eta(u) < -17 ) { cout << "OK: " << v1 << ".eta()" << u << " = " << v1.eta(u) << endl; } else { cout << "?????" << v1 << ".eta()" << u << " = " << v1.eta(u) << endl; problems=true; } cout << " { Eta with reference direction parallel to v }\n"; v1 = 3 * u; v2 = 5 * u; if ( v1.eta(v2) > 17 ) { cout << "OK: " << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl; } else { cout << "?????" << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl; problems=true; } v1 = -2 * u; if ( v1.eta(v2) < -17 ) { cout << "OK: " << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl; } else { cout << "?????" << v1 << ".eta()" << v2 << " = " << v1.eta(v2) << endl; problems=true; } } // anglebetween() void angledecomposition() { // Decomposition of angle between vectors cout << "\n --- Decomposition of angle between vectors\n\n"; cout << " { Vector and Unit, relative to the default ZHAT direction }\n"; // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12)/13. v1.set ( 4, -3, -5 ); cout << " ... v.theta is " << v1.theta()*180/PI << " degrees; u.theta is " << u.theta()*180/PI << " degrees " << endl; //if (close (v1.polarAngle(u), -PI/4) ) { if (close (v1.polarAngle(u), PI/4) ) { // JMM definition changed? cout << "OK: polar angle (wrt. Z axis) from v to u is " << v1.polarAngle(u)*180/PI << " degrees" << endl; } else { cout << "?????polar angle (wrt. Z axis) from v to u is " << v1.polarAngle(u)*180/PI << " degrees" << endl; problems=true; } v1.set ( -8, 6, 5 ) ; cout << " ... v.phi is " << v1.phi()*180/PI << "degrees;" << " u.phi is " << u.phi()*180/PI << " degrees " << endl; if (close (v1.azimAngle(u), -PI/2) ) { cout << "OK: azim angle (wrt. Z axis) from v to u is " << v1.azimAngle(u)*180/PI << " degrees" << endl; } else { cout << "?????azim angle (wrt. Z axis) from v to u is " << v1.azimAngle(u)*180/PI << " degrees" << endl; problems=true; } cout << "\n { V U Relative to a different reference direction }\n"; // Reminder -- u is (.6, .8, 0) and dir is (4, -3, 12)/13. v1.set ( .7, .8, .1 ); cout << " * * * polarAngle of " << v1 << "to" << u << "\n with respect to" << 3.2*dir << "is " << v1.polarAngle(u, dir) *180/PI << " degrees " << endl; cout << " * * * azimAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.azimAngle(u, dir) *180/PI << " degrees " << endl; if ( close(v1.polarAngle(u, dir),v1.polarAngle(u, 3.2*dir)) ) { cout << "OK: polar angle with respect to SpaceVector reference matches\n"; } else { cout << "????? polar angle with respect to SpaceVector reference matches\n"; problems=true; } if ( close(v1.azimAngle(u, dir),v1.azimAngle(u, 3.2*dir)) ) { cout << "OK: azim angle with respect to SpaceVector reference matches\n"; } else { cout << "????? azim angle with respect to SpaceVector reference matches\n"; problems=true; } v1.set ( -.6, -.8, 0 ); // -u if ( close(v1.polarAngle(u, dir) , 0) ) { cout << "OK: polarAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.polarAngle(u, dir)* 180/PI << " degrees " << endl; } else { cout << "????? polarAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.polarAngle(u, dir)* 180/PI << " degrees " << endl; problems=true; } cout << " Explanation: Angle between " << v1 << " and " << dir << " is " << v1.angle(dir) << "\n Angle between " << u << " and " << dir << " is " << u.angle(dir) << "\n"; if ( close(fabs(v1.azimAngle(u, dir)) , PI ) ) { cout << "OK: azimAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.azimAngle(u, dir)* 180/PI << " degrees" << endl; } else { cout << "????? azimAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.azimAngle(u, dir)* 180/PI << " degrees" << endl; problems=true; } v1.set ( -.6, -.8, .1 ); // Fairly close to -u cout << " * * * polarAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.polarAngle(u, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << v1 << "to" << u << "\n with respect to" << dir << "is " << v1.azimAngle(u, dir) *180/PI << " degrees" << endl; cout << "\n { V V Relative to a different reference direction }\n"; // Reminder -- vu is (6, 8, 0) and dir is (4, -3, 12)/13. v1.set ( .7, .8, .1 ); cout << " * * * polarAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.polarAngle(vu, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.azimAngle(vu, dir) *180/PI << " degrees" << endl; v1.set ( -.6, -.8, 0 ); // Fairly close to vu if ( close(v1.polarAngle(vu, dir) , 0) ) { cout << "OK: polarAngle of " << v1 << "to" << vu << " \n with respect to" << dir << "is " << v1.polarAngle(vu, dir)* 180/PI << " degrees" << endl; } else { cout << "????? polarAngle of " << v1 << "to" << vu << " \n with respect to" << dir << "is " << v1.polarAngle(vu, dir)* 180/PI << " degrees" << endl; problems=true; } if ( close(fabs(v1.azimAngle(vu, dir)) , PI ) ) { cout << "OK: azimAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.azimAngle(vu, dir)* 180/PI << " degrees" << endl; } else { cout << "????? azimAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.azimAngle(vu, dir)* 180/PI << " degrees" << endl; problems=true; } v1.set ( -.6, -.8, .1 ); cout << " * * * polarAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.polarAngle(vu, dir) *180/PI << " degrees" << endl; cout << " * * * azimAngle of " << v1 << "to" << vu << "\n with respect to" << dir << "is " << v1.azimAngle(vu, dir)*180/PI << " degrees" << endl; } // angledecomposition() void rotations2() { cout << " { Generic rotations using XYZ axes }\n"; v1.set (8, -5, 3); SpaceVector xv (4,0,0); UnitVector xu (1,0,0); v3 = v2 = v1; v2.rotate ( xu, 50*PI/180 ); v3.rotateX ( 50*PI/180 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( xv, -.4 ); v3.rotateX( -.4 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( -xv, .9 ); v3.rotateX( -.9 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } SpaceVector yv (0,5,0); UnitVector yu (0,-1,0); v3 = v2 = v1; v2.rotate ( yu, 50*PI/180 ); v3.rotateY ( -50*PI/180 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( yv, -.4 ); v3.rotateY( -.4 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( -yv, .9 ); v3.rotateY( -.9 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } SpaceVector zv (0,0,6); UnitVector zu (0,0,1); v3 = v2 = v1; v2.rotate ( zu, 50*PI/180 ); v3.rotateZ ( 50*PI/180 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( zv, -.4 ); v3.rotateZ( -.4 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } v3 = v2; v2.rotate ( -zv, .9 ); v3.rotateZ( -.9 ); if ( v2.isNear(v3) ) { cout << "OK: " << v2 << " ~ " << v3 << endl; } else { cout << "????? " << v2 << " not near " << v3 << endl; problems=true; } } // rotations2()