function combine_pointers, a, j
  N=N_elements(a)
  count = 0L
  for i=0L, N-1 Do count += a[i].Np_total
  Np = count
  f = dblarr(Np)
  count = 0L
  for i=0L, N-1 Do begin
     f[count:(count+a[i].Np_total-1)] = (*a[i].data[j])[*]
     count += a[i].Np_total
  endfor
return, f
end

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:'
      return
  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 = 1
;;   if interp then begin
;;         data_fields = ['x','y','z', data_field ]
;;         a = get_data(u_g, data_fields, left_index=start, right_index=righti, $
;;                      left_edge = min_left, right_edge=max_right)
;;         zero = -1e30
;;         zero_solution_under_subgrid, u_g, a, ZERO=zero
;;         f = combine_pointers(a, 3)
;;         ind = where(f ne zero)
;;         x = (combine_pointers(a,0))[ind]
;;         y = (combine_pointers(a,1))[ind]
;;         z = (combine_pointers(a,2))[ind]
;;         f = f[ind]
;;         a = 0.
;;         print,' interpolating on to grid ...'
;;         qhull, x,y,z, tet, /delaunay
;;         print, 'tetrahedra found. Now gridding ...'
;;         cube_data = qgrid3(x,y,z,f, tet, dimension=cube_size) > min(f)
;; ;        cube_data = grid3(x,y,z,f, ngrid=cube_size) > min(f)
;;         print, 'finished interpolating'
;;         goto, finish
;;      endif
;     endelse

  a = get_data(u_g, data_field, left_index=start, right_index=righti, $
               left_edge = min_left, right_edge=max_right)

  this_level = 0

        print, n_elements(uniq(1./delta_distance[0,*])), ' levels!'        
     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 = fix(dx/ndx)*(a[i].dims) ;< cube_size-1

           if interp then $
              sub_cube= congrid(temp, sub_cube_size[0],sub_cube_size[1],sub_cube_size[2], /center)$
           else $
              sub_cube= rebin(temp, sub_cube_size[0],sub_cube_size[1],sub_cube_size[2], sample=1)

           led = (u_g[i].Left_edge + dx*start[i,*])          ; left edge of data available (this grid)
           red = (led + dx*(a[i].dims))          ; 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
           lin = lin < rin
           rin = rin > lin
           rsi = helpi+din > 0
           
           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

;;         if  interp and ((u_g[(i+1)< (num_of_grids-1)].level ne  u_g[i].level) or (i eq num_of_grids-1))then begin
;;            smoothr = fix(2.*dx/ndx) < 50 > 2
;;            print, 'smooth radius:',smoothr
;; ;           cube_data = smooth(cube_data, smoothr , /edge_truncate)
;;        endif
     ENDFOR 
     



finish:
  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