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 : max. AMR level of refinement along the line-of-sight
scal_op = MaxLevelOperator()

# 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=[4.0E-2, 4.0E-2], \
			distance=2.0E-2, far_cut_depth=2.0E-2, map_max_size=512, log_sensitive=False)
	map = rt.process(scal_op, cam)
	scale = ro.info["unit_length"].express(C.Mpc)

	# Save map into HDF5 file
	mapname = "gas_rt_lmax_%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=("AMR level",1.0), axis_unit=("Mpc", scale), cmap="jet", discrete=True)
#	pil_img = save_HDF5_to_img(h5fname, cmap="jet", discrete=True)
	
	# Save into PNG image file
#	save_HDF5_to_plot(h5fname, map_unit=("AMR level",1.0), axis_unit=("Mpc", scale), img_path="./", cmap="jet", discrete=True)
#	save_HDF5_to_img(h5fname, img_path="./", cmap="jet", discrete=True)

#pylab.show()