Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision fbcbf50f90227240f64499afee215cd181250084)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision a693d612bf86be19ef9155b783060502007dcb88)
@@ -46,7 +46,10 @@
 
 #include <boost/assign.hpp>
+#include <boost/math/quaternion.hpp>
 
 #include "CodePatterns/Assert.hpp"
 #include "CodePatterns/Log.hpp"
+
+#include "LinearAlgebra/Line.hpp"
 
 #include "Fragmentation/Exporters/SphericalPointDistribution.hpp"
@@ -77,4 +80,73 @@
 }
 
+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
  */
@@ -126,4 +198,20 @@
     CPPUNIT_ASSERT_EQUAL( expected, remaining );
   }
+
+  // test with one point, full rotation
+  {
+    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
+    SphericalPointDistribution::Polygon_t polygon;
+    polygon += RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI);
+    SphericalPointDistribution::Polygon_t newpolygon =
+        SPD.get<2>();
+    SphericalPointDistribution::Polygon_t expected;
+    expected += RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
+    SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
 }
 
@@ -149,5 +237,5 @@
   }
 
-  // test with one point, just a flip of axis
+  // test with one point, just a flip of x and y axis
   {
     SphericalPointDistribution::Polygon_t polygon;
@@ -157,4 +245,9 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
     SphericalPointDistribution::Polygon_t remaining =
         SphericalPointDistribution::matchSphericalPointDistributions(
@@ -194,4 +287,9 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
     SphericalPointDistribution::Polygon_t remaining =
         SphericalPointDistribution::matchSphericalPointDistributions(
@@ -231,4 +329,9 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
     SphericalPointDistribution::Polygon_t remaining =
         SphericalPointDistribution::matchSphericalPointDistributions(
@@ -268,4 +371,9 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
     SphericalPointDistribution::Polygon_t remaining =
         SphericalPointDistribution::matchSphericalPointDistributions(
@@ -305,4 +413,9 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
     SphericalPointDistribution::Polygon_t remaining =
         SphericalPointDistribution::matchSphericalPointDistributions(
@@ -342,9 +455,14 @@
     SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-  }
-}
+    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
+        iter != expected.end(); ++iter) {
+      std::swap((*iter)[0], (*iter)[1]);
+      (*iter)[0] *= -1.;
+    }
+   SphericalPointDistribution::Polygon_t remaining =
+        SphericalPointDistribution::matchSphericalPointDistributions(
+            polygon,
+            newpolygon);
+    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+  }
+}
Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision fbcbf50f90227240f64499afee215cd181250084)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision a693d612bf86be19ef9155b783060502007dcb88)
@@ -22,4 +22,5 @@
 {
     CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
+    CPPUNIT_TEST ( QuaternionTest );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
     CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
@@ -34,4 +35,5 @@
       void setUp();
       void tearDown();
+      void QuaternionTest();
       void matchSphericalPointDistributionsTest_2();
       void matchSphericalPointDistributionsTest_3();
