// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.

#include "legendrePolynomials.h"



LegendrePolynomials::LegendrePolynomials(int o): n(o) {}



LegendrePolynomials::~LegendrePolynomials() {}



void LegendrePolynomials::f(double u, double *val) const
{
  val[0] = 1;

  for (int i=0;i<n;i++) {

    double a1i = i+1;
    double a3i = 2.*i+1;
    double a4i = i;

    val[i+1] = a3i * u * val[i];
    if (i>0) val[i+1] -= a4i * val[i-1];
    val[i+1] /= a1i;
  }
}



void LegendrePolynomials::fc(double u, double *val) const
{
  f(u, val);
  for (int i = 2; i < n+1; ++i) {
    if (i % 2) val[i] -= u;
    else val[i] -= 1;
  }
}



void LegendrePolynomials::df(double u, double *val) const
{

  // Indeterminate form for u == -1 and u == 1
  if ((u == 1.) || (u == -1.)) {

    for (int k=0;k<=n;k++) val[k] = 0.5*k*(k+1);
    if ((u == -1.) && (n >= 2)) for (int k=2;k<=n;k+=2) val[k] = -val[k];

    return;

  }

  // Now general case

  // Values of the Legendre polynomials
  std::vector<double> tmp(n+1);
  f(u,&(tmp[0]));

  // First value of the derivative
  val[0] = 0;
  double g2 = (1.-u*u);

  // Values of the derivative for orders > 1 computed from the values of the polynomials
  for (int i=1;i<=n;i++) {
    double g1 = - u*i;
    double g0 = (double) i;
    val[i] = (g1 * tmp[i] + g0 * tmp[i-1])/g2;
  }

}