// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */ 
/* License as published by the Free Software Foundation either */ 
/* version 3.0 of the License, or (at your option) any later version. */

/* This software is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/. */

#include "JacobiPoly.hpp"

double JacobiPoly::
evalJacobi(double a, double b, unsigned int n, double x)
{
 if (n == 0)
   return 1.0;
 else if (n == 1)
   return 0.5 * ( a - b + ( a + b + 2.0 ) * x );
 else
 {
   double apb = a + b;
   double jn_2 = 1.0;
   double jn_1 = 0.5 * ( a - b + ( apb + 2.0 ) * x );
   double j = 0.0;
   for (unsigned int k = 2; k < n+1; k++)
   {
      j = (((2.0*k+apb-1.0)*(a*a-b*b) + (2.0*k+apb-2.0)*(2.0*k+apb-1.0)*(2.0*k+apb)*x) * jn_1 - (2.0*(k-1.0+a)*(k-1.0+b)*(2.0*k+apb)) * jn_2)
                                                                                                          /(2.0*(k)*(k+apb)*(2.0*k+apb-2.0));
      jn_2 = jn_1;
      jn_1 = j;
    }
    return j;
  }
}

double JacobiPoly::
evalJacobiDeriv(double a, double b, unsigned int n, double x) 
{
  if (n == 0)
    return 0.0;
  else 
    return 0.5 * (a + b + n + 1) * evalJacobi(a+1, b+1, n-1, x);
}

double JacobiPoly::
evalJacobiSecondDeriv(double a, double b, unsigned int n, double x) 
{
  if (n == 0)
    return 0.0;
  else 
    return ((a - b + ((a + b + 2) * x))/(1 - x * x)) * evalJacobiDeriv(a, b, n, x);
}

void JacobiPoly::
legendreGaussQpointsWeights (unsigned int n, RealVector& qPoints, RealVector& qWeights)  
{
  qPoints.resize(n);
  qWeights.resize(n);
  int m,j,i;
  double z1,z,pp,p3,p2,p1;
  m=(n+1)/2;
  for (i=1;i<=m;i++) 
  {
    z=cos(3.141592654*(i-0.25)/(n+0.5));
    do 
    {
      p1=1.0;
      p2=0.0;
      for (j=1;j<=n;j++) 
      {
        p3=p2;
        p2=p1; 
	p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
      }
      pp=n*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
    } while (fabs(z-z1) > 1.e-12); 
    
    qPoints(i-1)=-z; 
    qPoints(n-i)=z; 
  
    qWeights(i-1)=2.0/((1.0-z*z)*pp*pp); 
    qWeights(n-i)=qWeights(i-1);
  }
}

void JacobiPoly::
jacobiGaussQpoints (double a, double b, unsigned int n, RealVector& qPoints)  
{
  double tolerance = 1e-12;
  unsigned int maxIter = 100;
  qPoints.resize(n);
  for (unsigned int k = 0; k < n; k++)
  {
    double guess = -cos((2.0*k+1.0) * 3.1415926535897932384626433832795029L / (2.0 * n));
    if (k>0)
      guess = 0.5 * (guess + qPoints(k-1));
    for (unsigned int j = 0; j <= maxIter; j++)
    {
      double s = 0.0;
      for (unsigned int i = 0; i < k; i++)
      {
         s = s + 1.0 / (guess - qPoints(i));
      }
      double jac = evalJacobi(a,b,n,guess);
      double jacDeriv = evalJacobiDeriv(a,b,n,guess);
      double delta = jac / (jacDeriv - jac * s);
      guess -= delta;
      if (fabs(delta) < tolerance)
      {
         qPoints(k) = guess;
         break;
      }
      assert (j != maxIter);
    }         
  }
}

