Index: src/Fragmentation/Exporters/SphericalPointDistribution.cpp
===================================================================
--- src/Fragmentation/Exporters/SphericalPointDistribution.cpp	(revision 36bd5961a16a97a5b7a28e6cac1b968b75d5347e)
+++ src/Fragmentation/Exporters/SphericalPointDistribution.cpp	(revision fe90ab76a0488cae2c9dc4e5dbf47938334d5839)
@@ -480,4 +480,5 @@
       for (IndexList_t::iterator iter = _indices.begin();
           (iter != _indices.end()) && (!_MCS.foundflag);) {
+
         // check whether we can stay in the current bin or have to move on to next one
         if (*_remainiter == 0) {
@@ -490,5 +491,6 @@
           }
         }
-        // advance in matching to same position
+
+        // advance in matching to current bin to fill in
         const size_t OldIndex = std::distance(_remainingweights.begin(), _remainiter);
         while (_matching.size() <= OldIndex) { // add empty lists of new bin is opened
@@ -498,4 +500,19 @@
         IndexTupleList_t::iterator filliniter = _matching.begin();
         std::advance(filliniter, OldIndex);
+
+        // check whether connection between bins' indices and candidate is satisfied
+        {
+          adjacency_t::const_iterator finder = _MCS.adjacency.find(*iter);
+          ASSERT( finder != _MCS.adjacency.end(),
+              "recurseMatchings() - "+toString(*iter)+" is not in adjacency list.");
+          if ((!filliniter->empty())
+              && (finder->second.find(*filliniter->begin()) == finder->second.end())) {
+            LOG(5, "DEBUG; Skipping index " << *iter
+                << " as is not connected to current set." << *filliniter << ".");
+            ++iter; // note that for loop does not contain incrementor
+            continue;
+          }
+        }
+
         // add index to matching
         filliniter->push_back(*iter);
@@ -532,4 +549,18 @@
 }
 
+SphericalPointDistribution::MatchingControlStructure::MatchingControlStructure(
+    const adjacency_t &_adjacency,
+    const VectorArray_t &_oldpoints,
+    const VectorArray_t &_newpoints,
+    const WeightsArray_t &_weights
+    ) :
+  foundflag(false),
+  bestL2(std::numeric_limits<double>::max()),
+  adjacency(_adjacency),
+  oldpoints(_oldpoints),
+  newpoints(_newpoints),
+  weights(_weights)
+{}
+
 /** Finds combinatorially the best matching between points in \a _polygon
  * and \a _newpolygon.
@@ -551,23 +582,23 @@
  */
 SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
-    const WeightedPolygon_t &_polygon,
-    Polygon_t &_newpolygon
+    const WeightedPolygon_t &_polygon
     )
 {
-  MatchingControlStructure MCS;
-  MCS.foundflag = false;
-  MCS.bestL2 = std::numeric_limits<double>::max();
   // transform lists into arrays
+  VectorArray_t oldpoints;
+  VectorArray_t newpoints;
+  WeightsArray_t weights;
   for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
       iter != _polygon.end(); ++iter) {
-    MCS.oldpoints.push_back(iter->first);
-    MCS.weights.push_back(iter->second);
-  }
-  MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
+    oldpoints.push_back(iter->first);
+    weights.push_back(iter->second);
+  }
+  newpoints.insert(newpoints.begin(), points.begin(), points.end() );
+  MatchingControlStructure MCS(adjacency, oldpoints, newpoints, weights);
 
   // search for bestmatching combinatorially
   {
     // translate polygon into vector to enable index addressing
-    IndexList_t indices(_newpolygon.size());
+    IndexList_t indices(points.size());
     std::generate(indices.begin(), indices.end(), UniqueNumber);
     IndexTupleList_t matching;
@@ -594,5 +625,5 @@
   // combine multiple points and create simple IndexList from IndexTupleList
   const SphericalPointDistribution::IndexList_t IndexList =
-      joinPoints(_newpolygon, MCS.newpoints, MCS.bestmatching);
+      joinPoints(points, MCS.newpoints, MCS.bestmatching);
 
   return IndexList;
@@ -781,24 +812,93 @@
 }
 
+void SphericalPointDistribution::initSelf(const int _NumberOfPoints)
+{
+  switch (_NumberOfPoints)
+  {
+    case 0:
+      points = get<0>();
+      adjacency = getConnections<0>();
+      break;
+    case 1:
+      points = get<1>();
+      adjacency = getConnections<1>();
+      break;
+    case 2:
+      points = get<2>();
+      adjacency = getConnections<2>();
+      break;
+    case 3:
+      points = get<3>();
+      adjacency = getConnections<3>();
+      break;
+    case 4:
+      points = get<4>();
+      adjacency = getConnections<4>();
+      break;
+    case 5:
+      points = get<5>();
+      adjacency = getConnections<5>();
+      break;
+    case 6:
+      points = get<6>();
+      adjacency = getConnections<6>();
+      break;
+    case 7:
+      points = get<7>();
+      adjacency = getConnections<7>();
+      break;
+    case 8:
+      points = get<8>();
+      adjacency = getConnections<8>();
+      break;
+    case 9:
+      points = get<9>();
+      adjacency = getConnections<9>();
+      break;
+    case 10:
+      points = get<10>();
+      adjacency = getConnections<10>();
+      break;
+    case 11:
+      points = get<11>();
+      adjacency = getConnections<11>();
+      break;
+    case 12:
+      points = get<12>();
+      adjacency = getConnections<12>();
+      break;
+    case 14:
+      points = get<14>();
+      adjacency = getConnections<14>();
+      break;
+    default:
+      ASSERT(0, "SphericalPointDistribution::initSelf() - cannot deal with the case "
+          +toString(_NumberOfPoints)+".");
+  }
+  LOG(3, "DEBUG: Ideal polygon is " << points);
+}
 
 SphericalPointDistribution::Polygon_t
-SphericalPointDistribution::matchSphericalPointDistributions(
-    const SphericalPointDistribution::WeightedPolygon_t &_polygon,
-    SphericalPointDistribution::Polygon_t &_newpolygon
-    )
+SphericalPointDistribution::getRemainingPoints(
+    const WeightedPolygon_t &_polygon,
+    const int _N)
 {
   SphericalPointDistribution::Polygon_t remainingpoints;
 
+  // initialze to given number of points
+  initSelf(_N);
   LOG(2, "INFO: Matching old polygon " << _polygon
-      << " with new polygon " << _newpolygon);
-
-  if (_polygon.size() == _newpolygon.size()) {
-    // same number of points desired as are present? Do nothing
-    LOG(2, "INFO: There are no vacant points to return.");
+      << " with new polygon " << points);
+
+  // check whether any points will remain vacant
+  int RemainingPoints = _N;
+  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
+      iter != _polygon.end(); ++iter)
+    RemainingPoints -= iter->second;
+  if (RemainingPoints == 0)
     return remainingpoints;
-  }
-
-  if (_polygon.size() > 0) {
-    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
+
+  if (_N > 0) {
+    IndexList_t bestmatching = findBestMatching(_polygon);
     LOG(2, "INFO: Best matching is " << bestmatching);
 
@@ -819,5 +919,5 @@
     newSet.indices.resize(NumberIds, -1);
     std::copy(beginiter, enditer, newSet.indices.begin());
-    std::copy(_newpolygon.begin(),_newpolygon.end(), std::back_inserter(newSet.polygon));
+    std::copy(points.begin(),points.end(), std::back_inserter(newSet.polygon));
 
     // determine rotation angles to align the two point distributions with
@@ -927,4 +1027,4 @@
     return remainingpoints;
   } else
-    return _newpolygon;
-}
+    return points;
+}
Index: src/Fragmentation/Exporters/SphericalPointDistribution.hpp
===================================================================
--- src/Fragmentation/Exporters/SphericalPointDistribution.hpp	(revision 36bd5961a16a97a5b7a28e6cac1b968b75d5347e)
+++ src/Fragmentation/Exporters/SphericalPointDistribution.hpp	(revision fe90ab76a0488cae2c9dc4e5dbf47938334d5839)
@@ -74,19 +74,14 @@
   }
 
-  /** Matches a given spherical distribution with another containing more
-   * points.
+  /** Returns vacant spots to fill to get a complete spherical point distribution from
+   * given points \a _polygon, containing then \a _N in total.
    *
    * This is a helper to determine points where to best insert saturation
    * hydrogens.
    *
-   * \param _polygon current occupied positions
-   * \param _newpolygon ideal distribution to match best with current occupied
-   *        positions
-   * \return remaining vacant positions relative to \a _polygon
-   */
-  static Polygon_t matchSphericalPointDistributions(
-      const WeightedPolygon_t &_polygon,
-      Polygon_t &_newpolygon
-      );
+   * \param _polygon already filled places to match
+   * \param _N desired total number fo points
+   */
+  Polygon_t getRemainingPoints(const WeightedPolygon_t &_polygon, const int _N);
 
   //!> default radius of the spherical distribution
@@ -131,4 +126,17 @@
 
 private:
+  //!> points for the ideal distribution
+  Polygon_t points;
+  //!> connection information between these ideal points
+  adjacency_t adjacency;
+
+  /** Initialize inner status (points and adjacency) to desired number of
+   * points.
+   *
+   * \param _N number of points
+   */
+  void initSelf(const int _N);
+
+private:
   //!> grant unit tests access to private parts
   friend class SphericalPointDistributionTest;
@@ -143,10 +151,17 @@
 
   struct MatchingControlStructure {
+    MatchingControlStructure(
+        const adjacency_t &_adjacency,
+        const VectorArray_t &_oldpoints,
+        const VectorArray_t &_newpoints,
+        const WeightsArray_t &_weights
+        );
     bool foundflag;
     double bestL2;
+    const adjacency_t &adjacency;
+    const VectorArray_t oldpoints;
+    const VectorArray_t newpoints;
+    const WeightsArray_t weights;
     IndexTupleList_t bestmatching;
-    VectorArray_t oldpoints;
-    VectorArray_t newpoints;
-    WeightsArray_t weights;
   };
 
@@ -160,8 +175,5 @@
       );
 
-  static IndexList_t findBestMatching(
-      const WeightedPolygon_t &_polygon,
-      Polygon_t &_newpolygon
-      );
+  IndexList_t findBestMatching(const WeightedPolygon_t &_polygon);
 
   static IndexList_t joinPoints(
Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision 36bd5961a16a97a5b7a28e6cac1b968b75d5347e)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp	(revision fe90ab76a0488cae2c9dc4e5dbf47938334d5839)
@@ -736,7 +736,7 @@
 }
 
-/** UnitTest for matchSphericalPointDistributions() with two points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_2()
+/** UnitTest for getRemainingPoints() with two points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_2()
 {
   SphericalPointDistribution SPD(1.);
@@ -745,12 +745,8 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair(Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
-        SPD.get<2>();
     SphericalPointDistribution::Polygon_t expected;
     expected += Vector(-1.,0.,0.);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 2);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -761,12 +757,8 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
-        SPD.get<2>();
     SphericalPointDistribution::Polygon_t expected;
     expected += Vector(0.,-1.,0.);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 2);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -777,12 +769,8 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,0.,-1.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
-        SPD.get<2>();
     SphericalPointDistribution::Polygon_t expected;
     expected += Vector(0.,0.,1.);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 2);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -794,20 +782,16 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
-    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 );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with three points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_3()
+        SPD.getRemainingPoints(polygon, 2);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with three points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_3()
 {
   SphericalPointDistribution SPD(1.);
@@ -817,12 +801,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<3>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 3);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -833,7 +814,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<3>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -843,7 +823,5 @@
     }
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 3);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -855,13 +833,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-0.5, sqrt(3)*0.5,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<3>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 3);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -878,7 +853,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-0.5, sqrt(3)*0.5,0.), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<3>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -887,7 +861,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 3);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -908,7 +880,5 @@
     SphericalPointDistribution::Polygon_t expected; // empty cause none are vacant
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 3);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -931,19 +901,17 @@
     SphericalPointDistribution::Polygon_t expected; // empty cause none are vacant
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with four points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_4()
+        SPD.getRemainingPoints(polygon, 3);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with four points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_4()
 {
   SphericalPointDistribution SPD(1.);
@@ -953,12 +921,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 4);
     //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
         CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -969,7 +934,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -979,7 +943,5 @@
     }
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 4);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -991,13 +953,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 4);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1013,13 +972,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 4);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1036,7 +992,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1045,7 +1000,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 4);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1062,14 +1015,11 @@
     polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 1);
     polygon += std::make_pair( Vector(-1./3.0, -M_SQRT2/3.0, M_SQRT2/sqrt(3)), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
     expected.pop_front(); // remove third point
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 4);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1087,7 +1037,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, -M_SQRT2/3.0, M_SQRT2/sqrt(3)), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<4>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1097,20 +1046,18 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with four points and weights
+        SPD.getRemainingPoints(polygon, 4);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with four points and weights
  * not all equal to one.
  */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_multiple()
+void SphericalPointDistributionTest::getRemainingPointsTest_multiple()
 {
   SphericalPointDistribution SPD(1.);
@@ -1126,7 +1073,5 @@
     expected += Vector(-0.5773502691896,-5.551115123126e-17,-0.8164965809277);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 4);
 //    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
@@ -1145,8 +1090,6 @@
     expected += Vector(-0.3535533905933,-0.3535533905933,-0.8660254037844);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
+        SPD.getRemainingPoints(polygon, 5);
+    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1165,7 +1108,5 @@
     expected += Vector(0.3535534025157,-0.3535533856548,0.8660254009332);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 5);
 //    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
@@ -1181,10 +1122,8 @@
         SPD.get<6>();
     SphericalPointDistribution::Polygon_t expected;
+    expected += Vector(0.,0.,1.);
     expected += Vector(0.,0.,-1.);
-    expected += Vector(0.,0.,1.);
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 6);
 //    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
@@ -1193,7 +1132,7 @@
 }
 
-/** UnitTest for matchSphericalPointDistributions() with five points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_5()
+/** UnitTest for getRemainingPoints() with five points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_5()
 {
   SphericalPointDistribution SPD(1.);
@@ -1203,12 +1142,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 5);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1219,7 +1155,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -1229,7 +1164,5 @@
     }
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 5);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1241,13 +1174,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 5);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1264,7 +1194,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180.*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180.*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1273,7 +1202,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180.*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 5);
     // the three remaining points sit on a plane that may be rotated arbitrarily
     // so we cannot simply check for equality between expected and remaining
@@ -1295,14 +1222,11 @@
     polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
     polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
     expected.pop_front(); // remove third point
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 5);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1320,7 +1244,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<5>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1330,19 +1253,17 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with six points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_6()
+        SPD.getRemainingPoints(polygon, 5);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with six points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_6()
 {
   SphericalPointDistribution SPD(1.);
@@ -1352,12 +1273,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 6);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1368,7 +1286,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -1378,7 +1295,5 @@
     }
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 6);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1390,13 +1305,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second spoint
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 6);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1413,7 +1325,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second spoint
@@ -1422,7 +1333,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 6);
     // the four remaining points sit on a plane that may have been rotated arbitrarily
     // so we cannot simply check for equality between expected and remaining
@@ -1444,14 +1353,11 @@
     polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
     polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
     expected.pop_front(); // remove third point
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 6);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1469,7 +1375,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<6>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1479,19 +1384,17 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with seven points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_7()
+        SPD.getRemainingPoints(polygon, 6);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with seven points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_7()
 {
   SphericalPointDistribution SPD(1.);
@@ -1501,12 +1404,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 7);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1517,7 +1417,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -1527,7 +1426,5 @@
     }
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 7);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1539,13 +1436,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 7);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1562,7 +1456,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1571,7 +1464,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 7);
     // the five remaining points sit on a plane that may have been rotated arbitrarily
     // so we cannot simply check for equality between expected and remaining
@@ -1593,14 +1484,11 @@
     polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
     polygon += std::make_pair( Vector(0.0, 1., 0.0), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
     expected.pop_front(); // remove third point
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 7);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1618,7 +1506,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(0.0, 1., 0.0), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<7>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1628,19 +1515,17 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
-
-/** UnitTest for matchSphericalPointDistributions() with eight points
- */
-void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_8()
+        SPD.getRemainingPoints(polygon, 7);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
+
+/** UnitTest for getRemainingPoints() with eight points
+ */
+void SphericalPointDistributionTest::getRemainingPointsTest_8()
 {
   SphericalPointDistribution SPD(1.);
@@ -1650,12 +1535,9 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 8);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1666,7 +1548,6 @@
     SphericalPointDistribution::WeightedPolygon_t polygon;
     polygon += std::make_pair( Vector(0.,1.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
@@ -1675,8 +1556,6 @@
       (*iter)[0] *= -1.;
     }
-   SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 8);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1688,13 +1567,10 @@
     polygon += std::make_pair( Vector(1.,0.,0.), 1);
     polygon += std::make_pair( Vector(-1.,0.,0.), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
-    expected.pop_front(); // remove first point
-    expected.pop_front(); // remove second point
-    SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+    expected.pop_front(); // remove first point
+    expected.pop_front(); // remove second point
+    SphericalPointDistribution::Polygon_t remaining =
+        SPD.getRemainingPoints(polygon, 8);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1711,7 +1587,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1720,7 +1595,5 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 8);
     // the six remaining points sit on two planes that may have been rotated arbitrarily
     // so we cannot simply check for equality between expected and remaining
@@ -1746,14 +1619,11 @@
     polygon += std::make_pair( Vector(-1., 0.0, 0.0), 1);
     polygon += std::make_pair( Vector(-1./3.0, 2.0*M_SQRT2/3.0, 0.0), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
     expected.pop_front(); // remove third point
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
+        SPD.getRemainingPoints(polygon, 8);
 //    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
@@ -1771,7 +1641,6 @@
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1., 0.0, 0.0), 47.6/180*M_PI), 1);
     polygon += std::make_pair(RotationAxis.rotateVector(Vector(-1./3.0, 2.0*M_SQRT2/3.0, 0.0), 47.6/180*M_PI), 1);
-    SphericalPointDistribution::Polygon_t newpolygon =
+    SphericalPointDistribution::Polygon_t expected =
         SPD.get<8>();
-    SphericalPointDistribution::Polygon_t expected = newpolygon;
     expected.pop_front(); // remove first point
     expected.pop_front(); // remove second point
@@ -1781,13 +1650,11 @@
       *iter = RotationAxis.rotateVector(*iter, 47.6/180*M_PI);
     SphericalPointDistribution::Polygon_t remaining =
-        SphericalPointDistribution::matchSphericalPointDistributions(
-            polygon,
-            newpolygon);
-//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
-    // also slightly perturbed
-    const double amplitude = 0.05;
-    perturbPolygon(polygon, amplitude);
-    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
-  }
-}
+        SPD.getRemainingPoints(polygon, 8);
+//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
+    // also slightly perturbed
+    const double amplitude = 0.05;
+    perturbPolygon(polygon, amplitude);
+    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, amplitude) );
+  }
+}
Index: src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp
===================================================================
--- src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision 36bd5961a16a97a5b7a28e6cac1b968b75d5347e)
+++ src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp	(revision fe90ab76a0488cae2c9dc4e5dbf47938334d5839)
@@ -65,12 +65,12 @@
     CPPUNIT_TEST ( getConnectionTest<12> );
     CPPUNIT_TEST ( getConnectionTest<14> );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_4 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_5 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_6 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_7 );
-    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_8 );
-//    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_multiple );
+    CPPUNIT_TEST ( getRemainingPointsTest_2 );
+    CPPUNIT_TEST ( getRemainingPointsTest_3 );
+    CPPUNIT_TEST ( getRemainingPointsTest_4 );
+    CPPUNIT_TEST ( getRemainingPointsTest_5 );
+    CPPUNIT_TEST ( getRemainingPointsTest_6 );
+    CPPUNIT_TEST ( getRemainingPointsTest_7 );
+    CPPUNIT_TEST ( getRemainingPointsTest_8 );
+    CPPUNIT_TEST ( getRemainingPointsTest_multiple );
     CPPUNIT_TEST_SUITE_END();
 
@@ -81,12 +81,12 @@
       void areEqualToWithinBoundsTest();
       void joinPointsTest();
-      void matchSphericalPointDistributionsTest_2();
-      void matchSphericalPointDistributionsTest_3();
-      void matchSphericalPointDistributionsTest_4();
-      void matchSphericalPointDistributionsTest_5();
-      void matchSphericalPointDistributionsTest_6();
-      void matchSphericalPointDistributionsTest_7();
-      void matchSphericalPointDistributionsTest_8();
-      void matchSphericalPointDistributionsTest_multiple();
+      void getRemainingPointsTest_2();
+      void getRemainingPointsTest_3();
+      void getRemainingPointsTest_4();
+      void getRemainingPointsTest_5();
+      void getRemainingPointsTest_6();
+      void getRemainingPointsTest_7();
+      void getRemainingPointsTest_8();
+      void getRemainingPointsTest_multiple();
 
 private:
