function NR_INTER, XA,YA, X,			$
		cspline=CSPLINE,		$
		polynom=POLYNOM,		$
		cubic=CUBIC, quadratic=QUADRATIC, $
		sortxa=SORTXA, plt=plt, oplt=oplt

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	NR_INTER
; PURPOSE:
;	Find approximate function values at positions X, based on some
;	interpolation scheme on positions XA and fnc-values YA
; CATEGORY:
;	Math: interpolation
; CALLING SEQUENCE:
;	Y = NR_INTER( XA,YA, X, /cspline, /quadratic, /cubic, /polynom,
;		/sortxa, /plt, /oplt)
; INPUTS:
;	XA	array		x-positions where fnc-values are known
;				the array XA MUST be sorted (this can be
;				enforced by using the SORTXA keyword
;	YA	array		fnc-values in positions X
;	X	array		x-positions where fnc-values are required
; OPTIONAL INPUT PARAMETERS:
;	CSPLINE		if set, natural cubic spline interpolation is used
;	QUADRATIC	if set, for each point in X the interpolation is based
;			on a quadratic polynomial based on three points in XA
;	CUBIC		if set, for each point in X the interpolation is base
;			on a cubic polynomial based on four points in XA
;	POLYNOM		if set, polynomial interpolation is used (i.e. the
;			unique polynomial of degree n_elements(XA)-1
;	SORTXA		if set, the array XA is explicitly sorted (YA is sorted
;			accordingly)
;	PLT		execute command PLOT,X,Y before returning
;	OPLT		execute command OPLOT,X,Y before returning
; OUTPUTS:
;	Y	array		interpolated fnc-values in positions X
; CALLS:
;	POLINT_ARR, SPLINT_ARR
; PROCEDURE:
; >	All interpolation procedures are based on chapter 3 of Numerical Recipes.
; >	If none of the keyword CPLINE,POLYNOM,CUBIC or QUADRATIC is set
;	linear interpolation is used
; >	Only one the keywords CPLINE,POLYNOM,CUBIC or QUADRATIC is used.
;	If more then one is set, the order of precedence is:
;	CSPLINE, POLYNOM, CUBIC, QUADRATIC
; >	If both the PLT and the OPLT keywords are specified, PLT take precedence
; MODIFICATION HISTORY:
;	JUN-1993, Paul Hick (UCSD)
;-

N = n_elements(XA)			; Check XA and YA dimensions
if N ne n_elements(YA) then message, 'XA and YA must have same dimension

if keyword_set(SORTX) then begin	; Sort XA and YA
    XS = sort(XA)  &  XA = XA(XS)  &  YA = YA(XS)
endif
Y = 0*X  &  M = n_elements(X)

if keyword_set (CSPLINE) then			$	; Natural cubic spline
    Y = splint_arr ( XA,YA, X )		else	$

if keyword_set (POLYNOM) then			$	; Polynomial
    Y = polint_arr ( XA,YA, X )		else	$

if keyword_set (CUBIC) then begin

    for i=0,M-1 do begin
	j = where( X(i) ge XA )
	j = j(n_elements(j)-1)
	j = (j > 0)

	if j eq 0 then begin
	    xs = [XA(0),XA(1),XA(2)]
	    ys = [YA(0),YA(1),YA(2)]
	endif else if j ge N-2 then begin
 	    xs = [XA(N-3),XA(N-2),XA(N-1)]
	    ys = [YA(N-3),YA(N-2),YA(N-1)]
	endif else begin
	    xs = [XA(j-1),XA(j),XA(j+1),XA(j+2)]
	    ys = [YA(j-1),YA(j),YA(j+1),YA(j+2)]
	endelse

	Y(i) = polint_arr( xs,ys, X(i) )

    endfor

endif else if keyword_set (QUADRATIC) then begin

    for i=0,M-1 do begin
	j = where( X(i) ge XA )
	j = j(n_elements(j)-1)
	j = (j > 0) < (N-2)

	if j lt N-2 then begin
	    xs = [XA(j),XA(j+1),XA(j+2)]
	    ys = [YA(j),YA(j+1),YA(j+2)]
	    yforward = polint_arr( xs,ys, X(i) )
	endif
	if j gt 0 then begin
 	    xs = [XA(j-1),XA(j),XA(j+1)]
	    ys = [YA(j-1),YA(j),YA(j+1)]
	    ybackward = polint_arr( xs,ys, X(i) )
	endif

	if j eq   0 then Y(i) = yforward  else	$
	if j eq N-2 then Y(i) = ybackward else	$
	   Y(i) = 0.5*(yforward+ybackward)

    endfor

endif else begin				; Linear interpolation

    for i=0,M-1 do begin
	j = where( X(i) ge XA )
	j = j(n_elements(j)-1)
	j = (j > 0) < (N-2)

 	xs = [XA(j),XA(j+1)]  &  ys = [YA(j),YA(j+1)]
	Y(i) = polint_arr( xs,ys, X(i) )

    endfor

endelse

if keyword_set ( plt) then  plot, X,Y	else	$	; Optional plotting
if keyword_set (oplt) then oplot, X,Y

return, Y  &  end