! Laphi.f90 Lagrange Phi functions, orthogonal polynomials
!         function prototype naming convention:
!           phi            basic Lagrange polynomial
!           phip           first derivative "p do prime"
!           phipp          second derivative,  etc
!
! example  phip(x, i, n, xg()) is
!                     ^__ nth-1 degree
!                  ^__ ith polynomial=set
!             ^___first derivative "prime"
!  xg() is x  grid values, not necessarily uniformly spaced.
!  xg(1) is minimum, left boundary.
!  xg(n) is maximum, right boundary.
!  All grid coordinates must be=increasing order.

module laphi

public :: phi
public :: phip
public :: phipp
public :: phippp
public :: phipppp

contains
  ! phi
  function phi(x, j, n, xg) result(y)
    implicit none
    integer, intent(in) :: j, n
    double precision,  intent(in) :: x
    double precision,  intent(in), dimension(*) :: xg
    double precision :: y, denom
    integer :: i

    denom  = 1.0
    do i=1,n
      if(i .ne. j) then
        denom = denom * (xg(j)-xg(i))
      end if
    end do
    ! evaluate function, n,j,x,xg,denom
    y = 1.0 ! accumulate product
    do i=1,n
      if(i .ne. j) then
        y = y * (x-xg(i))
      end if
    end do
    y = y/denom
  end function phi

  ! first derivative
  function phip(x, j, n, xg) result(yp)
    implicit none
    integer, intent(in) :: j, n
    double precision,  intent(in) :: x
    double precision,  intent(in), dimension(*) :: xg
    double precision :: yp, prod, denom
    integer :: i, ii

    denom = 1.0
    do i=1,n
      if(i .ne. j) then
        denom = denom * (xg(j)-xg(i))
      end if
    end do
    ! evaluate first derivative of function, n,j,x,xg,denom
    yp = 0.0 ! accumulate sum of products
    do i=1,n
      if(i .ne. j) then
        prod = 1.0
        do ii=1,n
          if(ii.ne.i .and. ii.ne.j) then
            prod = prod * (x-xg(ii))
          end if
        end do
        yp = yp + prod
      end if
    end do
    yp = yp/denom ! have first derivative
  end function phip

  ! second derivative 
  function phipp(x, j, n, xg) result(ypp)
    implicit none
    integer, intent(in) :: j, n
    double precision,  intent(in) :: x
    double precision,  intent(in), dimension(*) :: xg
    double precision :: ypp, prod, denom
    integer :: i, ii, iii

    denom = 1.0
    do i=1,n
      if(i .ne. j) then
        denom = denom * (xg(j)-xg(i))
      end if
    end do
    ! evaluate second derivative of function, n,j,x,xg,denom
    ypp = 0.0 ! accumulate sum of products
    do i=1,n
      if(i .ne. j) then
        do ii=1,n
          if(ii.ne.i .and. ii.ne.j) then
            prod = 1.0
            do iii=1,n
              if(iii.ne.i .and. iii.ne.ii .and. iii.ne.j) then
                prod = prod * (x-xg(iii))
              end if
            end do
            ypp = ypp + prod
          end if
        end do
      end if
    end do
    ypp = ypp/denom ! have second derivative
  end function phipp

  ! third derivative
  function phippp(x, j, n, xg) result(yppp)
    implicit none
    integer, intent(in) :: j, n
    double precision,  intent(in) :: x
    double precision,  intent(in), dimension(*) :: xg
    double precision :: yppp, prod, denom
    integer :: i, ii, iii, iiii

    denom = 1.0
    do i=1,n
      if(i .ne. j) then
        denom = denom * (xg(j)-xg(i))
      end if
    end do
    ! evaluate third derivative of function, n,j,x,xg,denom
    yppp = 0.0 ! accumulate sum of products
    do i=1,n
      if(i .ne. j) then
        do ii=1,n
          if(ii.ne.i .and. ii.ne.j) then
            do iii=1,n
              if(iii.ne.i .and. iii.ne.ii .and. iii.ne.j) then
                prod = 1.0
                do iiii=1,n
                  if(iiii.ne.i .and. iiii.ne.ii .and. iiii.ne.iii .and. &
                     iiii.ne.j) then
                    prod = prod * (x-xg(iiii))
                  end if
                end do
                yppp = yppp + prod
              end if
            end do
          end if
        end do
      end if
    end do
    yppp = yppp/denom  ! have third derivative
  end function phippp

  ! fourth derivative
  function phipppp(x, j, n, xg) result(ypppp)
    implicit none
    integer, intent(in) :: j, n
    double precision,  intent(in) :: x
    double precision,  intent(in), dimension(*) :: xg
    double precision :: ypppp, prod, denom
    integer :: i, ii, iii, iiii, iiiii

    denom = 1.0
    do i=1,n
      if(i .ne. j) then
        denom = denom * (xg(j)-xg(i))
      end if
    end do
    ! evaluate fourth derivative of function, n,j,x,xg,denom
    ypppp = 0.0 ! accumulate sum of products
    do i=1,n
      if(i .ne. j) then
        do ii=1,n
          if(ii.ne.i .and. ii.ne.j) then
            do iii=1,n
              if(iii.ne.i .and. iii.ne.ii .and. iii.ne.j) then
                do iiii=1,n
                  if(iiii.ne.i .and. iiii.ne.ii .and. iiii.ne.iii .and. &
                     iiii.ne.j) then
                    prod = 1.0
                    do iiiii=1,n
                      if(iiiii.ne.i .and. iiiii.ne.ii .and. iiiii.ne.iii &
                         .and. iiiii.ne.iiii .and. iiiii.ne.j) then
                        prod = prod * (x-xg(iiiii))
                      end if
                    end do
                    ypppp = ypppp + prod
                  end if
                end do
              end if
            end do
          end if
        end do
      end if
    end do
    ypppp = ypppp/denom ! have fourth derivative
  end function phipppp

end module laphi