from numpy import array
import pylab
from pymses.analysis.visualization import *
from pymses import RamsesOutput
from pymses.utils import constants as C

# Ramses data
ioutput = 193
ro = RamsesOutput("/data/Aquarius/output/", ioutput)

# Map operator : mass-weighted density map
up_func = lambda dset: (dset["rho"]**2)
down_func = lambda dset: (dset["rho"])
scal_op = FractionOperator(up_func, down_func)

# Map region
center = [ 0.567811, 0.586055, 0.559156 ]
axes = {"los": array([ -0.172935, 0.977948, -0.117099 ])}

# Map processing
rt = raytracing.RayTracer(ro, ["rho"])
for axname, axis in axes.items():
	cam  = Camera(center=center, line_of_sight_axis=axis, up_vector="z", region_size=[3.0E-2, 3.0E-2], \
			distance=2.0E-2, far_cut_depth=2.0E-2, map_max_size=512)
	map = rt.process(scal_op, cam)
	factor = ro.info["unit_density"].express(C.H_cc)
	scale = ro.info["unit_length"].express(C.Mpc)

	# Save map into HDF5 file
	mapname = "gas_rt_mw_%s_%5.5i"%(axname, ioutput)
	h5fname = save_map_HDF5(map, cam, map_name=mapname)

	# Plot map into Matplotlib figure/PIL Image
	fig = save_HDF5_to_plot(h5fname, map_unit=("H/cc",factor), axis_unit=("Mpc", scale), cmap="jet")
#	pil_img = save_HDF5_to_img(h5fname, cmap="jet")
	
	# Save into PNG image file
#	save_HDF5_to_plot(h5fname, map_unit=("H/cc",factor), axis_unit=("Mpc", scale), img_path="./", cmap="jet")
#	save_HDF5_to_img(h5fname, img_path="./", cmap="jet")

#pylab.show()