// laphi.java  Lagrange Phi functions, orthogonal polynomials    
//         function prototype naming convention:              
//           phi            basic Lagrange polynomial         
//           phip           first derivative "p for prime"    
//           phipp          second derivative,                
//           phippp         third derivative,                 
//           phipppp        fourth derivative,                
//                                                            
// example  phip(x, i, n, xg[]) is                            
//                        ^__ xg[0] through xg[n] "zeros"     
//                     ^__ nth degree                         
//                  ^__ ith polynomial in set  0 <= i <=n     
//             ^___first derivative "prime"                   
//  the ith polynomial is 1.0 at xg[i] and zero at xg[j] j!=i 
//  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                    

public class laphi
{ 

// 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 phipppp
} // end class laphi
// end laphi.java