/**
 * @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

#define _USE_MATH_DEFINES
#include <cmath>

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

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

using namespace VMG;

BSpline::BSpline(const vmg_float& width) :
  R(width)
{
  spline_nom.resize(BSPLINE_DEGREE/2+1);
  spline_denom.resize(BSPLINE_DEGREE/2+1);
  potential_nom.resize(BSPLINE_DEGREE/2+2);
  potential_denom.resize(BSPLINE_DEGREE/2+2);
  intervals.resize(BSPLINE_DEGREE/2+1);

  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*M_PI*POW(R,5));
  spline_denom[1] = Polynomial(0, 32.0*M_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, 320.0*M_PI*POW(R,5));
  potential_denom[1] = Polynomial(0, 640.0*M_PI*POW(R,5));
  potential_denom[2] = Polynomial(0, 4.0*M_PI);

  antid = 39.0 / (64.0 * M_PI * R);

#elif BSPLINE_DEGREE == 4

  spline_nom[0] = Polynomial(3,
			     8.0*POW(R,3),
			     0.0,
			     -48.0*R,
			     48.0);

  spline_nom[1] = Polynomial(3,
			     16.0*POW(R,3),
			     -48.0*POW(R,2),
			     48.0*R,
			     -16.0);

  spline_denom[0] = Polynomial(0, M_PI*POW(R,6));
  spline_denom[1] = Polynomial(0, M_PI*POW(R,6));

  potential_nom[0] = Polynomial(6,
				0.0,
				-21.0*POW(R,5),
				0.0,
				40.0*POW(R,3),
				0.0,
				-72.0*R,
				48.0);

  potential_nom[1] = Polynomial(6,
				1.0*POW(R,6),
				-48.0*POW(R,5),
				0.0,
				160.0*POW(R,3),
				-240.0*POW(R,2),
				144.0*R,
				-32.0);

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

  potential_denom[0] = Polynomial(0, 30.0*M_PI*POW(R,6));
  potential_denom[1] = Polynomial(0, 60.0*M_PI*POW(R,6));
  potential_denom[2] = Polynomial(0, 4.0*M_PI);

  antid = 7.0 / (10.0 * M_PI * R);

#elif BSPLINE_DEGREE == 5

  spline_nom[0] = Polynomial(4,
			     2875.0*POW(R, 4),
			     0.0,
			     -18750.0*POW(R,2),
			     0.0,
			     46875.0);

  spline_nom[1] = Polynomial(4,
			     1375.0*POW(R,4),
			     1250.0*POW(R,3),
			     -18750.0*POW(R,2),
			     31250.0*R,
			     -15625.0);

  spline_nom[2] = Polynomial(4,
			     15625.0*POW(R,4),
			     -62500.0*POW(R,3),
			     93750.0*POW(R,2),
			     -62500.0*R,
			     15625.0);

  spline_denom[0] = Polynomial(0, 256.0*M_PI*POW(R,7));
  spline_denom[1] = Polynomial(0, 128.0*M_PI*POW(R,7));
  spline_denom[2] = Polynomial(0, 512.0*M_PI*POW(R,7));

  potential_nom[0] = Polynomial(7,
				0.0,
				-8393.0*POW(R,6),
				0.0,
				20125.0*POW(R,4),
				0.0,
				-39375.0*POW(R,2),
				0.0,
				46875.0);

  potential_nom[1] = Polynomial(7,
				-POW(R,7),
				-20965.0*POW(R,6),
				0.0,
				48125.0*POW(R,4),
				21875.0*POW(R,3),
				-196875.0*POW(R,2),
				218750.0*R,
				78125.0);

  potential_nom[2] = Polynomial(7,
				874.0*POW(R,7),
				-21875.0*POW(R,6),
				0.0,
				109375.0*POW(R,4),
				-218750.0*POW(R,3),
				196875.0*POW(R,2),
				-87500.0*R,
				15625.0);

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

  potential_denom[0] = Polynomial(0, 10752.0*M_PI*POW(R,7));
  potential_denom[1] = Polynomial(0, 26880.0*M_PI*POW(R,7));
  potential_denom[2] = Polynomial(0, 21504.0*M_PI*POW(R,7));
  potential_denom[3] = Polynomial(0, 4.0*M_PI);

  antid = 1199.0 / (1536.0 * M_PI * R);

#elif BSPLINE_DEGREE == 6

  spline_nom[0] = Polynomial(5,
			     297.0*POW(R,5),
			     0.0,
			     -2430.0*POW(R,3),
			     0.0,
			     10935.0*R,
			     -10935.0);

  spline_nom[1] = Polynomial(5,
			     459.0*POW(R,5),
			     2025.0*POW(R,4),
			     -17010.0*POW(R,3),
			     36450.0*POW(R,2),
			     -32805.0*R,
			     10935.0);

  spline_nom[2] = Polynomial(5,
			     2187.0*POW(R,5),
			     -10935.0*POW(R,4),
			     21870.0*POW(R,3),
			     -21870.0*POW(R,2),
			     10935.0*R,
			     -2187.0);

  spline_denom[0] = Polynomial(0, 20.0*M_PI*POW(R,8));
  spline_denom[1] = Polynomial(0, 40.0*M_PI*POW(R,8));
  spline_denom[2] = Polynomial(0, 40.0*M_PI*POW(R,8));

  potential_nom[0] = Polynomial(8,
				0.0,
				-956.0*POW(R,7),
				0.0,
				2772.0*POW(R,5),
				0.0,
				-6804.0*POW(R,3),
				0.0,
				14580.0*R,
				-10935.0);

  potential_nom[1] = Polynomial(8,
				-5.0*POW(R,8),
				-5676.0*POW(R,7),
				0.0,
				12852.0*POW(R,5),
				28350.0*POW(R,4),
				-142884.0*POW(R,3),
				204120.0*POW(R,2),
				-131220.0*R,
				32805.0);

  potential_nom[2] = Polynomial(8,
				169.0*POW(R,8),
				-2916.0*POW(R,7),
				0.0,
				20412.0*POW(R,5),
				-51030.0*POW(R,4),
				61236.0*POW(R,3),
				-40824.0*POW(R,2),
				14580.0*R,
				-2187.0);

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

  potential_denom[0] = Polynomial(0, 1120.0*M_PI*POW(R,8));
  potential_denom[1] = Polynomial(0, 6720.0*M_PI*POW(R,8));
  potential_denom[2] = Polynomial(0, 2240.0*M_PI*POW(R,8));
  potential_denom[3] = Polynomial(0, 4.0*M_PI);

  antid = 239.0 / (280.0 * M_PI * R);

#endif /* BSPLINE_DEGREE */
}

vmg_float BSpline::EvaluateDerivative(const Vector& v, const int& dv)
{
  const vmg_float& v1 = v[dv];
  const vmg_float L = v.Length();
  vmg_float result = 0.0;

#if BSPLINE_DEGREE == 3

  if (L < intervals[0])
    result = v1*(243.0*L*L-135.0*R*R) / (80.0*M_PI*POW(R,5));
  else if (L < intervals[1])
    result = (81.0 * v1 * (POW(R,4)-6.0*POW(L,2)*POW(R,2)+8.0*POW(L,3)*R-8.0*POW(L,4))) / (128.0*M_PI*L*POW(R,5));
  else
    result = -1.0*v1 / (4*M_PI*POW(L,3));

#elif BSPLINE_DEGREE == 4

  if (L < intervals[0])
    result = 8.0*v1 * (-5.0*POW(R,3) + 18.0*R*POW(L,2) - 15.0*POW(L,3)) / (15.0*M_PI*POW(R,6));
  else if (L<intervals[1])
    result = v1 * (POW(R,6) - 320.0*POW(L,3)*POW(R,3) + 720.0*POW(L,4)*POW(R,2) - 576*POW(L,5)*R + 160.0*POW(L,6)) / (60.0*M_PI*POW(L,3)*POW(R,6));
  else
    result = -1.0*v1 / (4*M_PI*POW(L,3));

#elif BSPLINE_DEGREE == 5

  if (L < intervals[0])
    result = -125.0*v1 * (161.0*POW(R,4) - 630.0*POW(L,2)*POW(R,2) + 1125.0*POW(L,4)) / (5376.0*M_PI*POW(R,7));
  else if (L < intervals[1])
    result = v1 * (POW(R,7) + 96250.0*POW(L,3)*POW(R,4) + 65625.0*POW(L,4)*POW(R,3) - 787500.0*POW(L,5)*POW(R,2) + 1093750.0*POW(L,6)*R - 468750.0*POW(L,7)) / (26880.0*M_PI*POW(L,3)*POW(R,7));
  else if (L < intervals[2])
    result = v1 * (437.0*POW(R,7) - 109375.0*POW(L,3)*POW(R,4) + 328125.0*POW(L,4)*POW(R,3) - 393750.0*POW(L,5)*POW(R,2) + 218750.0*POW(L,6)*R - 46875.0*POW(L,7)) / (10752.0*M_PI*POW(L,3)*POW(R,7));
  else
    result = -1.0*v1 / (4*M_PI*POW(L,3));

#elif BSPLINE_DEGREE == 6

  if (L < intervals[0])
    result = 9.0*v1 * (616.0*POW(R,5) - 3024.0*POW(L,2)*POW(R,3) + 9720.0*POW(L,4)*R - 8505*POW(L,5)) / (1120.0*M_PI*POW(R,8));
  else if (L < intervals[1])
    result = -1.0*v1 * (5.0*POW(R,8) + 25704*POW(L,3)*POW(R,8) + 85050.0*POW(L,4)*POW(R,4) - 571536.0*POW(L,5)*POW(R,3) + 1020600*POW(L,6)*POW(R,2) - 787320.0*POW(L,7)*R + 229635.0*POW(L,8)) / (6720.0*M_PI*POW(L,3)*POW(R,8));
  else if (L < intervals[2])
    result = v1 * (169.0*POW(R,8) - 40824.0*POW(L,3)*POW(R,5) + 153090.0*POW(L,4)*POW(R,4) - 244944.0*POW(L,5)*POW(R,3) + 204120.0*POW(L,6)*POW(R,2) - 87480.0*POW(L,7)*R + 15309.0*POW(L,8)) / (2240.0*M_PI*POW(L,3)*POW(R,8));
  else
    result = -1.0*v1 / (4*M_PI*POW(L,3));

#endif

  return result;
}