source: src/samples/bspline.cpp@ ac6d04

/**
 * @file   bspline.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Nov 21 13:27:22 2011
 *
 * @brief  B-Splines for molecular dynamics.
 *
 */

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

#ifndef BSPLINE_DEGREE
#error BSPLINE_DEGREE not defined.
#endif

#if BSPLINE_DEGREE < 3 || BSPLINE_DEGREE > 6
#error Please choose a B-Spline degree between 3 and 6
#endif

#include "base/helper.hpp"
#include "base/math.hpp"
#include "samples/bspline.hpp"

#define POW(x,y) Helper::pow(x,y)

using namespace VMG;

BSpline::BSpline(const vmg_float& width) :
  spline_nom(BSPLINE_DEGREE/2+1),
  spline_denom(BSPLINE_DEGREE/2+1),
  potential_nom(BSPLINE_DEGREE/2+2),
  potential_denom(BSPLINE_DEGREE/2+2),
  intervals(BSPLINE_DEGREE/2+1),
  R(width)
{
  for (unsigned int i=0; i<intervals.size(); ++i)
    intervals[i] = R * ( -1.0 + 2.0 / static_cast<vmg_float>(BSPLINE_DEGREE) * (i + BSPLINE_DEGREE/2 + 1));

#if BSPLINE_DEGREE == 3

  spline_nom[0] = Polynomial(2, 81.0*POW(R,2), 0.0, -243.0);
  spline_nom[1] = Polynomial(2, 243.0*POW(R,2), -486.0*R, 243.0);

  spline_denom[0] = Polynomial(0, 16.0 * Math::pi * POW(R,5));
  spline_denom[1] = Polynomial(0, 32.0 * Math::pi * POW(R,5));

  potential_nom[0] = Polynomial(5,
                                0.0,
                                -195.0*POW(R,4),
                                0.0,
                                270.0*POW(R,2),
                                0.0,
                                -243.0);

  potential_nom[1] = Polynomial(5,
                                2.0*POW(R,5),
                                -405.0*POW(R,4),
                                0.0,
                                810.0*POW(R,2),
                                -810.0*R,
                                243.0);

  potential_nom[2] = Polynomial(0, -1.0);

  potential_denom[0] = Polynomial(0, 80.0*POW(R,5));
  potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
  potential_denom[2] = Polynomial(0, 1.0);

  antid = 13.0 / 48.0 * R;

#else
#error B-Spline degree not supported. Choose 3.
#endif /* BSPLINE_DEGREE */
}