# laphi.py  Lagrange Phi functions, orthogonal polynomials    
#           function 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                    

import math

def phi(x, j, n, xg):
  denom = 1.0
  for i in range(0,n+1):
    if i != j:
      denom = denom * (xg[j]-xg[i])

  y = 1.0  # accumulate product 
  for i in range (0,n+1):
    if i != j:
      y = y * (x-xg[i])

  return y/denom  # have phi, shape function

# first derivative 
def phip(x, j, n, xg):
  denom = 1.0
  for i in range(0,n+1):
    if i != j:
      denom = denom * (xg[j]-xg[i])

  yp = 0.0  # accumulate sum of products 
  for i in range(0,n+1):
    if i==j:
      continue

    prod = 1.0
    for ii in range(0,n+1):
      if ii==i or ii==j:
        continue

      prod = prod * (x-xg[ii])

    yp = yp + prod

  return yp/denom  # have first derivative 


# second derivative 
def phipp(x, j, n, xg):
  denom = 1.0
  for i in range(0,n+1):
    if i != j:
      denom = denom * (xg[j]-xg[i])

  ypp = 0.0  # accumulate sum of products 
  for i in range(0,n+1):
    if i==j:
      continue

    for ii in range(0,n+1):
      if ii==i or ii==j:
        continue

      prod = 1.0
      for iii in range (0,n+1):
        if iii==i or iii==ii or iii==j:
          continue

        prod = prod * (x-xg[iii])

      ypp = ypp + prod

  return ypp/denom  # have second derivative 


# third derivative 
def phippp(x, j, n, xg):
  denom = 1.0
  for i in range (0,n+1):
    if i != j:
      denom = denom * (xg[j]-xg[i])

  yppp = 0.0  # accumulate sum of products 
  for i in range (0,n+1):
    if i==j:
      continue

    for ii in range (0,n+1):
      if ii==i or ii==j:
        continue

      for iii in range (0,n+1):
        if iii==i or iii==ii or iii==j:
          continue

        prod = 1.0
        for iiii in range (0,n+1):
          if iiii==i or iiii==ii or iiii==iii or iiii==j:
            continue

          prod = prod * (x-xg[iiii])

        yppp = yppp + prod

  return yppp/denom  # have third derivative 


# fourth derivative 
def phipppp(x, j, n, xg):
  denom = 1.0
  for i in range (0,n+1):
    if i != j:
      denom = denom * (xg[j]-xg[i])

  ypppp = 0.0  # accumulate sum of products 
  for i in range (0,n+1):
    if i==j:
      continue

    for ii in range (0,n+1):
      if ii==i or ii==j:
        continue

      for iii in range (0,n+1):
        if iii==i or iii==ii or iii==j:
          continue

        for iiii in range (0,n+1):
          if iiii==i or iiii==ii or iiii==iii or iiii==j:
            continue

          prod = 1.0
          for iiiii in range (0,n+1):
            if iiiii==i or iiiii==ii or iiiii==iii or iiiii==iiii or iiiii==j:
              continue

            prod = prod * (x-xg[iiiii])

          ypppp = ypppp + prod

  return ypppp/denom  # have fourth derivative 

# end laphi.py