import yt
import numpy as np

# Create fake particle data for a million particles randomly distributed
ppx, ppy, ppz, ppm = np.random.random([4, 1e6])
data = {'particle_position_x': ppx,
        'particle_position_y': ppy,
        'particle_position_z': ppz,
        'particle_mass': ppm}

ds = yt.load_particles(data, length_unit=3.08e24, mass_unit=1.9891e36)

ad = ds.all_data()
print ad['deposit', 'all_cic']

slc = yt.SlicePlot(ds, 2, ('deposit', 'all_cic'))
slc.set_figure_size(4)

#level 5 corresponds to a 32^3 array
cg = ds.covering_grid(level=5, left_edge=[0,0,0], dims=[32, 32, 32])

arr = cg['deposit', 'all_density']
print arr.max()
print arr.min()
print arr.mean()
print arr.shape

arr.write_hdf5('cic_deposited_data.h5')

import h5py

f = h5py.File('cic_deposited_data.h5')

print f['array_data'][:].min()
print f['array_data'][:].max()
print f['array_data'][:].mean()
print f['array_data'].shape