FUNCTION grid_in_slice, grid_info, slice_ori, slice_coord, inside=inside
; 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 inside is set then the slice value has to be within cell centers
; not just the edge
  if N_ELEMENTS(where(grid_info[0].dim[*] gt 0)) eq 2 then begin
     print, 'grid_in_slice:  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_slice: slice_ori:', slice_ori,' not allowed !'
      RETURN, -1
  ENDIF
  const_sub =  where(axes eq max(axes))
  const_sub = const_sub(0)

  index_range = grid_info.End_index-grid_info.Start_index+1
  delta_distance = (grid_info.Right_edge - grid_info.Left_edge)/DOUBLE(index_range)

  if N_elements(inside) ne 0 then half = 0.5D else half = 0.D
  other_subs = where(axes eq 0)
  Left_edge  = grid_info.Left_edge(const_sub,*)+half*delta_distance[const_sub,*]
  o_Left_edge  = grid_info.Left_edge(other_subs,*)+half*delta_distance[other_subs,*]
  Right_edge = grid_info.Right_edge(const_sub,*)-half*delta_distance[const_sub,*]
  o_Right_edge  = grid_info.Right_edge(other_subs,*)-half*delta_distance[other_subs,*]

  help = ROUND(Left_Edge) < 0
  help = ((slice_ori(const_sub) ge Left_edge)       and $
          (Right_edge gt slice_ori(const_sub))      and $
          (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,*))    )
  grid_in_slice = help
;  print, 'there are ', total(grid_in_slice),' grids at least partially in slice'
  return, grid_in_slice
END