Index: src/Fragmentation/Exporters/SphericalPointDistribution.cpp
===================================================================
--- src/Fragmentation/Exporters/SphericalPointDistribution.cpp	(revision 8e4471912c0fd0ab5f57a3c97bdd1a097df54cb4)
+++ src/Fragmentation/Exporters/SphericalPointDistribution.cpp	(revision b3c0524231a367779d116b7b16afebad1d225967)
@@ -139,5 +139,5 @@
 
 SphericalPointDistribution::Polygon_t removeMatchingPoints(
-    const SphericalPointDistribution::Polygon_t &_returnpolygon,
+    const VectorArray_t &_points,
     const IndexList_t &_matchingindices
     )
@@ -147,15 +147,16 @@
   std::sort(indices.begin(), indices.end());
   LOG(4, "DEBUG: sorted matching is " << indices);
-  IndexArray_t::const_iterator valueiter = indices.begin();
-  SphericalPointDistribution::Polygon_t::const_iterator pointiter =
-      _returnpolygon.begin();
-  for (unsigned int i=0; i< _returnpolygon.size(); ++i, ++pointiter) {
-    // skip all those in values
-    if (*valueiter == i)
-      ++valueiter;
-    else
-      remainingreturnpolygon.push_back(*pointiter);
+  IndexArray_t remainingindices(_points.size(), -1);
+  std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
+  IndexArray_t::iterator remainiter = std::set_difference(
+      remainingindices.begin(), remainingindices.end(),
+      indices.begin(), indices.end(),
+      remainingindices.begin());
+  remainingindices.erase(remainiter, remainingindices.end());
+  LOG(4, "DEBUG: remaining indices are " << remainingindices);
+  for (IndexArray_t::const_iterator iter = remainingindices.begin();
+      iter != remainingindices.end(); ++iter) {
+    remainingreturnpolygon.push_back(_points[*iter]);
   }
-  LOG(4, "DEBUG: remaining indices are " << remainingreturnpolygon);
 
   return remainingreturnpolygon;
@@ -292,5 +293,5 @@
     // determine rotation angles to align the two point distributions with
     // respect to bestmatching
-    SphericalPointDistribution::Polygon_t rotated_newpolygon;
+    VectorArray_t rotated_newpolygon = remainingnew;
     Vector oldCenter;
     {
@@ -309,50 +310,69 @@
       if ((oldCenter - newCenter).NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
         // setup quaternion
-        Vector RotationAxis = newCenter;
-        RotationAxis.VectorProduct(oldCenter);
+        Vector RotationAxis = oldCenter;
+        RotationAxis.VectorProduct(newCenter);
+        Line Axis(zeroVec, RotationAxis);
         RotationAxis.Normalize();
-        const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
-//            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
-        boost::math::quaternion<double> q
-            (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
-        LOG(5, "DEBUG: RotationAxis is " << RotationAxis
-            << ", RotationAngle is " << RotationAngle);
-        LOG(5, "DEBUG: Quaternion describing rotation is " << q);
+        const double RotationAngle = oldCenter.Angle(newCenter); // /(M_PI/2.);
+        LOG(5, "DEBUG: Rotate coordinate system by " << RotationAngle
+            << " around axis " << RotationAxis);
+
+        // apply rotation angles
+        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
+            iter != rotated_newpolygon.end(); ++iter) {
+          Vector &current = *iter;
+          LOG(5, "DEBUG: Original point is " << current);
+          current =  Axis.rotateVector(current, RotationAngle);
+          LOG(5, "DEBUG: Rotated point is " << current);
+        }
+      }
+    }
+    // rotate triangle/line/point around itself to match orientation
+    if (MCS.bestmatching.size() > 1) {
+      if (oldCenter.NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
+        // construct RotationAxis and two points on its plane, defining the angle
+        const Line RotationAxis(zeroVec, oldCenter);
+        Vector oldPosition(rotated_newpolygon[*MCS.bestmatching.begin()]);
+        oldPosition.ProjectOntoPlane(RotationAxis.getDirection());
+        Vector newPosition(remainingold[*MCS.bestmatching.begin()]);
+        newPosition.ProjectOntoPlane(RotationAxis.getDirection());
+
+        // construct reference vector to determine direction of rotation
+        Vector dreiBein(oldPosition);
+        dreiBein.VectorProduct(oldCenter);
+        dreiBein.Normalize();
+        const double sign =
+            (dreiBein.ScalarProduct(newPosition) < 0.) ? -1. : +1.;
+        LOG(6, "DEBUG: oldCenter on plane is " << oldPosition
+            << ", newCenter in plane is " << newPosition
+            << ", and dreiBein is " << dreiBein);
+        const double RotationAngle = sign * oldPosition.Angle(newPosition);
+        LOG(5, "DEBUG: Rotating around self is " << RotationAngle
+            << " around axis " << RotationAxis);
+
 #ifndef NDEBUG
+        // check: first bestmatching in rotated_newpolygon and remainingnew
+        // should now equal
         {
-        // check that rotation works in DEBUG mode
-        boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
-        boost::math::quaternion<double> q_inverse =
-            boost::math::conj(q)/(boost::math::norm(q));
-        p = q * p * q_inverse;
-        boost::math::quaternion<double> identity(1,0,0,0);
-        ASSERT( boost::math::norm(q*q_inverse - identity) < std::numeric_limits<double>::epsilon()*1e4,
-            "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
-            +toString(q*q_inverse)+" into oldCenter "+toString(identity)+".");
-        boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
-        ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4,
-            "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
-            +toString(p)+" into oldCenter "+toString(comparison)+".");
+          const IndexList_t::const_iterator iter = MCS.bestmatching.begin();
+          Vector rotatednew = RotationAxis.rotateVector(
+              rotated_newpolygon[*iter],
+              RotationAngle);
+          LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
+              << " while old was " << remainingold[*iter]);
+          ASSERT( (rotatednew - remainingold[*iter]).Norm()
+              < std::numeric_limits<double>::epsilon()*1e4,
+              "matchSphericalPointDistributions() - orientation rotation does not work as expected.");
         }
 #endif
 
-        // rotate spherical distribution
-        rotated_newpolygon = rotatePolygon(_newpolygon, q);
-        LOG(5, "DEBUG: Rotated new polygon is " << rotated_newpolygon);
-      } else {
-        rotated_newpolygon = _newpolygon;
-      }
-    }
-    // rotate triangle/line/point around itself to match orientation
-    const Line RotationAxis(zeroVec, oldCenter);
-    const double RotationAngle =
-        oldCenter.Angle(remainingold[0])
-        - oldCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
-    LOG(5, "DEBUG: Rotate around self is " << RotationAngle
-        << " around axis " << RotationAxis);
-
-    for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_newpolygon.begin();
-        iter != rotated_newpolygon.end(); ++iter) {
-      RotationAxis.rotateVector(*iter, RotationAngle);
+        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
+            iter != rotated_newpolygon.end(); ++iter) {
+          Vector &current = *iter;
+          LOG(6, "DEBUG: Original point is " << current);
+          current = RotationAxis.rotateVector(current, RotationAngle);
+          LOG(6, "DEBUG: Rotated point is " << current);
+        }
+      }
     }
 
Index: src/Fragmentation/Exporters/unittests/Makefile.am
===================================================================
--- src/Fragmentation/Exporters/unittests/Makefile.am	(revision 8e4471912c0fd0ab5f57a3c97bdd1a097df54cb4)
+++ src/Fragmentation/Exporters/unittests/Makefile.am	(revision b3c0524231a367779d116b7b16afebad1d225967)
@@ -21,5 +21,4 @@
 
 TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
-XFAIL_TESTS += SphericalPointDistributionUnitTest
 check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS) 
 noinst_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)
Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision 8e4471912c0fd0ab5f57a3c97bdd1a097df54cb4)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision b3c0524231a367779d116b7b16afebad1d225967)
@@ -82,72 +82,4 @@
 }
 
-void SphericalPointDistributionTest::QuaternionTest()
-{
-  Vector oldCenter(0.,1.,0.);
-  Vector newCenter(1.,0.,0.);
-
-  {
-    // setup quaternion
-    Vector RotationAxis = newCenter;
-    RotationAxis.VectorProduct(oldCenter);
-    RotationAxis.Normalize();
-    const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
-  //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
-    boost::math::quaternion<double> q
-        (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
-    LOG(5, "DEBUG: RotationAxis is " << RotationAxis
-        << ", RotationAngle is " << RotationAngle);
-    LOG(5, "DEBUG: Quaternion describing rotation is " << q);
-    boost::math::quaternion<double> q_inverse =
-        boost::math::conj(q)/(boost::math::norm(q));
-    LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
-    boost::math::quaternion<double> identity(1,0,0,0);
-    const boost::math::quaternion<double> unity = q*q_inverse;
-    LOG(5, "DEBUG: q * q^-1 is " << unity);
-    CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
-
-    // check that rotation works
-    boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
-    LOG(5, "DEBUG: Original newCenter is " << p);
-    p = p * q_inverse;
-    p = q * p;
-    LOG(5, "DEBUG: Rotated newCenter is " << p);
-    boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
-    LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
-    CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
-  }
-
-  // rotating with angle = 0 flips the vector unwantedly
-  {
-    // setup quaternion
-    Vector RotationAxis = newCenter;
-    RotationAxis.VectorProduct(oldCenter);
-    RotationAxis.Normalize();
-    const double RotationAngle = 0.;
-  //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
-    boost::math::quaternion<double> q
-        (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
-    LOG(5, "DEBUG: RotationAxis is " << RotationAxis
-        << ", RotationAngle is " << RotationAngle);
-    LOG(5, "DEBUG: Quaternion describing rotation is " << q);
-    boost::math::quaternion<double> q_inverse =
-        boost::math::conj(q)/(boost::math::norm(q));
-    LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
-    boost::math::quaternion<double> identity(1,0,0,0);
-    const boost::math::quaternion<double> unity = q*q_inverse;
-    LOG(5, "DEBUG: q * q^-1 is " << unity);
-    CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
-
-    // check that rotation works
-    boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
-    boost::math::quaternion<double> comparison(0., -newCenter[0], newCenter[1], newCenter[2]);
-    LOG(5, "DEBUG: Original newCenter is " << p);
-    p = p * q_inverse;
-    p = q * p;
-    LOG(5, "DEBUG: Rotated newCenter is " << p);
-    LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
-    CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
-  }
-}
 
 /** UnitTest for matchSphericalPointDistributions() with two points
Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision 8e4471912c0fd0ab5f57a3c97bdd1a097df54cb4)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision b3c0524231a367779d116b7b16afebad1d225967)
@@ -22,5 +22,4 @@
 {
     CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
-    CPPUNIT_TEST ( QuaternionTest );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
@@ -35,5 +34,4 @@
       void setUp();
       void tearDown();
-      void QuaternionTest();
       void matchSphericalPointDistributionsTest_2();
       void matchSphericalPointDistributionsTest_3();
