;+ 
;NAME: 
;	INTERP.PRO 
; 
;LAST CHANGED: ----------------------------------- 15-nov-00
; 
;PURPOSE: 
;	Interpolates values
;	may be obsolete and might be eliminated
; 
;CATEGORY: 
;	STAND ALONE: utilities 
; 
;CALLING SEQUENCE: 
;	Result = INTERP(Xi,XX,YY,MISSING=MI) 
; 
;CALLED FROM AXIS: 
;	many procedures 
; 
;ROUTINES 
;	function INTERP, Xi, XX, YY, MISSING=MI
;	interpolates values
; 
;COMMON BLOCKS: none 
; 
;MODIFICATION HISTORY: 
; (08-nov-00 cgz) AXIS standard documentation 
; (08-nov-00 cgz) migrated from MATH.PRO
; (15-nov-00 cgz) Replaced call of NCOLS function by explicit (SIZE(var))(1)
;- 
; ******************************************************************************* 

;like igor's interp, returns yy(xx) interpolated at xi
;if xi is a constant (array), returns a constant (array)
;requires: xx should be monotonic increasing, xi doesn't have to be

function interp,xi, xx, yy, missing=mi

; is this function ever used as part of AXIS ?????
; let's find out..... if this is ever called, create a file
print, string(7B),'INTERP.PRO is indeed used',string(7B) 
get_lun,lun 
openw,lun,'keep_interp.txt' 
printf,lun,'INTERP.PRO is indeed used' 
printf,lun,'So you can NOT delete it' 
close,lun 
free_lun,lun 

  on_error,2
;convert xi to an array if it isn't one
  si=size(xi)
  if si(0) eq 0 then a=[xi] else a=xi
  a=a+0.0; convert to floating
  si=size(a)
  mn=min(xx)
  mx=max(xx)
  ;convert xi values to index values
  iout=a
  j=0
  for i=0,si(1)-1 do begin
    if (a(i) ge mn) and (a(i) le mx) then begin
      while a(i) gt xx(j+1) do j=j+1
      iout(i)=j+(a(i)-xx(j))/(xx(j+1)-xx(j)) 
    endif else begin ;out of range
;      if (a(i) lt mn) then iout(i)=-1 else iout(i)=ncols(xx)
      if (a(i) lt mn) then iout(i)=-1 else iout(i)=(size(xx))(1)
    endelse
    if i lt si(1)-1 then if a(i+1) lt a(i) then j=0
  endfor
  ;if si(1) eq 1 then return,iout(0) else return,iout
  if keyword_set(mi) then $
    return,interpolate(yy, iout,missing=mi) $
  else $
    return,interpolate(yy, iout)

end

