C
C PURPOSE:
C       Interpolate the function f(x) at the point xtarget. That is given the
C       set of ordinal points (x1,y1), (x2,y2), (x3,y3), ..., (xn,yn) this
C       routine will return the value f(xtarget), where xtarget lies in the
C       range (x1,yx) --> (xn,yn). To accoplish this goal polynomial
C       interpolation is used via Neville's algorithm.
C
C INPUTS:
C       porder  = integer, number of points used in interpolation
C       xvals   = real, dimension(porder), ordinal component of interpolation
C                 nodes
C       yvals   = real, dimension(porder), range of inerpolation nodes
C       xtarget = real, x location of interpolation
C
C OUTPUTS:
C       dy      = real, error estimate
C       yinterp = real, y value of interpolation (the answer)
C
C DEPENDENCIES:
C       None
C
C NOTES:
C       This routine was addapted (taken) from the second edition of 'Numerical
C       Recipes in Fortran 77, The Art of Scientific Computing', by Press,
C       Teukolsky, Vetterling, and Flannery, copyright 1992
C
C       VARIABLES:
C               ccor    = real, dim(porder), c correction values
C               dcor    = real, dim(porder), d correction values
C               nsindex = integer, index of the closest table entry
C
C DATE:
C       February 13, 2020 @ 1:03 PM
C
C AUTHOR:
C       Luke Cota
C

      subroutine poly_interpolate(porder, xvals, yvals, xtarget, dy_out,
     &                            yinterp_out)
        implicit none
        integer :: porder, icur, mcur, nsindex
        real    :: xtarget, yinterp, diff, difft, den
        real    :: dy, ho, hp, weight
        real    :: xvals(porder), yvals(porder)
        real    :: ccor(porder), dcor(porder)
C        double precision :: xtarget, yinterp, diff, difft, den
C        double precision :: dy, ho, hp, weight
C        double precision :: xvals(porder), yvals(porder)
C        double precision :: ccor(porder), dcor(porder)
        real :: yinterp_out, dy_out

        nsindex = 1
        diff = abs(xtarget - xvals(1))

        do icur = 1, porder
          difft = abs(xtarget - xvals(icur))
          if (difft .lt. diff) then
            nsindex = icur
            diff = difft
          end if
          ccor(icur) = yvals(icur)
          dcor(icur) = yvals(icur)
        end do

        yinterp = yvals(nsindex)
        nsindex = nsindex - 1

        do mcur = 1, porder - 1
          do icur = 1, porder - mcur
            ho = xvals(icur) - xtarget
            hp = xvals(icur + mcur) - xtarget
            weight = ccor(icur + 1) - dcor(icur)
            den = ho - hp
            if (den .eq. 0.0) then
              print *, 'Error in poly_interpolate routine: '
              print *, 'den = 0, this indicates two of the input'
              print *, 'nodes are the same.'
              stop
            end if
            den = weight / den
            dcor(icur) = hp * den
            ccor(icur) = ho * den
          end do
          if (2 * nsindex .lt. porder - mcur) then
            dy = ccor(nsindex + 1)
          else
            dy = dcor(nsindex)
            nsindex = nsindex - 1
          end if
          yinterp = yinterp + dy
        end do

        yinterp_out = real(yinterp)
        dy_out = real(dy_out)
        return
      end subroutine poly_interpolate