Changeset 42a101

Timestamp:

May 28, 2010, 8:25:27 AM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

ccf826

Parents:

643e76

Message:

Removed RotateVector() function in favor of Line::rotateVector() method.

Line::rotateVector() is able to rotate any vector around any line in space.

Location:

src

Files:

• Legacy/oldmenu.cpp (modified) (3 diffs)
• Line.cpp (modified) (1 diff)
• Line.hpp (modified) (1 diff)
• Plane.cpp (modified) (1 diff)
• Space.cpp (modified) (1 diff)
• Space.hpp (modified) (1 diff)
• unittests/LineUnittest.cpp (modified) (1 diff)
• unittests/LineUnittest.hpp (modified) (2 diffs)
• unittests/vectorunittest.cpp (modified) (1 diff)
• unittests/vectorunittest.hpp (modified) (1 diff)
• vector.cpp (modified) (2 diffs)
• vector.hpp (modified) (1 diff)
• vector_ops.cpp (modified) (1 diff)
• vector_ops.hpp (modified) (1 diff)

src/Legacy/oldmenu.cpp

@@ -26 +26 @@
 #include "vector_ops.hpp"
 #include "Plane.hpp"
+#include "Line.hpp"
 
 #include "UIElements/UIFactory.hpp"
@@ -184 +185 @@
           // rotate vector around first angle
           first->x = x;
-          first->x = RotateVector(first->x,z,b - M_PI);
+          first->x = Line(zeroVec,z).rotateVector(first->x,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
           // remove the projection onto the rotation plane of the second angle
@@ -200 +201 @@
           // rotate another vector around second angle
           n = y;
-          n = RotateVector(n,x,c - M_PI);
+          n = Line(zeroVec,x).rotateVector(n,c - M_PI);
           Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 

src/Line.cpp

@@ -175 +175 @@
 }
 
+/** Rotates the vector by an angle of \a alpha around this line.
+ * \param rhs Vector to rotate
+ * \param alpha rotation angle in radian
+ */
+Vector Line::rotateVector(const Vector &rhs, double alpha) const{
+  Vector helper = rhs;
+
+  // translate the coordinate system so that the line goes through (0,0,0)
+  helper -= *origin;
+
+  // partition the vector into a part that gets rotated and a part that lies along the line
+  pair<Vector,Vector> parts = helper.partition(*direction);
+
+  // we just keep anything that is along the axis
+  Vector res = parts.first;
+
+  // the rest has to be rotated
+  Vector a = parts.second;
+  // we only have to do the rest, if we actually could partition the vector
+  if(!a.IsZero()){
+    // construct a vector that is orthogonal to a and direction and has length |a|
+    Vector y = a;
+    // direction is normalized, so the result has length |a|
+    y.VectorProduct(*direction);
+
+    res += cos(alpha) * a + sin(alpha) * y;
+  }
+
+  // translate the coordinate system back
+  res += *origin;
+  return res;
+}
+
 Line makeLineThrough(const Vector &x1, const Vector &x2){
   if(x1==x2){

src/Line.hpp

@@ -33 +33 @@
   Vector getIntersection(const Line& otherLine) const;
 
+  Vector rotateVector(const Vector &rhs, double alpha) const;
+
 private:
   std::auto_ptr<Vector> origin;

src/Plane.cpp

@@ -42 +42 @@
 /**
  * Constructs a plane from two direction vectors and a offset.
- * If no offset is given a plane through origin is assumed
  */
 Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException) :

src/Space.cpp

@@ -19 +19 @@
 }
 
+Vector Space::getVectorToPoint(const Vector &origin) const{
+  Vector support = getClosestPoint(origin);
+  return support-origin;
+}
+
 bool Space::isContained(const Vector &point) const{
   return (distance(point)) < MYEPSILON;

src/Space.hpp

@@ -19 +19 @@
   virtual double distance(const Vector &point) const=0;
   virtual Vector getClosestPoint(const Vector &point) const=0;
+  virtual Vector getVectorToPoint(const Vector &point) const;
   virtual bool isContained(const Vector &point) const;
   virtual bool hasZero() const;

src/unittests/LineUnittest.cpp

@@ -221 +221 @@
   CPPUNIT_ASSERT_THROW(la3->getIntersection(*lp1),SkewException);
 }
+
+void LineUnittest::rotationTest(){
+  Vector fixture;
+
+  // rotate zero Vector along the axes lines by various degrees
+  fixture = la1->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la2->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la3->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+
+  fixture = la1->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la2->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la3->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+
+  // rotate vectors on the axis around their lines
+  fixture = la1->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = la1->rotateVector(e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->rotateVector(e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->rotateVector(e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // more vectors on the axis
+  fixture = la1->rotateVector(2*e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
+  fixture = la2->rotateVector(2*e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
+  fixture = la3->rotateVector(2*e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
+
+  fixture = la1->rotateVector(2*e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
+  fixture = la2->rotateVector(2*e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
+  fixture = la3->rotateVector(2*e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
+
+  // negative factors
+  fixture = la1->rotateVector(-1*e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(-1*e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(-1*e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+
+  fixture = la1->rotateVector(-1*e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(-1*e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(-1*e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+
+
+
+  // now the real rotations
+  // e2 around e1
+  fixture = la1->rotateVector(e2,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la1->rotateVector(e2,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  fixture = la1->rotateVector(e2,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la1->rotateVector(e2,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+
+  // e3 around e2
+  fixture = la2->rotateVector(e3,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = la2->rotateVector(e3,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(e3,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  fixture = la2->rotateVector(e3,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // e1 around e3
+  fixture = la3->rotateVector(e1,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la3->rotateVector(e1,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(e1,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la3->rotateVector(e1,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+
+
+  // and some rotation around the plane lines
+
+  // Vectors on the line
+  fixture = lp1->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = lp1->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+
+  fixture = lp2->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = lp2->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = lp3->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = lp3->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+
+  // the real stuff
+  fixture = lp1->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,1,0));
+  fixture = lp2->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(0,1,1));
+  fixture = lp3->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,0,1));
+
+  fixture = lp1->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = lp2->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = lp3->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+}

src/unittests/LineUnittest.hpp

@@ -21 +21 @@
   CPPUNIT_TEST ( isContainedTest );
   CPPUNIT_TEST ( intersectionTest );
+  CPPUNIT_TEST ( rotationTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -31 +32 @@
   void isContainedTest();
   void intersectionTest();
+  void rotationTest();
 
 private:

src/unittests/vectorunittest.cpp

@@ -230 +230 @@
 };
 
-/** UnitTest for vector rotations.
- */
-void VectorTest::VectorRotationTest()
-{
-  fixture = Vector(-1.,0.,0.);
-
-  // zero vector does not change
-  fixture = RotateVector(zero,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  fixture = RotateVector(zero, two, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
-
-  // vector on axis does not change
-  fixture = RotateVector(unit,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // rotations
-  fixture = RotateVector(otherunit, unit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( Vector(0.,-1.,0.), fixture );
-
-  fixture = RotateVector(otherunit, unit, 2. * M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(otherunit,unit, 0);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(Vector(0.,0.,1.), notunit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-}
-
 /**
  * UnitTest for Vector::IsInParallelepiped().

src/unittests/vectorunittest.hpp

@@ -28 +28 @@
     CPPUNIT_TEST ( NormalsTest );
     CPPUNIT_TEST ( LineIntersectionTest );
-    CPPUNIT_TEST ( VectorRotationTest );
     CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();

src/vector.cpp

@@ -213 +213 @@
 {
   Vector tmp;
-  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
-  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
-  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
+  tmp[0] = x[1]* y[2] - x[2]* y[1];
+  tmp[1] = x[2]* y[0] - x[0]* y[2];
+  tmp[2] = x[0]* y[1] - x[1]* y[0];
   (*this) = tmp;
 };
@@ -232 +232 @@
   *this -= tmp;
 };
-
-/** Calculates the minimum distance vector of this vector to the plane.
- * \param *out output stream for debugging
- * \param *PlaneNormal normal of plane
- * \param *PlaneOffset offset of plane
- * \return distance to plane
- * \return distance vector onto to plane
- */
-Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
-{
-  Vector temp = (*this) - PlaneOffset;
-  temp.MakeNormalTo(PlaneNormal);
-  temp.Scale(-1.);
-  // then add connecting vector from plane to point
-  temp += (*this)-PlaneOffset;
-  double sign = temp.ScalarProduct(PlaneNormal);
-  if (fabs(sign) > MYEPSILON)
-    sign /= fabs(sign);
-  else
-    sign = 0.;
-
-  temp.Normalize();
-  temp.Scale(sign);
-  return temp;
-};
-
 
 /** Calculates the minimum distance of this vector to the plane.

src/vector.hpp

@@ -44 +44 @@
 
   double DistanceSquared(const Vector &y) const;
-  Vector GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
   double DistanceToSpace(const Space& space) const;
   double PeriodicDistance(const Vector &y, const double * const cell_size) const;

src/vector_ops.cpp

@@ -111 +111 @@
   return true;
 };
-
-/** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
- * \param *axis rotation axis
- * \param alpha rotation angle in radian
- */
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha)
-{
-  Vector a,y;
-  Vector res;
-  // normalise this vector with respect to axis
-  a = vec;
-  a.ProjectOntoPlane(axis);
-  // construct normal vector
-  try {
-    y = Plane(axis,a,0).getNormal();
-  }
-  catch (MathException &excp) {
-    // The normal vector cannot be created if there is linar dependency.
-    // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
-    return vec;
-  }
-  y.Scale(vec.Norm());
-  // scale normal vector by sine and this vector by cosine
-  y.Scale(sin(alpha));
-  a.Scale(cos(alpha));
-  res = vec.Projection(axis);
-  // add scaled normal vector onto this vector
-  res += y;
-  // add part in axis direction
-  res += a;
-  return res;
-};

src/vector_ops.hpp

@@ -10 +10 @@
 
 bool LSQdistance(Vector &res,const Vector **vectors, int num);
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha);
 
 #endif /* VECTOR_OPS_HPP_ */