FUNCTION interp_pruned_x, x, ys, ac, ABS_ERROR = abs_error, $
	PLOT_PRUNE=plot_prune
;
; This procedure is used to "prune" a (possibly multi-Y-column)
; lookup table by returning the indices of the xs array which,
; if they alone are extracted from the table and used with
; linear interpolation, will lead to interpolated values with
; errors less than some specified level (for each Y-column).
;
; EXAMPLE of USE
;
; hetgcal> meg = rdb_read(!DDHETGCAL+'/cip/hetgmegD1996-11-01greffN0002.rdb')
; hetgcal> acs = [0.001, 0.003, 0.01, 0.03, 0.1]
; hetgcal> for ia=0,4 do xs_needed = interp_pruned_x(meg.energy, $
;	[[meg.op1], [meg.oz]], acs(ia)*[1.,1.])
;
; Need  791 points out of  5001 to get accuracy of  0.100000 %
;
; Need  441 points out of  5001 to get accuracy of  0.300000 %
;
; Need  250 points out of  5001 to get accuracy of  1.00000 %
;
; Need  142 points out of  5001 to get accuracy of  3.00000 %
;
; Need  74 points out of  5001 to get accuracy of  10.0000 %
;
; 
; INPUTS
;
;   x - the "X" array of the lookup table
;
;  ys - One or more arrays that correspond to looked-up values
;       w.r.t. x, that is ys(*,0), ..., y(*,N-1) are each arrays
;       with the same lenght as x.  ys is conveniently formed as:
;        ys = [ [lookup1], [lookup2], ..., [lookupN-1] ]
;
;  as - are the desired accuracy for each of the y's in relative
;       terms, that is:
;
;        | (pruned_Y(m) - actual_Y(m) ) / actual_Y(m) | < as(m)
;
; RETURNED VALUE
;  The indices of the rows in the table that are needed to achieve the
;  desired accuracy are returned.

; KEYWORDS
;
;  ABS_ERROR - if this keyword is set then the as values are 
;              interpreted as error values in an absolute sense,
;              that is:
;
;        | (pruned_Y(m) - actual_Y(m) ) | < as(m)
;
; FUNCTION REQUIRED
;    interp_prune_err(ixl, ixn, x, ys) - calculates and returns the
; max error for each table when interpolated from indices ixl to ixn.
;
; ALGORITHM overview: ixl points at the index of that last point
; which we've included in the pruned x list and ixn points at
; the candidate next point.  The errors for the tables if 
; linear interpolation is used between ixl and ixn is calculated.
; If all of the tables have errors less than desired then
; the value of ixn is incremented by one and we repeat, otherwise
; (errors are too big in at least one table) we add the point
; ixn-1 to the pruned-x-list and set ixl = ixn-1 and repeat.
;
; --------------------------------------------------------------
; 1999-01-16 dd Initial version created.
; --------------------------------------------------------------
;

; Make sure these values are set so they can be passed...
if NOT(KEYWORD_SET(ABS_ERROR)) then abs_error = 0
if NOT(KEYWORD_SET(PLOT_PRUNE)) then plot_prune = 0

nxs = n_elements(x)

; Start by setting ixl to 0 and xs_needed to [0]:
xs_needed = [0]
ixl=0

; Here is the general loop to keep increasing ixn until the
; error is too large, then set ixl to that ixn-1 and continue...
; until ixl is equal to the last x value (always accepted).
while ixl LT nxs-1 do begin
  ; Show the progress...
  if KEYWORD_SET(PLOT_PRUNE) then print, $
	'  '+STRCOMPRESS(ixl)+'/'+STRCOMPRESS(nxs)

  ; catch the case that ixl is the next-to-last point, if so
  ; accept it... and next iteration we'll be done
  if ixl EQ nxs-2 then begin
    ; add the last point
    xs_needed = [xs_needed, (ixl+1)]
    ; and move ixl:
    ixl = ixl+1
  end else begin
    ; OK, now we can set ixn to ixl+1
    ixn = ixl + 1
    ; and start increasing ixn as long as the errors are OK
    err_OK = 1
    while err_OK do begin
      ixn = ixn + 1
      ; if ixn points outside the array then we're done -
      ; signal an error and ixl will be set to ixn-1
      if ixn GE nxs then begin
        err_OK = 0
      end else begin
        max_errors = interp_prune_err(ixl, ixn, x, ys, $
		PLOT_PRUNE=plot_prune, ABS_ERROR=abs_error)
        err_OK = ( MAX(max_errors GT ac) EQ 0 )
      end  
    end
    ; at this point ixn points at a point which gives too
    ; large an error, so set ixl = ixn-1 and continue...
    ixl = ixn-1
    xs_needed = [xs_needed, (ixl)]
  end

end ; of going through the x's.

if KEYWORD_SET(ABS_ERROR) then begin
  print, ' Need ' + STRCOMPRESS(n_elements(xs_needed)) + ' points out of ' + $
	STRCOMPRESS(nxs)+ ' to get accuracy of '+STRCOMPRESS(ac(0))
end else begin
  print, ' Need ' + STRCOMPRESS(n_elements(xs_needed)) + ' points out of ' + $
	STRCOMPRESS(nxs)+ ' to get accuracy of '+STRCOMPRESS(100.*ac(0))+' %'
end
print, ''

RETURN, xs_needed
END