;+
; SURF_CLIP: clip a 2D Z array, and two matching 1D arrays for x
;        and y.
;
; The function returns a 2D array that is truncated in x and y.
;
; XRANGE, YRANGE -> the range for x and y in data coordinates
; XCLIP, YCLIP -> return the truncated x and y arrays
;
; author: Martin Schultz, Harvard University (1998)
;
; EXAMPLE:   x=findgen(100)
;            y=findgen(50)
;            z=dist(100,50)
;            zc=surf_clip(z,x,y,xrange=[20,80],yrange=[20,40], $ ;          
;                     xclip=xc,yclip=yc)
;            surface,zc,xc,yc
;-
function surf_clip,z,x,y,xrange=xrange,yrange=yrange, $
              xclip=xclip,yclip=yclip

    ; clip z,x, and y arrays for surface plot
    ; x and y must be monotone

    if (n_elements(xrange) eq 0) then xrange=[min(x,max=xm),xm]
    if (n_elements(yrange) eq 0) then yrange=[min(y,max=ym),ym]

    ; find minimum and maximum index in x and y

    ix0 = min(where(x ge xrange[0]))  ; <<
    ix1 = max(where(x le xrange[1]))
    iy0 = min(where(y ge yrange[0]))
    iy1 = max(where(y le yrange[1]))  ; <<

    ; NOTE: you can restrain these indices to valid values as :
    ; ix0 = min(...) > 0    ; ix0 always at least 0
    ; ix1 = max(...) < (n_elements(x)-1)  ; prevent subscript range err

    if (ix0<ix1<iy0<iy1 ge 0) then begin
         xclip = x[ix0:ix1]   ; clipped x array
         yclip = y[iy0:iy1]   ; clipped y array
         return,z[ix0:ix1,iy0:iy1]   ; clipped z array
    endif else message,'Invalid range'

    return,-1

end