function set_parent_start_end_indices, grid_info
  print, 'start set_parent_start_end_indices'
  b =  0
  count = histogram(grid_info.level,reverse_indices=r)
  total_cell_count =  0L
  num_of_levels = N_elements(count) 

  for i=1, num_of_levels-1 DO BEGIN
     if count[i]-1 eq 0 or (r[i]-r[i-1]) le 0 then continue
     b = r[r[i-1]:(r[i]-1)]
     lcells = 0L
     FOR j = 0L,N_elements(b)-1 DO BEGIN
        cg = b[j]
; set parent start and end indeces 
        IF (grid_info[cg].parent GT 0) THEN BEGIN
           pi = where(grid_info[b].num eq grid_info[cg].parent, tcount)
           if tcount gt 0 then begin
              pi  = b[pi]
              p = grid_info[pi]
              c = grid_info[cg]
              pd = p.dim-1
              pdx = (p.right_edge-p.left_edge)/DOUBLE(pd+1L)
              grid_info[cg].parent_s_index = ROUND((c.left_edge-p.left_edge)/pdx)  < pd
              grid_info[cg].parent_e_index = (ROUND((c.right_edge-p.left_edge)/pdx)-1) < pd
           endif else print, 'grid ', grid_info[cg].parent, ' was not found'
        ENDIF
     ENDFOR
  ENDFOR

  print, 'end set_parent_start_end_indices'
return, grid_info
end

function find_and_set_parent_grids, grid_info
  print, 'start find_and_set_parent_grids'
if N_elements(grid_info) eq 1 then return, grid_info
nlb = histogram(grid_info.level, reverse_indices=r)
nl = N_elements(nlb)

if nl eq 1 then return, grid_info ; only one level. Nothing to do.
if max(grid_info.parent_e_index) gt 0 then return, grid_info ; parents already set

for i=1L,nl-1 do begin          ; loop over levels (not lowest level grid ones which have no parents)
   gwnp = 0
   if (nlb[i-1] eq 0) or (nlb[i] eq 0) then begin
      print, 'no grids to consider. Strange.'
      continue
   endif
   potp = grid_info[r[r[i-1]:r[i]-1]] ; potential parent grids
   parenttimes = potp[UNIQ(sort(potp.time))].time
   ind = where(parenttimes le grid_info[r[r[i]]].time, count)
   if count gt 0 then plevelt = max(parenttimes[ind]) else plevelt=-1   
   indi = where(potp.time eq plevelt)
   if indi[0] ge 0 then pp = potp[indi] else break
   tl = grid_info[r[r[i]:r[i+1]-1]]
   for j=0L,nlb[i]-1 do begin       ; loop over grids on this level
;      cg = grid_info[r[r[i]+j]]     ; current grid
      cg = tl[j]     ; current grid
      cond =  (cg.left_edge[0]  ge pp.left_edge[0]) AND    $
        (cg.left_edge[1]  ge pp.left_edge[1]) AND    $
        (cg.left_edge[2]  ge pp.left_edge[2]) AND    $
        (cg.right_edge[0] le pp.right_edge[0] ) AND  $
        (cg.right_edge[1] le pp.right_edge[1] ) AND  $
        (cg.right_edge[2] le pp.right_edge[2] )
      ind = where(cond eq 1, count)
;      ind = ind(where(ind ge 0))
 ;    if count gt 1 then begin
 ;        print, cg.num,': grid with ', strcompress(N_elements(ind)),$
 ;          ' parents. Using first one ...'
 ;        ind = ind[0]
 ;     endif
      if count eq 1 then begin  ; this is the one
          p = pp[ind]
          ci = r[r[i]+j]
          grid_info[ci].parent = p.num 
          pd = p.dim-1
          pdx = (p.right_edge-p.left_edge)/DOUBLE(pd+1L)
          grid_info[ci].parent_s_index = ROUND((cg.left_edge-p.left_edge)/pdx)  < pd
          grid_info[ci].parent_e_index = (ROUND((cg.right_edge-p.left_edge)/pdx)-1) < pd
      endif
  endfor                        ; end loop over grids on this level

  print, i, ' level found parents'
endfor
  print, 'end find_and_set_parent_grids'
return, grid_info
end


function all_perms_le, a1, b1, a2, b2, a3, b3
return, ((reform(a1)#(1./b1)) le 1) + ((reform(a2) # (1./b2)) le 1) + ((reform(a3) # (1./b3)) le 1) 
end

function experimental_fast_find_and_set_parent_grids, grid_info
  print, 'start find_and_set_parent_grids'
if N_elements(grid_info) eq 1 then return, grid_info
nlb = histogram(grid_info.level, reverse_indices=r)
nl = N_elements(nlb)

if nl eq 1 then return, grid_info ; only one level. Nothing to do.
if max(grid_info.parent_e_index) gt 0 then return, grid_info ; parents already set

for i=1L,nl-1 do begin          ; loop over levels (not lowest level grid ones which have no parents)
   gwnp = 0
   if (nlb[i-1] eq 0) or (nlb[i] eq 0) then begin
      print, 'no grids to consider. Strange.'
      continue
   endif
   pp = grid_info[r[r[i-1]:r[i]-1]] ; potential parent grids with the right level
   npotp = N_elements(pp)                            ; this is how many potential parents we have
   bi  = r[r[i]:r[i+1]-1]
;   fi  = LINDGEN(N_elements(bi))
   tl = grid_info[bi]

   tll = tl.left_edge
   tlr = tl.right_edge

   cond =      all_perms_le( pp.left_edge[0] , tll[0,*],    $
                             pp.left_edge[1] , tll[1,*],    $
                             pp.left_edge[2] , tll[2,*])
   cond += $
     transpose(all_perms_le(tlr[0,*],  pp.right_edge[0], $
                            tlr[1,*],  pp.right_edge[1], $
                            tlr[2,*],  pp.right_edge[2])) eq 6
   mi = array_indices([npotp, nlb[i]], lindgen(npotp*nlb[i]),/DIMENSIONS)
   hi = histogram(cond , reverse_indices=s)
   count = hi[1]
   ind  = s[s[1]:s[2]-1]
   ci = bi[mi[1,ind]]
   pi = mi[0,ind]
   grid_info[ci].parent = pp[pi].num 
;if i eq 1 then breaki
  print, 'level',strcompress(string(i))+' found parents'
endfor
  print, 'end find_and_set_parent_grids'
return, grid_info
end


function fast_find_and_set_parent_grids, grid_info
  print, 'start find_and_set_parent_grids'
if N_elements(grid_info) eq 1 then return, grid_info
nlb = histogram(grid_info.level, reverse_indices=r)
nl = N_elements(nlb)

if nl eq 1 then return, grid_info ; only one level. Nothing to do.
if max(grid_info.parent_e_index) gt 0 then return, grid_info ; parents already set

for i=1L,nl-1 do begin          ; loop over levels (not lowest level grid ones which have no parents)
   gwnp = 0
   if (nlb[i-1] eq 0) or (nlb[i] eq 0) then begin
      print, 'no grids to consider. Strange.'
      continue
   endif
   pp = grid_info[r[r[i-1]:r[i]-1]] ; potential parent grids with the right level
   npotp = N_elements(pp)                            ; this is how many potential parents we have
   bi  = r[r[i]:r[i+1]-1]
;   fi  = LINDGEN(N_elements(bi))
   tl = grid_info[bi]

   tll = tl.left_edge
   tlr = tl.right_edge

   fnp = 0L                     ; found this many parents
   for j=0L,npotp-1 do begin       ; loop over all potential parents
       cp = pp[j]                    ; current potential parent
       cond =  (tll[0,*]  ge cp.left_edge[0]) AND    $
               (tll[1,*]  ge cp.left_edge[1]) AND    $
               (tll[2,*]  ge cp.left_edge[2]) AND    $
               (tlr[0,*]  le cp.right_edge[0] ) AND  $
               (tlr[1,*]  le cp.right_edge[1] ) AND  $
               (tlr[2,*]  le cp.right_edge[2] )
;      ind = where(cond eq 1, count, complement=comp)
      hi = histogram(cond , reverse_indices=s)
      count = hi[1]
      if count ge 1 then begin  ; these grids have cp = pp[j]  as parent
          ind  = s[s[1]:s[2]-1]
          fnp += count
          ci = bi[ind]
          grid_info[bi[s[s[1]:s[2]-1]]].parent = cp.num 
          pd = cp.dim
          pdx = (cp.right_edge-cp.left_edge)/DOUBLE(pd)
          for k=0,2 do begin
              grid_info[ci].parent_s_index[k] = ROUND((grid_info[ci].left_edge[k] - cp.left_edge[k])/pdx[k])  < (pd[k]-1)
              grid_info[ci].parent_e_index[k] = ROUND((grid_info[ci].right_edge[k]- cp.left_edge[k])/pdx[k])-1 < (pd[k]-1)
          endfor
      endif

;      if comp[0] eq -1 then break  ; no grids left that are looking for a parent
  endfor                        ; end loop over grids on this level

  print, 'level',strcompress(string(i))+' found parents'
endfor
  print, 'end find_and_set_parent_grids'
return, grid_info
end