//
// CDDL HEADER START
//
// The contents of this file are subject to the terms of the Common Development
// and Distribution License Version 1.0 (the "License").
//
// You can obtain a copy of the license at
// http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
// specific language governing permissions and limitations under the License.
//
// When distributing Covered Code, include this CDDL HEADER in each file and
// include the License file in a prominent location with the name LICENSE.CDDL.
// If applicable, add the following below this CDDL HEADER, with the fields
// enclosed by brackets "[]" replaced with your own identifying information:
//
// Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
//
// CDDL HEADER END
//

//
// Copyright (c) 2014, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Ryan S. Elliott
//    Mingjian Wen
//

#include "EAM_QuinticHermiteSpline.hpp"
#include "EAM_Implementation.hpp"

//************************************************************************************
void EAM_Implementation::SplineInterpolate(double const * const dat,
                                           double const delta,
                                           int const n,
                                           double * const coe)
{  // setup convenient pointers (spline) into the coefficients (coe) array
  double ** const spline = new double *[n];  // deleted at end of function
  for (int i = 0; i < n; ++i) { spline[i] = &coe[i * NUMBER_SPLINE_COEFF]; }

  for (int m = 0; m < n; ++m) { spline[m][F_CONSTANT] = dat[m]; }

  // linear term and quadratic term are estimated using 7 points finite
  // difference
  // linear coeff
  spline[0][F_LINEAR]
      = (-11 * dat[0] + 18 * dat[1] - 9 * dat[2] + 2 * dat[3]) / 6;
  spline[1][F_LINEAR]
      = (-3 * dat[0] - 10 * dat[1] + 18 * dat[2] - 6 * dat[3] + dat[4]) / 12;
  spline[2][F_LINEAR] = dat[0] / 20 - dat[1] / 2 - dat[2] / 3 + dat[3]
                        - dat[4] / 4 + dat[5] / 30;
  spline[n - 3][F_LINEAR] = -dat[n - 6] / 30 + dat[n - 5] / 4 - dat[n - 4]
                            + dat[n - 3] / 3 + dat[n - 2] / 2 - dat[n - 1] / 20;
  spline[n - 2][F_LINEAR] = (-dat[n - 5] + 6 * dat[n - 4] - 18 * dat[n - 3]
                             + 10 * dat[n - 2] + 3 * dat[n - 1])
                            / 12;
  spline[n - 1][F_LINEAR]
      = (-2 * dat[n - 4] + 9 * dat[n - 3] - 18 * dat[n - 2] + 11 * dat[n - 1])
        / 6;

  for (int m = 3; m < n - 3; ++m)
  {
    spline[m][F_LINEAR] = -dat[m - 3] / 60 + 3 * dat[m - 2] / 20
                          - 3 * dat[m - 1] / 4 + 3 * dat[m + 1] / 4
                          - 3 * dat[m + 2] / 20 + dat[m + 3] / 60;
  }

  // quadratic coeff
  spline[0][F_QUADRATIC] = (2 * dat[0] - 5 * dat[1] + 4 * dat[2] - dat[3]) / 2;
  spline[1][F_QUADRATIC]
      = ((11 * dat[0] - 20 * dat[1] + 6 * dat[2] + 4 * dat[3] - dat[4]) / 12)
        / 2;
  spline[2][F_QUADRATIC] = (-dat[0] / 12 + 4 * dat[1] / 3 - 5 * dat[2] / 2
                            + 4 * dat[3] / 3 - dat[4] / 12)
                           / 2;
  spline[n - 3][F_QUADRATIC]
      = (-dat[n - 5] / 12 + 4 * dat[n - 4] / 3 - 5 * dat[n - 3] / 2
         + 4 * dat[n - 2] / 3 - dat[n - 1] / 12)
        / 2;
  spline[n - 2][F_QUADRATIC] = ((-dat[n - 5] + 4 * dat[n - 4] + 6 * dat[n - 3]
                                 - 20 * dat[n - 2] + 11 * dat[n - 1])
                                / 12)
                               / 2;
  spline[n - 1][F_QUADRATIC]
      = (-dat[n - 4] + 4 * dat[n - 3] - 5 * dat[n - 2] + 2 * dat[n - 1]) / 2;

  for (int m = 3; m < n - 3; ++m)
  {
    spline[m][F_QUADRATIC]
        = (dat[m - 3] / 90 - 3 * dat[m - 2] / 20 + 3 * dat[m - 1] / 2
           - 49 * dat[m] / 18 + 3 * dat[m + 1] / 2 - 3 * dat[m + 2] / 20
           + dat[m + 3] / 90)
          / 2;
  }

  // cubic, quartic, and quintic coeff
  for (int m = 0; m < n - 1; ++m)
  {
    spline[m][F_CUBIC]
        = 10 * spline[m + 1][F_CONSTANT] - 4 * spline[m + 1][F_LINEAR]
          + spline[m + 1][F_QUADRATIC] - 10 * spline[m][F_CONSTANT]
          - 6 * spline[m][F_LINEAR] - 3 * spline[m][F_QUADRATIC];

    spline[m][F_QUARTIC]
        = -15 * spline[m + 1][F_CONSTANT] + 7 * spline[m + 1][F_LINEAR]
          - 2 * spline[m + 1][F_QUADRATIC] + 15 * spline[m][F_CONSTANT]
          + 8 * spline[m][F_LINEAR] + 3 * spline[m][F_QUADRATIC];

    spline[m][F_QUINTIC]
        = 6 * spline[m + 1][F_CONSTANT] - 3 * spline[m + 1][F_LINEAR]
          + spline[m + 1][F_QUADRATIC] - 6 * spline[m][F_CONSTANT]
          - 3 * spline[m][F_LINEAR] - spline[m][F_QUADRATIC];
  }
  spline[n - 1][F_CUBIC] = 0;
  spline[n - 1][F_QUARTIC] = 0;
  spline[n - 1][F_QUINTIC] = 0;

  // derivatives
  for (int m = 0; m < n - 1; ++m)
  {
    spline[m][DF_CONSTANT] = spline[m][F_LINEAR] / delta;
    spline[m][DF_LINEAR] = 2.0 * spline[m][F_QUADRATIC] / delta;
    spline[m][DF_QUADRATIC] = 3.0 * spline[m][F_CUBIC] / delta;
    spline[m][DF_CUBIC] = 4.0 * spline[m][F_QUARTIC] / delta;
    spline[m][DF_QUARTIC] = 5.0 * spline[m][F_QUINTIC] / delta;
  }

  for (int m = 0; m < n - 1; ++m)
  {
    spline[m][D2F_CONSTANT] = spline[m][DF_LINEAR] / delta;
    spline[m][D2F_LINEAR] = 2.0 * spline[m][DF_QUADRATIC] / delta;
    spline[m][D2F_QUADRATIC] = 3.0 * spline[m][DF_CUBIC] / delta;
    spline[m][D2F_CUBIC] = 4.0 * spline[m][DF_QUARTIC] / delta;
  }

  delete[] spline;
}