FUNCTION interp_prune_err, ixl, ixn, x, ys, ABS_ERROR=abs_error, $
	PLOT_PRUNE=plot_prune
;
; See interp_pruned_x.pro
;
; - - - - - - - - - - - - -
;
; EXAMPLE of USE
;
;hetgcal> meg = rdb_read(!DDHETGCAL+'/cip/hetgmegD1996-11-01greffN0002.rdb')
;hetgcal> print, interp_prune_err(3500,3800,meg.energy, $
;                                               [[meg.oz],[meg.op1]] ,/PLOT)
;     0.525078     0.763524
; - - - - - - - - - - - - -
; 1999-01-16 dd Initial creation.
; - - - - - - - - - - - - -

; Check some dimensions
ysize = SIZE(ys)
; Number of y tables supplied...
; Y should be 2 dimensional but may be 1-D if
; only a single Y-column is being used.
if ysize(0) EQ 1 then ntabs=1 else ntabs = ysize(2)

; returned values:
max_errors = FLTARR(ntabs)

; Plot what we're doing?!?
if KEYWORD_SET(PLOT_PRUNE) then begin
  !p.multi= [0,1,ntabs]
end

; Check for no points between ixl and ixn - lin interp is exact -
; and ixn LE ixl:
if ixn LE ixl+1 then RETURN, 0.0*FLTARR(ntabs)

; Select the x values that are under consideration (between
; ixl and ixn):
ixs = indgen(ixn-(ixl+1)) + ixl+1  ; starts at ixl+1

; Assign the linear x-fraction for these points:
x_fracs = ( x(ixs) - x(ixl) )/( x(ixn) - x(ixl) )

; OK, let's evaluate the errors for each table

for it=0,ntabs-1 do begin
  y_interps = x_fracs * ( ys(ixn, it) - ys(ixl, it) ) + ys(ixl, it)
  errors = ABS(y_interps - ys(ixs, it))
  if KEYWORD_SET(PLOT_PRUNE) then begin
    plot, x(ixl:ixn), ys(ixl:ixn, it), PSYM=4
    oplot, x(ixl:ixn), [ys(ixl,it),y_interps,ys(ixn,it)], PSYM=-3
  end
  if NOT(KEYWORD_SET(ABS_ERROR)) then begin
    ; relative errors
    ; the error is the abs error divided by the tabulated
    ; value, ydenom.  This has a problem when the tabulated value 
    ; is 0.0 !  Catch and handle these cases:
    ydenom = ABS(ys(ixs, it))
    zeros = where(ydenom EQ 0.0, nzs)
    if nzs GE 1 then begin
      ; set ydenom equal to 1.0 at any place that errors is 0.0 -
      ; this will produce a relative error of zero at these
      ; places even if the value (denom) is 0.0:
      zerrs = where(errors EQ 0.0, nzerrs)
      if nzerrs GE 1 then ydenom(zerrs) = 1.0
      ; now, any ydenom which is (still) zero is matched
      ; with a non-zero error - hence the relative error
      ; is huge - call it 10^3 - and set the denom to error/10^3:
      zeros = where(ydenom EQ 0.0, nzs)      
      if nzs GE 1 then begin
        ydenom(zeros) = errors(zeros)/1000.0
      end
    end
    ; convert to relative errors by dividing by (the non-zero!)
    ; tabulated Y values...
    errors = errors/ydenom
  end
  max_errors(it) = MAX(errors)
end

RETURN, max_errors
END