FUNCTION grid_in_projection, grid_info, slice_ori, slice_coord, depth
; determine whether the specified grid  is at least partially in the
; slice  ; defined by the 3 component vector slice_ori 
; return for each grid 1 if it is and 0 if it is not in slice
  if N_ELEMENTS(where(grid_info[0].dim[*] gt 0)) eq 2 then begin
     print, 'grid_in_projection:  2 D data: all grids in slice'
     return, lindgen(N_elements(grid_info))
  endif

  axes = FIX(slice_ori gt 1.e-30)
  IF (Fix(TOTAL(axes)) gt 1) then BEGIN
      print, 'grid_in_projection: slice_ori:', slice_ori,' not allowed !'
      RETURN, -1
  ENDIF
  const_sub =  where(axes eq max(axes))
  const_sub = const_sub(0)
;  print, 'const_sub', const_sub
  other_subs = where(axes eq 0)
  Left_edge  = grid_info.Left_edge(const_sub,*)
  o_Left_edge  = grid_info.Left_edge(other_subs,*)
  Right_edge = grid_info.Right_edge(const_sub,*)
  o_Right_edge  = grid_info.Right_edge(other_subs,*)

  help = ROUND(Left_Edge) < 0
  help = ((slice_coord(0) lt o_Right_edge(0,*))     and $
          (slice_coord(1) lt o_Right_edge(1,*))     and $
          (slice_coord(2) gt o_Left_edge(0,*))      and $
          (slice_coord(3) gt o_Left_edge(1,*))      and $
          (depth[0]       le Right_edge[0, *])          and $
          (depth[1]       ge Left_edge[0, *]) )

;  print, 'there are ', total(grid_in_projection),' grids at least partially in projection'
  return, help
END

function grid_in_bounds, grid_info, bounds ; find intersections with bounds[0:2] : bounds[3:5]
    return, TOTAL([grid_info.right_edge, -grid_info.left_edge] + $
                  CMREPLICATE(bounds*[-1,-1,-1,1,1,1], N_ELEMENTS(grid_info)) gt 0, 1, /INT) eq 6
end