C+
C NAME:
C	nrRatInt
C PURPOSE:
C	
C CATEGORY:
C	Math: interpolation and extrapolation
C CALLING SEQUENCE:
	subroutine nrRatInt(N,XA,YA,X,Y,DY)
C INPUTS:
C	N	integer
C	XA
C	N	integer		# points
C	XA(N)	real		x-values with known fnc-values
C	YA(N)	real		known fnc-values in XA
C	X	real		x-value where fnc-value is required
C OUTPUTS:
C	Y	real		interpolated fnc-value
C	DY	real		error estimate
C CALLS:
C	Say
C PROCEDURE:
C       See Numerical Recipes, par. 3.2, p. 85
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    integer	N
	    real	XA(N)
	    real	YA(N)
	    real	X
	    real	Y
	    real	DY

	parameter	(TINY = 1.e-25)
	parameter       (NMAX = 20)
	real		C(NMAX)
	real		D(NMAX)
	character	cFlt2Str*14

	if (N .gt. NMAX) call Say('nrRatInt','E','Stop','Dimension XA must be less/equal NMAX')

        NS = 1
        HH = abs(X-XA(1))
        do I=1,N
            H = abs(X-XA(I))
            if (H .eq. 0) then
                Y = YA(I)
                DY = 0.
                return
            else if (H .lt. HH) then
                NS = I
                HH = H
            end if
            C(I) = YA(I)
            D(I) = YA(I)+TINY
        end do
        Y = YA(NS)
        NS = NS-1
        do M=1,N-1
            do I=1,N-M
                W = C(I+1)-D(I)
                H = XA(I+M)-X
                T = (XA(I)-X)*D(I)/H
                DD = T-C(I+1)
                if (DD .eq. 0) then
		    call Say('nrRatInt','W','Pole','function has pole at '//cFlt2Str(X,4))
		    Y = 1.e30
		    return
		end if
                DD = W/DD
                D(I) = C(I+1)*DD
                C(I) = T*DD
            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