/*
 * SphericalPointDistribution.hpp
 *
 *  Created on: May 29, 2014
 *      Author: heber
 */


#ifndef SPHERICALPOINTDISTRIBUTION_HPP_
#define SPHERICALPOINTDISTRIBUTION_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/IteratorAdaptors.hpp"
#include "CodePatterns/toString.hpp"

#include <algorithm>
#include <cmath>
#include <limits>
#include <list>
#include <vector>
#include <map>

#include "LinearAlgebra/Line.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"

/** contains getters for the VSEPR model for specific number of electrons.
 *
 * This struct contains specialized functions returning a list of Vectors
 * (points in space) to match the VSEPR model for the given number of electrons.
 *
 * This is implemented via template specialization of the function get().
 *
 * These specializations are taken from the python script \b CreateVspeShapes.py
 * by Christian Neuen, 07th May 2009.
 */
struct SphericalPointDistribution
{
  /** Cstor for SphericalPointDistribution, allows setting radius of sphere
   *
   * \param _BondLength desired radius of sphere
   */
  SphericalPointDistribution(const double _Bondlength = 1.) :
    Bondlength(_Bondlength),
    SQRT_3(sqrt(3.0))
  {}

  //!> typedef for the list of points
  typedef std::list<Vector> Polygon_t;

  /** General getter function for the distribution of points on the surface.
   *
   * \warn this function needs to be specialized!
   *
   * \return Polygon_t with points on the surface centered at (0,0,0)
   */
  template <int N> Polygon_t get()
  {
    ASSERT(0, "SphericalPointDistribution::get() - not specialized for "+toString(N)+".");
  }

  /** Matches a given spherical distribution with another containing more
   * points.
   *
   * The idea is to produce a matching from all points in \a _polygon to those
   * in \a _newpolygon in such a way that their distance difference is minimal.
   * As we just look at numbers of points determined by valency, i.e.
   * independent of the number of atoms, we simply go through each of the possible
   * mappings. We stop when the L1 error is below a certain \a threshold,
   * otherwise we pick the matching with the lowest L2 error.
   *
   * 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 Polygon_t &_polygon,
      const Polygon_t &_newpolygon
      );


  //!> default radius of the spherical distribution
  const double Bondlength;
  //!> precalculated value for root of 3
  const double SQRT_3;
};

typedef std::list<unsigned int> IndexList_t;
typedef std::vector<unsigned int> IndexArray_t;
typedef std::vector<Vector> VectorArray_t;
typedef std::vector<double> DistanceArray_t;

DistanceArray_t calculatePairwiseDistances(
    const std::vector<Vector> &_points,
    const IndexList_t &_indices
    )
{
  DistanceArray_t result;
  for (IndexList_t::const_iterator firstiter = _indices.begin();
      firstiter != _indices.end(); ++firstiter) {
    for (IndexList_t::const_iterator seconditer = firstiter;
        seconditer != _indices.end(); ++seconditer) {
      if (firstiter == seconditer)
        continue;
      const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
      result.push_back(distance);
    }
  }
  return result;
}

// class generator: taken from www.cplusplus.com example std::generate
struct c_unique {
  int current;
  c_unique() {current=0;}
  int operator()() {return ++current;}
} UniqueNumber;

/** Returns squared L2 error of the given \a _Matching.
 *
 * We compare the pair-wise distances of each associated matching
 * and check whether these distances each match between \a _old and
 * \a _new.
 *
 * \param _old first set of points (fewer or equal to \a _new)
 * \param _new second set of points
 * \param _Matching matching between the two sets
 * \return pair with L1 and squared L2 error
 */
std::pair<double, double> calculateErrorOfMatching(
    const std::vector<Vector> &_old,
    const std::vector<Vector> &_new,
    const IndexList_t &_Matching)
{
  std::pair<double, double> errors( std::make_pair( 0., 0. ) );

  if (_Matching.size() > 1) {
    // convert matching into two vectors to calculate distance among another

    // calculate all pair-wise distances
    IndexList_t keys(_Matching.size());
    std::generate (keys.begin(), keys.end(), UniqueNumber);
    const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
    const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);

    ASSERT( firstdistances.size() == keys.size()*(keys.size()+1)/2,
        "calculateL2ErrorOfMatching() - unexpected pair-wise distance array size.");
    ASSERT( firstdistances.size() == seconddistances.size(),
        "calculateL2ErrorOfMatching() - mismatch in pair-wise distance array sizes.");
    DistanceArray_t::const_iterator firstiter = firstdistances.begin();
    DistanceArray_t::const_iterator seconditer = seconddistances.begin();
    for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
        ++firstiter, ++seconditer) {
      const double gap = *firstiter - *seconditer;
      // L1 error
      if (errors.first < gap)
        errors.first = gap;
      // L2 error
      errors.second += gap*gap;
    }
  }

  return errors;
}

SphericalPointDistribution::Polygon_t removeMatchingPoints(
    const SphericalPointDistribution::Polygon_t &_points,
    const IndexList_t &_matchingindices
    )
{
  SphericalPointDistribution::Polygon_t remainingpoints;
  IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
  std::sort(indices.begin(), indices.end());
  IndexArray_t::const_iterator valueiter = indices.begin();
  SphericalPointDistribution::Polygon_t::const_iterator pointiter =
      _points.begin();
  for (unsigned int i=0; i< _points.size(); ++i, ++pointiter) {
    // skip all those in values
    if (*valueiter == i)
      ++valueiter;
    else
      remainingpoints.push_back(*pointiter);
  }

  return remainingpoints;
}

/** Rotates a given polygon around x, y, and z axis by the given angles.
 *
 * Essentially, we concentrate on the three points of the polygon to rotate
 * to the correct position. First, we rotate its center via \a angles,
 * then we rotate the "triangle" around itself/\a _RotationAxis by
 * \a _RotationAngle.
 *
 * \param _polygon polygon whose points to rotate
 * \param _angles vector with rotation angles for x,y,z axis
 * \param _RotationAxis
 * \param _RotationAngle
 */
SphericalPointDistribution::Polygon_t rotatePolygon(
    const SphericalPointDistribution::Polygon_t &_polygon,
    const std::vector<double> &_angles,
    const Line &_RotationAxis,
    const double _RotationAngle)
{
  SphericalPointDistribution::Polygon_t rotated_polygon = _polygon;
  RealSpaceMatrix rotation;
  ASSERT( _angles.size() == 3,
      "rotatePolygon() - not exactly "+toString(3)+" angles given.");
  rotation.setRotation(_angles[0] * M_PI/180., _angles[1] * M_PI/180., _angles[2] * M_PI/180.);

  // apply rotation angles
  for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_polygon.begin();
      iter != rotated_polygon.end(); ++iter) {
    *iter = rotation * (*iter);
    _RotationAxis.rotateVector(*iter, _RotationAngle);
  }

  return rotated_polygon;
}

struct MatchingControlStructure {
  bool foundflag;
  double bestL2;
  IndexList_t bestmatching;
  VectorArray_t oldpoints;
  VectorArray_t newpoints;
};

/** Recursive function to go through all possible matchings.
 *
 * \param _MCS structure holding global information to the recursion
 * \param _matching current matching being build up
 * \param _indices contains still available indices
 * \param _matchingsize
 */
void recurseMatchings(
    MatchingControlStructure &_MCS,
    IndexList_t &_matching,
    IndexList_t _indices,
    unsigned int _matchingsize)
{
  //!> threshold for L1 error below which matching is immediately acceptable
  const double L1THRESHOLD = 1e-2;
  if (!_MCS.foundflag) {
    if (_matching.size() < _matchingsize) {
      // go through all indices
      for (IndexList_t::iterator iter = _indices.begin();
          iter != _indices.end();) {
        // add index to matching
        _matching.push_back(*iter);
        // remove index but keep iterator to position (is the next to erase element)
        IndexList_t::iterator backupiter = _indices.erase(iter);
        // recurse with decreased _matchingsize
        recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
        // re-add chosen index and reset index to new position
        _indices.insert(backupiter, _matching.back());
        iter = backupiter;
        // remove index from _matching to make space for the next one
        _matching.pop_back();
      }
      // gone through all indices then exit recursion
      _MCS.foundflag = true;
    } else {
      // calculate errors
      std::pair<double, double> errors = calculateErrorOfMatching(
          _MCS.oldpoints, _MCS.newpoints, _matching);
      if (errors.first < L1THRESHOLD) {
        _MCS.bestmatching = _matching;
        _MCS.foundflag = true;
      }
      if (_MCS.bestL2 > errors.second) {
        _MCS.bestmatching = _matching;
        _MCS.bestL2 = errors.second;
      }
    }
  }
}

SphericalPointDistribution::Polygon_t
SphericalPointDistribution::matchSphericalPointDistributions(
    const SphericalPointDistribution::Polygon_t &_polygon,
    const SphericalPointDistribution::Polygon_t &_newpolygon
    )
{
  SphericalPointDistribution::Polygon_t remainingpoints;
  VectorArray_t remainingold(_polygon.begin(), _polygon.end());
  VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());

  if (_polygon.size() > 0) {
    MatchingControlStructure MCS;
    MCS.foundflag = false;
    MCS.bestL2 = std::numeric_limits<double>::max();
    MCS.oldpoints.insert(MCS.oldpoints.begin(), _polygon.begin(),_polygon.end() );
    MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );

    // search for bestmatching combinatorially
    {
      // translate polygon into vector to enable index addressing
      IndexList_t indices(_newpolygon.size());
      std::generate(indices.begin(), indices.end(), UniqueNumber);
      IndexList_t matching;

      // walk through all matchings
      const unsigned int matchingsize = _polygon.size();
      ASSERT( matchingsize <= indices.size(),
          "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
      recurseMatchings(MCS, matching, indices, matchingsize);
    }

    // determine rotation angles to align the two point distributions with
    // respect to bestmatching
    std::vector<double> angles(3);
    Vector newCenter;
    {
      // calculate center of triangle/line/point consisting of first points of matching
      Vector oldCenter;
      IndexList_t::const_iterator iter = MCS.bestmatching.begin();
      unsigned int i = 0;
      for (; (i<3) && (i<MCS.bestmatching.size()); ++i, ++iter) {
        oldCenter += remainingold[i];
        newCenter += remainingnew[*iter];
      }
      oldCenter *= 1./(double)i;
      newCenter *= 1./(double)i;

      Vector direction(0.,0.,0.);
      for(size_t i=0;i<NDIM;++i) {
        // create new rotation axis
        direction[i] = 1.;
        const Line axis (zeroVec, direction);
        // calculate rotation angle for this axis
        const double alpha = direction.Angle(oldCenter) - direction.Angle(newCenter);
        // perform rotation
        axis.rotateVector(newCenter, alpha);
        // store angle
        angles[i] = alpha;
        // reset direction component for next iteration
        direction[i] = 0.;
      }
    }
    const Line RotationAxis(zeroVec, newCenter);
    const double RotationAngle =
        newCenter.Angle(remainingold[0])
        - newCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);

    // rotate _newpolygon
    SphericalPointDistribution::Polygon_t rotated_newpolygon =
        rotatePolygon(_newpolygon, angles, RotationAxis, RotationAngle);

    // remove all points in matching and return remaining ones
    return removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
  } else
    return _newpolygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength,0.,0.));
  ASSERT( polygon.size() == 2,
      "SphericalPointDistribution::get<2>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength*0.5, SQRT_3*0.5,0.));
  polygon.push_back( Vector(-Bondlength*0.5, -SQRT_3*0.5,0.));
  ASSERT( polygon.size() == 3,
      "SphericalPointDistribution::get<3>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength/3.0, Bondlength*2.0*M_SQRT2/3.0,0.));
  polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
  polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
  ASSERT( polygon.size() == 4,
      "SphericalPointDistribution::get<4>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  polygon.push_back( Vector(0.0, Bondlength, 0.0));
  polygon.push_back( Vector(0.0, -Bondlength*0.5, Bondlength*SQRT_3*0.5));
  polygon.push_back( Vector(0.0, -Bondlength*0.5, -Bondlength*SQRT_3*0.5));
  ASSERT( polygon.size() == 5,
      "SphericalPointDistribution::get<5>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  polygon.push_back( Vector(0.0, Bondlength, 0.0));
  polygon.push_back( Vector(0.0, -Bondlength, 0.0));
  polygon.push_back( Vector(0.0, 0.0, Bondlength));
  polygon.push_back( Vector(0.0, 0.0, -Bondlength));
  ASSERT( polygon.size() == 6,
      "SphericalPointDistribution::get<6>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  polygon.push_back( Vector(0.0, Bondlength, 0.0));
  polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*0.4), Bondlength*sin(M_PI*0.4)));
  polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*0.8), Bondlength*sin(M_PI*0.8)));
  polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*1.2), Bondlength*sin(M_PI*1.2)));
  polygon.push_back( Vector(0.0, Bondlength*cos(M_PI*1.6), Bondlength*sin(M_PI*1.6)));
  ASSERT( polygon.size() == 7,
      "SphericalPointDistribution::get<7>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  polygon.push_back( Vector(-Bondlength/3.0, Bondlength*2.0*M_SQRT2/3.0, 0.0));
  polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
  polygon.push_back( Vector(-Bondlength/3.0, -Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
  polygon.push_back( Vector(Bondlength/3.0, -Bondlength*2.0*M_SQRT2/3.0, 0.0));
  polygon.push_back( Vector(Bondlength/3.0, Bondlength*M_SQRT2/3.0, -Bondlength*M_SQRT2/SQRT_3));
  polygon.push_back( Vector(Bondlength/3.0, Bondlength*M_SQRT2/3.0, Bondlength*M_SQRT2/SQRT_3));
  ASSERT( polygon.size() == 8,
      "SphericalPointDistribution::get<8>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), Bondlength*sin(0.4*M_PI), 0.0));
  polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), -Bondlength*sin(0.4*M_PI), 0.0));
  polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), 0.0, Bondlength*sin(0.4*M_PI)));
  polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), 0.0, -Bondlength*sin(0.4*M_PI)));
  polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(0.25*M_PI), Bondlength*sin(0.8*M_PI)*cos(0.25*M_PI)));
  polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(0.75*M_PI), Bondlength*sin(0.8*M_PI)*cos(0.75*M_PI)));
  polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(1.25*M_PI), Bondlength*sin(0.8*M_PI)*cos(1.25*M_PI)));
  polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), Bondlength*sin(0.8*M_PI)*sin(1.75*M_PI), Bondlength*sin(0.8*M_PI)*cos(1.75*M_PI)));
  ASSERT( polygon.size() == 9,
      "SphericalPointDistribution::get<9>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  const double temp_Bondlength = Bondlength*0.5*SQRT_3;
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, 0.0, temp_Bondlength));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, 0.0, -temp_Bondlength));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(0.25*M_PI), temp_Bondlength*cos(0.25*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(0.75*M_PI), temp_Bondlength*cos(0.75*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(1.25*M_PI), temp_Bondlength*cos(1.25*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*sin(1.75*M_PI), temp_Bondlength*cos(1.75*M_PI)));
  ASSERT( polygon.size() == 10,
      "SphericalPointDistribution::get<10>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  {
    const double temp_Bondlength = Bondlength*sin(0.4*M_PI);
    polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength, 0.0));
    polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(0.4*M_PI),  temp_Bondlength*sin(0.4*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(0.8*M_PI),  temp_Bondlength*sin(0.8*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(1.2*M_PI),  temp_Bondlength*sin(1.2*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.4*M_PI), temp_Bondlength*cos(1.6*M_PI),  temp_Bondlength*sin(1.6*M_PI)));
  }
  {
    const double temp_Bondlength = Bondlength*sin(0.8*M_PI);
    polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(0.2*M_PI),  temp_Bondlength*sin(0.2*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(0.6*M_PI),  temp_Bondlength*sin(0.6*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.0*M_PI),  temp_Bondlength*sin(1.0*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.4*M_PI),  temp_Bondlength*sin(1.4*M_PI)));
    polygon.push_back( Vector(Bondlength*cos(0.8*M_PI), temp_Bondlength*cos(1.8*M_PI),  temp_Bondlength*sin(1.8*M_PI)));
  }
  ASSERT( polygon.size() == 11,
      "SphericalPointDistribution::get<11>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  const double temp_Bondlength = Bondlength*0.5*SQRT_3;
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(0.4*M_PI),  temp_Bondlength*sin(0.4*M_PI)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(0.8*M_PI),  temp_Bondlength*sin(0.8*M_PI)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(1.2*M_PI),  temp_Bondlength*sin(1.2*M_PI)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(1.6*M_PI),  temp_Bondlength*sin(1.6*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(0.4*M_PI),  -temp_Bondlength*sin(0.4*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(0.8*M_PI),  -temp_Bondlength*sin(0.8*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(1.2*M_PI),  -temp_Bondlength*sin(1.2*M_PI)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, -temp_Bondlength*cos(1.6*M_PI),  -temp_Bondlength*sin(1.6*M_PI)));
  ASSERT( polygon.size() == 12,
      "SphericalPointDistribution::get<12>() - polygon has wrong size.");
  return polygon;
}

template <>
SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>()
{
  Polygon_t polygon;
  polygon.push_back( Vector(Bondlength,0.,0.));
  polygon.push_back( Vector(-Bondlength, 0.0, 0.0));
  const double temp_Bondlength = Bondlength*0.5*SQRT_3;
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength, 0.0));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/3.0), temp_Bondlength*sin(M_PI/3.0)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(2.0*M_PI/3.0), temp_Bondlength*sin(2.0*M_PI/3.0)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, -temp_Bondlength, 0.0));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(4.0*M_PI/3.0), temp_Bondlength*sin(4.0*M_PI/3.0)));
  polygon.push_back( Vector(temp_Bondlength/SQRT_3, temp_Bondlength*cos(5.0*M_PI/3.0), temp_Bondlength*sin(5.0*M_PI/3.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/6.0),  temp_Bondlength*sin(M_PI/6.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(M_PI/3.0+M_PI/6.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(2.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(2.0*M_PI/3.0 +M_PI/6.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(3.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(3.0*M_PI/3.0 +M_PI/6.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(4.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(4.0*M_PI/3.0 +M_PI/6.0)));
  polygon.push_back( Vector(-temp_Bondlength/SQRT_3, temp_Bondlength*cos(5.0*M_PI/3.0 +M_PI/6.0),  temp_Bondlength*sin(5.0*M_PI/3.0 +M_PI/6.0)));
  ASSERT( polygon.size() == 14,
      "SphericalPointDistribution::get<14>() - polygon has wrong size.");
  return polygon;
}

#endif /* SPHERICALPOINTDISTRIBUTION_HPP_ */