C+
C NAME:
C	nrInterpol
C PURPOSE:
C	Find approximate function values at positions X, based on some
C	interpolation scheme on positions XA and fnc-values YA
C CATEGORY:
C	Math: interpolation
C CALLING SEQUENCE:
	subroutine nrInterpol(N,XA,YA, M,X,Y,Y2A, ID_IN)
C INPUTS:
C	N	integer		# XA-positions
C	XA(N)	real		x-positions where fnc-values are known
C				the array XA MUST be sorted (this can be
C				enforced by the value of ID)
C	YA(N)	real		fnc-values in positions X
C	M	integer		# X-positions
C	X	real		x-positions where fnc-values are required
C				(does not have to be sorted, but it does
C				speed things up)
C	Y2A(N)	real		scratch array
C	ID	integer		determines the type of interpolation
C
C		If ID is negative, the arrays XA,YA are explicitly sorted to
C		ensure that XA is in ascending order
C
C		abs(ID) = 0,1	linear interpolation
C			= 2	for each point in X the interpolation is based
C				on a quadratic polynomial based on three
C				points in XA (piecewise quadratic)
C			= 3	for each point in X the interpolation is base
C				on a cubic polynomial based on four points in XA
C				(piecewise cubic)
C			= 4	natural cubic spline interpolation is used
C			= 5	polynomial interpolation is used (i.e. the
C				unique polynomial of degree n_elements(XA)-1
C OUTPUTS:
C	X	real		may have been sorted (if ID negative)
C	Y	real		interpolated fnc-values in positions X
C CALLS:
C	Sort2R4, nrHunt, nrRatInt, nrPolInt, nrSpline, nrSplInt
C SIDE EFFECTS:
C	Extrapolation is possible but will produce funny results.
C PROCEDURE:
C >	All interpolation procedures are based on chapter 3 of Numerical
C	Recipes (see nrSplInt,nrSpline,nrPolInt for more details)
C >	If the ID value is out of range then ID=1 (linear interpolation,
C	no sorting is assumed)
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    integer	N
	    real	XA(N)
	    real	YA(N)
	    integer	M
	    real	X(M)
	    real	Y(M)
	    real	Y2A(N)
	    integer	ID_IN

	ID = ID_IN
	if (ID .lt. 0) then
	    ID = -ID
	    call Sort2R4(1,N,1,N,XA,YA)
	end if

	if (ID .eq. 9) then
	    do I=1,M
		call nrRatInt(N,XA,YA, X(I),Y(I),DY)
	    end do

	else if (ID .eq. 8) then		! Piecewise cubic
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		if (ILO .le. 1) then
		    call nrRatInt(3,XA(1),YA(1), X(I),Y(I),DY)
		else if (ILO .ge. N-1) then
 		    call nrRatInt(3,XA(N-2),YA(N-2), X(I),Y(I),DY)
		else
		    call nrRatInt(4,XA(ILO-1),YA(ILO-1), X(I),Y(I),DY)
		end if
	    end do

	else if (ID .eq. 7) then		! Piecewise quadratic
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		ILO = min(ILO,N-1)
		if (ILO .lt. N-1) call nrRatInt(3,XA(ILO),YA(ILO), X(I),Yforward,DY)
		if (ILO .gt. 1) call nrRatInt(3,XA(ILO-1),YA(ILO-1), X(I),Ybackward,DY)

		if (ILO .eq. 1) then
		    Y(i) = Yforward
		else if (ILO .eq. N-1) then
		    Y(i) = Ybackward
		else
		    Y(i) = 0.5*(Yforward+Ybackward)
		end if
	    end do

	else if (ID .eq. 6) then		! Linear rational interpolation
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		call nrRatInt(2,XA(ILO),YA(ILO), X(I),Y(I),DY)
	    end do

	else if (ID .eq. 5) then		! Mth order polynomial
	    do I=1,M
		call nrPolInt(N,XA,YA, X(I),Y(I),DY)
	    end do

	else if (ID .eq. 4) then		! Natural cubic spline

	    call nrSpline(N,XA,YA,1.e30,1.e30,Y2A)
	    do I=1,M
		call nrSplInt (N,XA,YA,Y2A, X(I),Y(I))
	    end do

	else if (ID .eq. 3) then		! Piecewise cubic
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		if (ILO .le. 1) then
		    call nrPolInt(3,XA(1),YA(1), X(I),Y(I),DY)
		else if (ILO .ge. N-1) then
 		    call nrPolInt(3,XA(N-2),YA(N-2), X(I),Y(I),DY)
		else
		    call nrPolInt(4,XA(ILO-1),YA(ILO-1), X(I),Y(I),DY)
		end if
	    end do

	else if (ID .eq. 2) then		! Piecewise quadratic
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		ILO = min(ILO,N-1)
		if (ILO .lt. N-1) call nrPolInt(3,XA(ILO),YA(ILO), X(I),Yforward,DY)
		if (ILO .gt. 1) call nrPolInt(3,XA(ILO-1),YA(ILO-1), X(I),Ybackward,DY)

		if (ILO .eq. 1) then
		    Y(i) = Yforward
		else if (ILO .eq. N-1) then
		    Y(i) = Ybackward
		else
		    Y(i) = 0.5*(Yforward+Ybackward)
		end if
	    end do

	else					! Linear interpolation
	    ILO = 1
	    do I=1,M
		call nrHunt(N,XA,X(I),ILO)
		call nrPolInt(2,XA(ILO),YA(ILO), X(I),Y(I),DY)
	    end do
	end if

	return
	end