pro construct_3D_data, grid_info, cube_center, cube_length,$
             data_dir,data_field,cube_size, cube_data, interpolate=interp
; using grid_info read in data (found in data_dir directory) of all
; necessary grids for specified cube ;
; cube_length: side length [0..1]
; cube_center: [x,y,z]
  if N_elements(interp) eq 0 then interp = 0
  num_dim = 3
; rough check of input data
  cube_data = fltarr(cube_size, cube_size, cube_size)
  cube_data[*,*,*] = 0.
  
  IF (N_ELEMENTS(grid_info) lt 1) THEN BEGIN
      print, 'construct_3D_data: no grid information in grid_info:',$
        grid_info
      STOP
  ENDIF
  dim_grid_info = size(grid_info)
  num_of_grids  = dim_grid_info(1)

; determine grids that have sufficient res and are at least partially
; in the cube 
; rough first pass
  u_g = grid_in_cube(grid_info, cube_center, cube_length, cube_size)
  
  index_range = u_g.End_index-u_g.Start_index+1
  delta_distance = DOUBLE(u_g.Right_edge - u_g.Left_edge)/FLOAT(index_range)

;adjust coordinates according to finest level that will be included
  highest_level = (where(u_g.level eq max(u_g.level)))[0] ; do this in units of
  fdx = delta_distance[0,highest_level]
  cube_length = 2D^(ROUND(alog(cube_length)/alog(2)))
  print, 'cube length:', cube_length
;  cube_length = cube_length-(cube_length mod fdx[0]/cube_size)
  print, 'old center:', cube_center
  cube_center = cube_center-(cube_center mod fdx[0])
  print, 'new center:', cube_center

  u_g = grid_in_cube(grid_info, cube_center, cube_length, cube_size)
  index_range = u_g.End_index-u_g.Start_index+1
  delta_distance = DOUBLE(u_g.Right_edge - u_g.Left_edge)/FLOAT(index_range)
  num_of_grids  = n_elements(u_g)

; min left edge and scaling
  min_left  = [cube_center[0],cube_center[1],cube_center[2]]-cube_length/2.D > 0.D
  max_right = [cube_center[0],cube_center[1],cube_center[2]]+cube_length/2.D < 1.D
  ndx = min(max_right-min_left)/DOUBLE(cube_size)

  vertexc = 0
  if vertexc eq 0 then begin 
     a = get_data(u_g, data_field, left_index=start, right_index=righti, $
                  left_edge = min_left, right_edge=max_right,interpolate=interp)

     this_level = 0
     if interp eq 1 then cor = 1 else cor = 0
     if interp then print, 'interpolating 3d cube' else print, 'neighrest neighbor sampling'
     FOR i=0L,num_of_grids-1 DO BEGIN
; check this
        dx = delta_distance[*,i]
;        breaki
        if max(dx/ndx) lt 100. then begin  
           temp = (*a[i].data[0])
           sub_cube_size = dx/ndx*(a[i].dims-cor)
           case interp of
              0: BEGIN
                 sub_cube= rebin(temp, sub_cube_size[0],sub_cube_size[1],sub_cube_size[2], /sample)
              END
              1: BEGIN
                 x = (findgen(sub_cube_size[0])+0.5)/sub_cube_size[0] * a[i].dims[0]
                 y = (findgen(sub_cube_size[1])+0.5)/sub_cube_size[1] * a[i].dims[1]
                 z = (findgen(sub_cube_size[2])+0.5)/sub_cube_size[2] * a[i].dims[2]

                 sub_cube= interpolate(temp, x, y, z, /grid)
              END
           endcase


           led = (u_g[i].Left_edge + dx*start[i,*])          ; left edge of data available (this grid)
           red = (led + dx*(a[i].dims-cor))          ; left edge of data available (this grid)
;           if total((red-min_left) le 0) ne 0 then break
           lin = FLOOR((led - min_left)/ndx)                        ; index in cube
           rin = FLOOR((red - min_left)/ndx)-1

           helpi = [0,0,0]
           din = rin- (lin > 0) 
           for j=0,2 do begin
              if lin[j] lt 0 then begin
                 helpi[j] = abs(lin[j]) < (sub_cube_size[j]-1)
                 lin[j] = lin[j] > 0
              endif
              if rin[j] gt cube_size-1 then begin
                  din[j] = cube_size-1-lin[j]
              endif
           endfor
           rin = lin+din
           rsi = helpi+din
           
           cube_data[lin[0]:rin[0],      $
                     lin[1]:rin[1],      $
                     lin[2]:rin[2]       ] = sub_cube[ helpi[0]:rsi[0], $
                                                       helpi[1]:rsi[1], $
                                                       helpi[2]:rsi[2]]
           
           next_level = u_g[(i+1)<num_of_grids-1].level 
        endif
     ENDFOR
     
  endif else begin              ; vertex centered data and 
     if not data_is_vertexc then begin 
;        start  = start  - 1   > 0         ; 
;        righti = righti + 1   <
        

     endif
  endelse

  dmin = min(cube_data, max=dmax) ;  add a little bit to the minimum value to avoid zeros
  print, 'min and max of data:', dmin, dmax
  if dmin eq 0. then cube_data = temporary(cube_data) > (dmin + 1e-30*(dmax-dmin)) 
;breaki
 print, 'done extracting cube'

END