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