# Import a HartreePotential object.
potential = nlread('hartree_pot.hdf5', HartreePotential)[0]

# Calculate the average along z axis.
v_z = numpy.apply_over_axes(numpy.mean, potential[:,:,:], [0,1]).flatten()

# Add unit.
v_z = v_z * potential.unit()

# Get the volume element of the grid.
dX, dY, dZ = potential.volumeElement()
dz = dZ.norm()

# Print out the result.
shape = potential.shape()
for i in range(shape[2]):
    print(dz*i, v_z[i])