/* laphi.c  Lagrange Phi functions, orthogonal polynomials    */
/*         function prototype naming convention:              */
/*           phi            basic Lagrange polynomial         */
/*           phip           first derivative "p for prime"    */
/*           phipp          second derivative,  etc           */
/*                                                            */
/* example  phip(x, i, n, xg[]) is                            */
/*                     ^__ nth degree                         */
/*                  ^__ ith polynomial in set                 */
/*             ^___first derivative "prime"                   */
/*  xg[] is x  grid values, not necessarily uniformly spaced. */
/*  xg[0] is minimum, left boundary.                          */
/*  xg[n] is maximum, right boundary.                         */
/*  All grid coordinates must be in increasing order.         */
/* laphi.c  family of Lagrange functions and derivatives      */ 
/*          L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)   */
/*          for points x_0, x_1, ... , x_n                    */

#include "laphi.h"

/* phi */
double phi(double x, int j, int n, double xg[])
{
  double y, prod;
  double denom = 1.0;
  int i;

  for(i=0; i<=n; i++)
  {
    if(i != j) denom = denom * (xg[j]-xg[i]);
  }

  /* evaluate function, n,j,x,xg,denom */
  y = 1.0; /* accumulate product */
  for(i=0; i<=n; i++)
  {
    if(i != j) y = y * (x-xg[i]);
  }
  y = y/denom;
  return y;
}

/* first derivative */
double phip(double x, int j, int n, double xg[])
{
  double yp, prod;
  double denom = 1.0;
  int i, ii;

  for(i=0; i<=n; i++)
  {
    if(i != j) denom = denom * (xg[j]-xg[i]);
  }
  /* evaluate first derivative of function, n,j,x,xg,denom */
  yp = 0.0; /* accumulate sum of products */
  for(i=0; i<=n; i++)
  {
    if(i==j) continue;
    prod = 1.0;
    for(ii=0; ii<=n; ii++)
    {
      if(ii==i || ii==j) continue;
      prod = prod * (x-xg[ii]);
    }
    yp = yp + prod;
  }
  yp = yp/denom; /* have first derivative */
  return yp;
}

/* second derivative */
double phipp(double x, int j, int n, double xg[])
{
  double ypp, prod;
  double denom = 1.0;
  int i, ii, iii;

  for(i=0; i<=n; i++)
  {
    if(i != j) denom = denom * (xg[j]-xg[i]);
  }
  /* evaluate second derivative of function, n,j,x,xg,denom */
  ypp = 0.0; /* accumulate sum of products */
  for(i=0; i<=n; i++)
  {
    if(i==j) continue;
    for(ii=0; ii<=n; ii++)
    {
      if(ii==i || ii==j) continue;
      prod = 1.0;
      for(iii=0; iii<=n; iii++)
      {
        if(iii==i || iii==ii || iii==j) continue;
        prod = prod * (x-xg[iii]);
      }
      ypp = ypp + prod;
    }
  }
  ypp = ypp/denom; /* have second derivative */
  return ypp;
}

/* third derivative */
double phippp(double x, int j, int n, double xg[])
{
  double yppp, prod;
  double denom = 1.0;
  int i, ii, iii, iiii;

  for(i=0; i<=n; i++)
  {
    if(i != j) denom = denom * (xg[j]-xg[i]);
  }
  /* evaluate third derivative of function, n,j,x,xg,denom */
  yppp = 0.0; /* accumulate sum of products */
  for(i=0; i<=n; i++)
  {
    if(i==j) continue;
    for(ii=0; ii<=n; ii++)
    {
      if(ii==i || ii==j) continue;
      for(iii=0; iii<=n; iii++)
      {
        if(iii==i || iii==ii || iii==j) continue;
        prod = 1.0;
        for(iiii=0; iiii<=n; iiii++)
        {
          if(iiii==i || iiii==ii || iiii==iii || iiii==j) continue;
          prod = prod * (x-xg[iiii]);
        }
        yppp = yppp + prod;
      }
    }
  }
  yppp = yppp/denom;  /* have third derivative */
  return yppp;
}

/* fourth derivative */
double phipppp(double x, int j, int n, double xg[])
{
  double ypppp, prod;
  double denom = 1.0;
  int i, ii, iii, iiii, iiiii;

  for(i=0; i<=n; i++)
  {
    if(i != j) denom = denom * (xg[j]-xg[i]);
  }

  /* evaluate fourth derivative of function, n,j,x,xg,denom */
  ypppp = 0.0; /* accumulate sum of products */
  for(i=0; i<=n; i++)
  {
    if(i==j) continue;
    for(ii=0; ii<=n; ii++)
    {
      if(ii==i || ii==j) continue;
      for(iii=0; iii<=n; iii++)
      {
        if(iii==i || iii==ii || iii==j) continue;
        for(iiii=0; iiii<=n; iiii++)
	{
          if(iiii==i || iiii==ii || iiii==iii || iiii==j) continue;
          prod = 1.0;
          for(iiiii=0; iiiii<=n; iiiii++)
          {
            if(iiiii==i || iiiii==ii || iiiii==iii || iiiii==iiii ||
               iiiii==j) continue;
            prod = prod * (x-xg[iiiii]);
          }
          ypppp = ypppp + prod;
        }
      }
    }
  }
  ypppp = ypppp/denom; /* have fourth derivative */
  return ypppp;
}

/* end laphi.c */