C+
C NAME:
C	nrPolInt
C PURPOSE:
C	Function value approximation by 1-dimensional polynomial interpolation
C CATEGORY:
C	Math: polynomial interpolation
C CALLING SEQUENCE:
	subroutine nrPolInt(N,XA,YA,X,Y,DY)
C INPUTS:
C	N	integer		order of interpolating polynomial
C				max. permitted value is NMAX=10
C	XA(N)	real		x-coordinates
C	YA(N)	real		y-coordinates (fnc-values)
C	X	real		x-coordinate (where interpolated fnc value is required)
C OUTPUTS:
C	Y	real		y-coordinate (interpolated fnc-value)
C	DY	real		error estimate of Y
C CALLS:
C	Say
C RESTRICTIONS:
C	The maximum order is set by parameter NMAX (=10)
C PROCEDURE:
C	See Numerical Recipes, par. 3.1, p. 81 (Neville's algorithm)
C MODIFICATION HISTORY:
C	JUN-1992, Paul Hick (UCSD)
C-
	    integer	N
	    real	XA(N)
	    real	YA(N)
	    real	X
	    real	Y
	    real	DY

	parameter	(NMAX=100)
	real		C(NMAX)
	real		D(NMAX)

	if (N .gt. NMAX) call Say('nrPolInt','E','Stop','maximum dimension XA,YA is NMAX=100')

	NS = 1
	DIF = abs(X-XA(1))
	do I=1,N
	    DIFT = abs(X-XA(I))
	    if (DIFT .lt. DIF) then
		NS = I
		DIF = DIFT
	    end if
	    C(I) = YA(I)
	    D(I) = YA(I)
	end do

	Y = YA(NS)
	NS = NS-1

	do M=1,N-1
	    do I=1,N-M
		HO = XA(I)-X
		HP = XA(I+M)-X
		W = C(I+1)-D(I)
		DEN = HO-HP
		if (DEN .eq. 0) call Say('nrPolInt','E','Stop','identical input XA values')
		DEN = W/DEN
		D(I) = HP*DEN
		C(I) = HO*DEN
	    end do
	    if (2*NS .lt. N-M) then
		DY = C(NS+1)
	    else
		DY = D(NS)
		NS = NS-1
	    end if
	    Y = Y+DY
	end do

	return
	end