pro pack_particles, field_names, a
; for gadget or partiview data repackage particles
  readall        = {readall,                         $
                    file_name: "",                   $
                    Np_total: 0L,                    $
                    dims: [0,0,0],                   $
                    li: [0,0,0],                     $  ; left index
                    ri: [0,0,0],                     $  ; right index
                    N_dfields: 0B,                   $
                    DATA_FIELDS: STRARR(100),        $
                    Data:     PTRARR(100)   }

  cnt = 0
  FOR i=0,N_elements(fields)-1 DO BEGIN
     field = fields[i]
     all[cnt].file_name   = gident
     all[cnt].np_total    = N_elements(data)
     s = size(data)
     if (s[0] eq 2) then s[3] = 1
     all[cnt].li[*]       = this_start[*]
     all[cnt].ri[*]       = (s[1:3]+this_start-1)[*]
     all[cnt].dims        = all[cnt].ri - all[cnt].li + 1
     all[cnt].N_Dfields   = N_df_requested
     all[cnt].data_fields[i] = fields[i]
     if fields[i] eq 'particle massdensity' then data *= product(dd)
     if fields[i] eq 'mass' then data *= product(dd)
     all[cnt].Data[i]     = ptr_new(parti.d[i], /no_copy)
  ENDFOR                        ; loop over fields
  cnt += 1

end

function return_null_data, VariableNames
  readall        = {readall,                         $
                    file_name: "",                   $
                    Np_total: 1L,                    $
                    dims: [1,1,1],                   $
                    li: [0,0,0],                     $  ; left index
                    ri: [1,1,1],                     $  ; right index
                    N_dfields: 0B,                   $
                    DATA_FIELDS: STRARR(100),        $
                    Data:     PTRARR(100)   }
  result = REPLICATE(readall, 1)
  result[0].data_fields = VariableNames
  fake_value = [-1.0]
  for field = 0, N_ELEMENTS(VariableNames) do $
     result[0].data[field] = PTR_NEW(fake_value)
  return, result
end

function voronoi_volume, x,y,z, fractional=frac
    qhull, x, y, z, tri, /delaunay
    print, 'qhull done'
    ntri = n_elements(tri[0,*])
    vol_tri = fltarr(ntri)
    for itri = 0L, ntri-1 do $
      vol_tri[itri] = abs(determ([[x[tri[*,itri]]], $
                                      [y[tri[*,itri]]], $
                                      [z[tri[*,itri]]], $
                                      [replicate(1,4)]])) / 6.
    print, 'computed vertices'
    h = histogram(tri, reverse_indices=ri)
    vol = x
    vol[*] = 0.
    rr = ri/4
    for i=0L,N_elements(x)-1 DO vol[i]=total(vol_tri[rr[ri[i]:(ri[i+1]-1)]])/4.
    print, 'computed volume'
    if N_elements(frac) gt 0 then begin
       volume = total(vol)
       vol /= volume            ; normalize so we return fractional volume
    endif
return, vol
end

function get_particle_data,  usegi, VariableNames, DERIVED_FIELDS=derived_fields,$
                   STRIDE=stride, redshift=zred, center=tcenter
@common_blocks.inc
if N_elements(tcenter) le 1 then tcenter = center
  Nderived = 20
  id = 0L
  dfields = PTRARR(Nderived) 
  dfields[id++] = PTR_NEW(["abs(part. velocity)", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
  dfields[id++] = PTR_NEW(["particle radial velocity","particle_position_x", "particle_position_y", "particle_position_z", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
  dfields[id++] = PTR_NEW(["particle tangential velocity", "particle_position_x", "particle_position_y", "particle_position_z", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
  dfields[id++] = PTR_NEW(["particle radius", "particle_position_x", "particle_position_y", "particle_position_z"])
  dfields[id++] = PTR_NEW(["particle cyl. radius", "particle_position_x", "particle_position_y", "particle_position_z"])
  dfields[id++] = PTR_NEW(["particle age", "creation_time"])
  dfields[id++] = PTR_NEW(["enclosed mass in particles","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
  dfields[id++] = PTR_NEW(["particle voronoi density","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
  dfields[id++] = PTR_NEW(["particle voronoi radius","particle_position_x", "particle_position_y", "particle_position_z"])
  dfields[id++] = PTR_NEW(["particle voronoi volume","particle_position_x", "particle_position_y", "particle_position_z"])
; the voronoi calculations are done in read_particle after all
; particles have been read.
  dfields = dfields[0:id-1]

  if N_elements(stride) eq 1 then stride = LONARR(N_elements(usegi))+stride
  if  KEYWORD_SET(derived_fields) then return, dfields $ ; return field_names & stop
  else begin
     Nvar = N_ELEMENTS(VariableNames)
     data_format = determine_data_format_from_file_name(usegi[0].baryon_file)
     if data_format lt 0 then return, 0
; determine all data fields to read and read data
     field_names = identify_unique_field_names(VariableNames, dfields)

     if data_format eq 4 then read_all_grid_hdf,  $
        usegi, field_names, a, /VERBOSE, STRIDE=stride 
     if ((data_format eq 5) OR (data_format eq 6)) then $
        read_all_grid_hdf5, usegi,  field_names, a, $
                            /VERBOSE, STRIDE=stride
     if ((data_format eq 1) or (data_format eq 2)) then $
        pack_particles, a
     
     ax = double(slice_ori eq 0.)
     print, field_names
     for i=0, Nvar-1 do begin   ; if derived field was requested compute it here
        match = is_derived_field_number(VariableNames[i], dfields)
        if match ge 0 then begin
           for j=0L,N_elements(a)-1 DO BEGIN ; loop over grids
               if ptr_valid(a[j].data[0]) then begin
                   case match of 
                       -1: junk = ''
                       0: a[j].data[i] = ptr_new(sqrt(fn(a[j], "particle_velocity_x")^2+ $
                                                      fn(a[j], "particle_velocity_y")^2+ $
                                                      fn(a[j], "particle_velocity_z")^2))
                       1: a[j].data[i] = ptr_new(radial_velocity(tcenter, $
                           fn(a[j], "particle_position_x"),fn(a[j], "particle_position_y"),fn(a[j], "particle_position_z"), $
                           fn(a[j], "particle_velocity_x"),fn(a[j], "particle_velocity_y"),fn(a[j], "particle_velocity_z"),/cons_vel))

                       2: a[j].data[i] = ptr_new(tangential_velocity(tcenter, $
                           fn(a[j], "particle_position_x"),fn(a[j], "particle_position_y"),fn(a[j], "particle_position_z"), $
                           fn(a[j], "particle_velocity_x"),fn(a[j], "particle_velocity_y"),fn(a[j], "particle_velocity_z")))

                       3: a[j].data[i] = ptr_new(sqrt((tcenter[0]-fn(a[j], "particle_position_x"))^2 + $
                                                      (tcenter[1]-fn(a[j], "particle_position_y"))^2 + $
                                                      (tcenter[2]-fn(a[j], "particle_position_z"))^2))

                       4: a[j].data[i] = ptr_new(sqrt(ax[0]*(tcenter[0]-fn(a[j], "particle_position_x"))^2 + $
                                                      ax[1]*(tcenter[1]-fn(a[j], "particle_position_y"))^2 + $
                                                      ax[2]*(tcenter[2]-fn(a[j], "particle_position_z"))^2))

                       5: a[j].data[i] = ptr_new(times[time_index] - fn(a[j], "creation_time"))
                       6: a[j].data[i] = $  ; allocate array
                          ptr_new(fn(a[j], "particle massdensity")) ; 
                       7: a[j].data[i] = $
                          ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
                       8: a[j].data[i] = $
                          ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
                       9: a[j].data[i] = $
                          ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
                   endcase
;                   print, max(*a[j].data[i])
;                   a[j].data_fields[i] = (*dfields[match])[0]
               endif
           endfor               ; loop over grids
        endif                   ; is derived variable
     endfor                     ; loop over variables
  endelse                       ; if derived_field_names
return, a
end

pro store_particle_data, a, names
@common_blocks.inc
  for j=0, N_elements(names)-1 do begin
     cni = where((*parti.names) eq names[j])
     cname = names[cni]
     if total((a.np_total)) eq 0 then return
     temp = fltarr(total(a.np_total))
     if parti.read[cni] ne 1 then begin
        runi=0L
        for i=0L, N_elements(a)-1 do begin
            if ptr_valid(a[i].data[j]) then begin
                temp[runi:(runi+a[i].np_total-1)] = (*a[i].data[j])[*]
                runi += a[i].np_total
            endif
        endfor
        parti.d[cni] = ptr_new(temp, /no_copy)
     end
  endfor
  return
end

function fp, name, ind
@common_blocks.inc
  if ind[0] eq -1 then return, [0]
  cni = (where((*parti.names) eq name[0]))[0]
  if N_elements(ind) gt 0 then temp = (*parti.d[cni])[ind] $
  else temp = (*parti.d[cni])
return, temp
end

pro read_particle_data, fields=rfields
@common_blocks.inc

if N_elements(rfields) gt 0 then fields = rfields $
else fields = [0,1,2, parti.index]

; also read the particle type field if we are going to select on it
if parti.type gt 0 then begin
   typeind = where(*parti.names eq 'particle_type')
   if typeind[0] gt -1 then fields=[fields[*], typeind] else $
      print, 'read_particle_data: no field found with name "particle_type"'
endif

fields = remove_additional_occurance(fields)
nfields = [0]
for i=0, N_elements(fields)-1 do if parti.read[fields[i]] lt 1 then nfields = [nfields,fields[i]]
if N_elements(nfields) gt 1 then begin
   nfields = nfields[1:*]
   print, 'fields not read yet:', nfields
   ind= where(grid_info.num_particle gt 0)
   if (ind[0] ge 0) then begin
      a = get_particle_data(grid_info[ind], (*parti.names)[nfields],  $
                            STRIDE=parti.every)
   endif else begin
      a = return_null_data((*parti.names)[nfields])
   endelse
   store_particle_data, a, (*parti.names)[nfields]
   parti.read[fields] = 1
endif

ind = where(strmatch((*parti.names)[nfields] , '*voronoi*') eq 1, count)
if count gt 0 then begin
   ind = nfields[ind[0]]
print, 'set voronoi fields #:', ind
   if strmatch((*parti.names)[ind], '*density') then $ ; voronoi density
   *parti.d[ind] = ((*parti.d[ind])/voronoi_volume(*parti.d[0], *parti.d[1], *parti.d[2])) else $
      if strmatch((*parti.names)[ind], '*radius') then $ ; voronoi radius
         *parti.d[ind] = (voronoi_volume(*parti.d[0],*parti.d[1],*parti.d[2])^(1./3.)) else $ 
      if strmatch((*parti.names)[ind], '*volume') then $ ; voronoi radius
         *parti.d[ind] = (voronoi_volume(*parti.d[0],*parti.d[1],*parti.d[2])) $ 
      else print, 'read_particle_data: did not find the definition ...'
print, (*parti.names)[ind]
endif

ind = where(strmatch((*parti.names)[nfields] , 'enclosed mass in part*') eq 1, count)
if count gt 0 then begin
   ind = nfields[ind[0]]
   r = sqrt((*parti.d[0])^2+(*parti.d[1])^2+(*parti.d[2])^2)
   is = sort(r)
   m = TOTAL(((*parti.d[ind])[is]), /cumulative)
   n = n_elements(r)
   ii = (lindgen(n))[is]
   
   *parti.d[ind] = m[ii[is]]

endif

end

function render_particle_data, HISTOGRAM=histogram
@common_blocks.inc
if N_ELEMENTS(*parti.names) eq 0 or (*parti.names)[0] eq '' then return, -1
names = (*parti.names)

axes = (slice_ori gt 1.e-30)
const_sub =  where(axes eq max(axes))
const_sub = const_sub(0)
other_subs = where(axes eq 0)
dim_one = other_subs[0]
dim_two = other_subs[1]
slice_width = slice_size[2] - slice_size[0]

;ind = where(*parti.names eq 'particle mass', count)
;if count gt 0 then read_particle_data, fields = [0, 1, 2, ind[0], parti.index]
;else 
read_particle_data

depthp = (slice_value+depth_value/2.) 
depthm = (slice_value-depth_value/2.) 
;if depthm lt 0. then depthp = depth_value
;if depthp gt 1. then depthm = 1.-depth_value

;select region
nind= where((*parti.d[dim_one] gt slice_size[0]) and $
            (*parti.d[dim_one] lt slice_size[2]) and $
            (*parti.d[dim_two] gt slice_size[1]) and $
            (*parti.d[dim_two] lt slice_size[3]) and $
            (*parti.d[const_sub] le depthp) and $
            (*parti.d[const_sub] ge depthm) )

sin = sort(fp(names[parti.index],nind))
ind = nind[sin]

typeind = where(*parti.names eq 'particle_type')
if parti.type gt 0 then begin
   th = histogram((*parti.d[typeind[0]])[ind])
   print, STRCOMPRESS(FIX(total(th gt 0))), ' different types of particles found.'
    if typeind[0] gt -1 then begin
      nind = where((*parti.d[typeind[0]])[ind] eq parti.type)
      if nind[0] gt -1 then ind = ind[nind] else print, 'no matching particle types found: plotting all instead'
   end else print, 'no field with name "particle_type" read'
endif
npoints = N_elements(ind)
print, 'including, ', STRCOMPRESS(npoints), ' points'

devicePart1 = FIX((fp([names[dim_one]],ind)-slice_size[0])*(xy_sl_size-1)/slice_width)+1
devicePart2 = FIX((fp([names[dim_two]],ind)-slice_size[1])*(xy_sl_size-1)/slice_width)+1

if npoints lt 10 then print, devicePart1, devicePart2

;pointimage = hist2d(devicePart1, devicePart2, $
;                      min1=1.,min2=1.,max1=xy_sl_size,max2=xy_sl_size, $
;                    binsize1=parti.size,binsize2=parti.size, density=rho)
;pointimage = rho

himage = hist2d(devicePart1, devicePart2, fp(names[parti.index],ind), binsize1=parti.size,binsize2=parti.size,$
                min1=1.,min2=1.,max1=xy_sl_size,max2=xy_sl_size, density=rho) ;/pointimage

print, size(himage)
ind  = where(rho gt 0, count)
if not shouldbe_histogram(names[parti.index]) then begin
   if count gt 0 then himage[ind] /= float(rho[ind])
print, min(rho), max(rho)
endif

secmin = min(himage)
ind  = where(himage gt secmin, count)
;print, secmin


if shouldbe_logged(names[parti.index]) or parti.index le 2 then begin
;   himage = himage > 0.5*secmin
   himage = alog10(himage > secmin)
endif

s = size(himage)
dix = s[1]
diy = s[2]
locs = lindgen(dix*diy)
nind = ARRAY_INDICES([dix,diy], locs, /DIMENSIONS)

  dx_stuetz = fltarr(dix*diy)
  dx_stuetz[*]= float(slice_size[2]-slice_size[0])/(dix-1)
  x_stuetz = ((0.5+nind[0,*])*dx_stuetz+slice_size[0])[*]
  y_stuetz = ((0.5+nind[1,*])*dx_stuetz+slice_size[1])[*]
  z_stuetz = himage[*]
  n_stuetz = N_elements(z_stuetz)
if N_elements(z_stuetz) ne dix*diy then breaki
;  if parti.size lt 4 then $
;     pimage = (rebin(himage, dix*parti.size,diy*parti.size,  sample=1-(interpolate_i < 1)))[0:xy_sl_size-1, 0:xy_sl_size-1] $
;  else  
  pimage = render_image()

return, pimage
end