pro read_positions, pimage, flag
@common_blocks.inc

;; Initialize arrays and counters
pimage = fltarr(xy_sl_size, xy_sl_size)

nparticles = TOTAL(grid_info.num_particle)
if nparticles le 0 then return


nGrids     = n_elements(grid_info)
inGrids    = intarr(nGrids)
count      = 0
nParticles = 0L
slab       = dblarr(6)

;; Setup slab
slice_width = slice_size[2] - slice_size[0]
other_sub = where(slice_ori ne 0)
other_sub = other_sub[0]
slice_sub = where(indgen(3) ne other_sub)
slab_coord = [center[other_sub] - slice_width/2, $
              center[other_sub] + slice_width/2]

dim_one = -1
count = 0
for dim = 0, 2 do begin
    if (other_sub eq dim) then begin
        slab[2*dim] = slab_coord[0]
        slab[2*dim+1] = slab_coord[1]
    endif else begin
        slab[2*dim] = slice_size[count]
        slab[2*dim+1] = slice_size[count+2]
        dim_one = (dim_one eq -1) ? dim : dim_one
        dim_two = dim
        count = count + 1
    endelse
endfor
count = 0


index_range    = grid_info.End_index[0]-grid_info.Start_index[0]+1
delta_distance = (DOUBLE(grid_info.Right_edge[0] - $
                         grid_info.Left_edge[0])/DOUBLE(index_range))
del_pix =   DOUBLE(xy_sl_size)*delta_distance/slice_width
  
grid_flag  = grid_in_projection(grid_info, DOUBLE(slice_ori), DOUBLE(slice_size)) 
inGrids  = where((del_pix GE 0.1) AND (del_pix LE xy_sl_size/3) AND (grid_flag gt 0) $
                 AND (grid_info.num_Particle GE 1) $
                AND (grid_info.left_edge[other_sub] le slab[2*other_sub+1]) $
                AND (grid_info.right_edge[other_sub] ge slab[2*other_sub])   ) 

inNum = N_ELEMENTS(inGrids)
nparticles = TOTAL(grid_info[inGrids].num_particle)
if verbose then print, inNum, ' Grids have ', nparticles,' relevant particles. Reading those.'

;return

;; Allocate memory
partPos = dblarr(nParticles,3)
colors  = fltarr(nParticles)
arrayPos = 0L

;; Read particle positions from files
for i = 0, inNum-1 do begin
    exi = file_test(data_dir+grid_info[inGrids[i]].particle_file, /READ)
    if exi then begin
        for dim = 0, 2 do begin
            read_hdf, data_dir+grid_info[inGrids[i]].particle_file, $
              grid_info[inGrids[i]].num_baryon_fields+dim+3, trash
            partPos[arrayPos:(arrayPos+grid_info[inGrids[i]].num_particle-1), dim] = trash
        endfor  ;; end: dims
;        read_hdf, data_dir+grid_info[inGrids[i]].particle_file, $
;          grid_info[inGrids[i]].num_baryon_fields+6, trash
;        colors[arrayPos:(arrayPos+grid_info[inGrids[i]].num_particle-1)] = trash
        arrayPos = arrayPos + grid_info[inGrids[i]].num_particle
    endif
endfor  ;; end: grids
if verbose then print, 'read ',nParticles, ' from ', inNum, ' grids.'
;; Throw away particles outside of slab
partInside = where(partPos[*,0] gt slab[0] and partPos[*,0] lt slab[1] and $
                   partPos[*,1] gt slab[2] and partPos[*,1] lt slab[3] and $
                   partPos[*,2] gt slab[4] and partPos[*,2] lt slab[5], cnt)
if (cnt gt 0) then begin
    partPosTemp = partPos[partInside,*]
    partPos = 0.
    colors = temporary(colors[partInside])
endif else return

if verbose then print, 'Using ', cnt,' of them.'
devicePart1 = (partPosTemp[*,dim_one]-slice_size[0])*xy_sl_size/slice_width
devicePart2 = (partPosTemp[*,dim_two]-slice_size[1])*xy_sl_size/slice_width

colors = alog10(colors > 1e-10)
if n_elements(flag) gt 1 then begin ; overlay scale to image data
iplot, partpostemp[*,0], partpostemp[*,1], partpostemp[*,2], LINESTYLE=6, sym_index=3
    mic = min(colors, max=mac)
    if mic eq mac then mac = flag[0]
    pimage[devicePart1, devicePart2] = (colors-mic)/(mac-mic)*(flag[1]-flag[0])+flag[0]
end else  pimage[devicePart1, devicePart2] = colors
;; Display particles
;maxLevel = max(colors)
;for i = 0L, cnt-1 do $
;  plots, devicePart1[i], devicePart2[i], /device, psym=3, $
;  color=round((maxLevel-colors[i])*255./maxLevel)
plot, [slice_size[0], slice_size[2]], [slice_size[1], slice_size[3]], $
  /nodata, ystyle=1, xstyle=1
return
end