pipeline/disk_postprocess.py

# coding: utf-8
import sys
import os

import pymses
import numpy as np
import matplotlib as mpl

if os.environ.get("DISPLAY", "") == "":
    print("No display found. Using non-interactive Agg backend")
    mpl.use("Agg")

import pylab as P
import glob as glob

try:
    import cPickle as pickle
except:
    print("cPickle not found, using pickle")
    import pickle
import tables

from pymses.sources.ramses import output
from pymses.analysis import Camera, raytracing, slicing, splatting
from pymses.filters import CellsToPoints
from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator

# extension for out files
out_ext = ".jpeg"

P.rcParams["image.cmap"] = "plasma"
P.rcParams["savefig.dpi"] = 400


def make_image_disk(
    path,
    num,
    path_out=None,
    order="<",
    save_data=True,
    force=False,
    tag="",
    vel_red=20,
    map_size=512,
    put_title=True,
    cpu=False,
    level=False,
    pos_star=np.array([1.0, 1.0, 1.0]),
    interactive=False,
    fft=False,
):
    """
    Make several useful image of an output of a simulation

    Parameters
    ----------
    path      path of the Ramses output
    num       Ramses output number
    path_out  path of the pipeline outputb
    order     '<' or '>' order used by pymses for reading ramses output
    force     if set, erase any existing pipeline output files
    tag       string to add to the output name
    vel_red   number of point where velocity should be plot in the slices
    map_size  size of the map
    cpuamr    plot also levels and cpus at each step
    """

    ro = pymses.RamsesOutput(path, num, order=order)
    amr = ro.amr_source(["rho", "vel", "P"])
    rad = 0.5
    center = [0.5, 0.5, 0.5]

    return make_image_aux(
        amr,
        ro,
        center,
        rad,
        num,
        path,
        force=force,
        path_out=path_out,
        save_data=save_data,
        map_size=map_size,
        vel_red=vel_red,
        tag=tag,
        cpu=cpu,
        level=level,
        put_title=put_title,
        pos_star=pos_star,
        interactive=interactive,
        fft=fft,
    )


def make_image_aux(
    amr,
    ro,
    center,
    radius,
    num,
    path,
    force=False,
    path_out=None,
    save_data=True,
    map_size=512,
    vel_red=20,
    tag="",
    cpu=False,
    level=False,
    pos_star=np.array([1.0, 1.0, 1.0]),
    put_title=True,
    interactive=False,
    fft=False,
):
    """
    Make several useful image of an output of a simulation, auxillary function

    Parameters
    ----------
    amr
    ro        pymses.RamsesOutput object

    center    3D array for coordinates center
    num       output number
    path_out  path of the pipeline output
    force     if set, erase any existing pipeline output files
    tag       string to add to the output name
    vel_red   number of point where velocity should be plot in the slices
    map_size  size of the map
    """

    lbox = ro.info["boxlen"]  # boxlen in codeunits
    G = 1.0  # Gravitational constant

    # Plotting parameters
    ntick = 6
    title_ax = {"x": "x (code)", "y": "y (code)", "z": "z (code)"}
    im_extent = [
        (-radius + center[0]) * lbox,
        (radius + center[0]) * lbox,
        (-radius + center[1]) * lbox,
        (radius + center[1]) * lbox,
    ]

    time = ro.info["time"]  # time in codeunits
    title = "t=" + str(time)[0:5] + " (code)"

    # Determining outpout directory
    if path_out is not None:
        directory = path_out
    else:
        directory = path

    # Checking for existing files
    name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext
    if len(glob.glob(name)) == 1 and not force:
        return

    # Prepare saving data
    if save_data:
        maps_disk = {"time": time, "im_extent": im_extent}

    rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"])

    rt = None
    if fft:
        rt = splatting.SplatterProcessor(amr, ro.info, rho_op)
    else:
        rt = raytracing.RayTracer(amr, ro.info, rho_op)

    axes_los = ["x", "y", "z"]  # Line of sight axes
    axes_h = ["y", "x", "x"]  # Horizontal axes
    axes_v = ["z", "z", "y"]  # Vertical axes

    ax_nb = {"x": 0, "y": 1, "z": 2}

    for i, ax_los in enumerate(axes_los):
        ax_h = axes_h[i]
        ax_v = axes_v[i]

        cam = Camera(
            center=center,
            line_of_sight_axis=ax_los,
            region_size=[2.0 * radius, 2.0 * radius],
            distance=radius,
            far_cut_depth=radius,
            up_vector=ax_v,
            map_max_size=map_size,
        )

        if interactive:
            P.figure()
        else:
            P.close()

        # Levels
        if level:
            amr.set_read_levelmax(20)
            level_op = MaxLevelOperator()
            rt_level = raytracing.RayTracer(amr, ro.info, level_op)
            datamap = rt_level.process(cam, surf_qty=True)
            map_level = datamap.map.T
            # if save_data:
            #   maps_disk['levels_' + ax_los] = map_level

            levels_ar = np.arange(ro.info["levelmin"], ro.info["levelmax"] + 1)

            # Computing linewidths
            lw = np.ones(levels_ar.size) * 2
            lvl_th = 8  # Level threeshold for reducing linewidths
            lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th] ** (
                lvl_th - levels_ar[levels_ar >= lvl_th]
            )
            lw[levels_ar < lvl_th] = 1.0

            cont = P.contour(
                map_level,
                extent=im_extent,
                origin="lower",
                colors="k",
                linewidths=lw,
                levels=levels_ar,
            )
            cont.levels = cont.levels + 1

            P.clabel(
                cont,
                levels_ar[levels_ar < lvl_th + 2][1::2],
                inline=1,
                fontsize=8.0,
                fmt="%1d",
            )

        # Column density
        datamap = rt.process(cam, surf_qty=True)
        dmap_col = datamap.map.T * lbox
        map_col = np.log10(dmap_col)

        if save_data:
            maps_disk["coldens_" + ax_los] = dmap_col

        im = P.imshow(map_col, extent=im_extent, origin="lower")

        P.locator_params(axis=ax_h, nbins=ntick)
        P.locator_params(axis=ax_v, nbins=ntick)

        if put_title:
            P.title(title)

        P.xlabel(title_ax[ax_h])
        P.ylabel(title_ax[ax_v])

        cbar = P.colorbar(im)
        cbar.set_label(r"$log(N)$ (code)")
        name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05")
        name_im = name + out_ext

        if interactive:
            P.figure()
        else:
            P.savefig(name_im)
            P.close()

        # Rho slice
        datamap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0)
        dmap_rho = (datamap_rho.map).T
        map_rho = np.log10(dmap_rho)

        vh_op = ScalarOperator(
            lambda dset: dset["vel"][:, ax_nb[ax_h]], ro.info["unit_velocity"]
        )
        dmap_vh = slicing.SliceMap(amr, cam, vh_op, z=0.0)
        map_vh_red = dmap_vh.map[
            ::vel_red, ::vel_red
        ]  # take only a subset of velocities
        map_vh_red = map_vh_red.T

        vv_op = ScalarOperator(
            lambda dset: dset["vel"][:, ax_nb[ax_v]], ro.info["unit_velocity"]
        )
        dmap_vv = slicing.SliceMap(amr, cam, vv_op, z=0.0)
        map_vv_red = dmap_vv.map[::vel_red, ::vel_red]
        map_vv_red = map_vv_red.T

        # if save_data:
        #     maps_disk['rho_' + ax_los] = dmap_rho
        #     maps_disk['v' + ax_h + '_' + ax_los] = (dmap_vh.map).T
        #     maps_disk['v' + ax_v + '_' + ax_los] = (dmap_vv.map).T

        im = P.imshow(map_rho, extent=im_extent, origin="lower")
        P.locator_params(axis=ax_h, nbins=ntick)
        P.locator_params(axis=ax_v, nbins=ntick)
        nh = map_vh_red.shape[0]
        nv = map_vv_red.shape[1]
        vec_h = (
            np.arange(nh) * 2.0 / nh * radius - radius + center[0] + radius / nh
        ) * lbox
        vec_v = (
            np.arange(nv) * 2.0 / nv * radius - radius + center[1] + radius / nv
        ) * lbox
        hh, vv = np.meshgrid(vec_h, vec_v)
        max_v = np.max(np.sqrt(map_vh_red ** 2 + map_vv_red ** 2))

        Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units="width")
        P.quiverkey(
            Q,
            0.7,
            0.95,
            max_v,
            r"$" + str(max_v)[0:4] + "$ (code)",
            labelpos="E",
            coordinates="figure",
        )
        if put_title:
            P.title(title)
        P.xlabel(title_ax[ax_h])
        P.ylabel(title_ax[ax_v])
        cbar = P.colorbar(im)
        cbar.set_label(r"$log(n)$ (code)")

        name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05")
        name_im = name + out_ext
        if interactive:
            P.figure()
        else:
            P.savefig(name_im)
            P.close()

        P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"])
        dmap_P = (slicing.SliceMap(amr, cam, P_op, z=0.0)).map.T

        dmap_T = dmap_P / dmap_rho
        map_T = np.log10(dmap_T)

        # if save_data:
        #     maps_disk['T_' + ax_los] = dmap_T

        im = P.imshow(map_T, extent=im_extent, origin="lower")

        P.locator_params(axis="x", nbins=ntick)
        P.locator_params(axis="y", nbins=ntick)

        if put_title:
            P.title(title)
        P.xlabel(title_ax[ax_h])
        P.ylabel(title_ax[ax_v])
        cbar = P.colorbar(im)
        cbar.set_label(r"$log(T) \, (K)$")

        name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05")
        name_im = name + out_ext
        if interactive:
            P.figure()
        else:
            P.savefig(name_im)
            P.close()

        # Toomre parameter

        if ax_los == "z":

            def omega_rho_func(dset):
                pos = dset.get_cell_centers()
                pos = pos - (pos_star / lbox)

                xx = pos[:, :, 0]
                yy = pos[:, :, 1]
                rc = np.sqrt(xx ** 2 + yy ** 2)  # cylindrical radius
                vx = dset["vel"][:, :, 0]
                vy = dset["vel"][:, :, 1]
                omega_rho = 1.0 / rc ** 2
                omega_rho = omega_rho * dset["rho"]
                vyx = vy * xx
                vxy = vx * yy
                omega_rho = omega_rho * (vyx - vxy)
                return omega_rho

            omega_op = FractionOperator(
                omega_rho_func, lambda dset: dset["rho"], 1.0 / ro.info["unit_time"]
            )
            cs_op = FractionOperator(
                lambda dset: dset["P"],
                lambda dset: dset["rho"],
                ro.info["unit_velocity"],
            )
            rt_omega = raytracing.RayTracer(amr, ro.info, omega_op)

            if fft:
                rt_cs = splatting.SplatterProcessor(amr, ro.info, cs_op, surf_qty=False)
            else:
                rt_cs = raytracing.RayTracer(amr, ro.info, cs_op)

            dmap_omega = rt_omega.process(cam)
            dmap_cs = rt_cs.process(cam)
            map_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col)

            # if save_data:
            #     maps_disk['Q_' + ax_los] = map_Q

            im = P.imshow(
                map_Q, extent=im_extent, origin="lower", norm=mpl.colors.LogNorm()
            )

            P.locator_params(axis="x", nbins=ntick)
            P.locator_params(axis="y", nbins=ntick)

            if put_title:
                P.title(title)
            P.xlabel(title_ax[ax_h])
            P.ylabel(title_ax[ax_v])
            cbar = P.colorbar(im)
            cbar.set_label(r"$Q$")

            name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05")
            name_im = name + out_ext
            if interactive:
                P.figure()
            else:
                P.savefig(name_im)
                P.close()

        if cpu:
            cpu_op = ScalarOperator(
                lambda dset: dset.icpu * (np.ones(dset["P"].shape)),
                ro.info["unit_pressure"],
            )
            rt_cpu = raytracing.RayTracer(amr, ro.info, cpu_op)
            datamap = rt_cpu.process(cam, surf_qty=True)
            map_cpu = datamap.map.T

            # if save_data:
            #     maps_disk['cpu_' + ax_los] = map_cpu

            im = P.imshow(map_cpu, extent=im_extent, origin="lower")
            P.locator_params(axis="x", nbins=ntick)
            P.locator_params(axis="y", nbins=ntick)

            if put_title:
                P.title(title)

            P.xlabel(title_ax[ax_h])
            P.ylabel(title_ax[ax_v])
            cbar = P.colorbar(im)
            cbar.set_label(r"cpu")

            name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05")
            name_im = name + out_ext

            if interactive:
                P.figure()
            else:
                P.savefig(name_im)
                P.close()

    if save_data:
        name_save = (
            directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save"
        )
        f = open(name_save, "w")
        pickle.dump(maps_disk, f)
        f.close()
    return maps_disk


def disk_prop(
    path_in,
    num,
    path_out=None,
    force=False,
    nb_bin=20,
    rad_ext=1.0,
    mass_star=1.0,
    pos_star=np.array([1.0, 1.0, 1.0]),
    binning="log",
):
    """
    Compute properties of a disk in the plane (0,x,y)
    with a protostar at the center of the box

    The region of the disk is defined by its scale height

    Parameters
    ----------

    path_in   path of the input data files (output of ramses)
    num       id of the output
    path_out  optional path to the output files
    force     if set, redo ouptut even if already done

    nb_bin    Number of radial bins
    rad_ext   Outer radius of the disk

    pos_star  position of the central protostar
    mass_star mass of the central protostar
    """

    # Set th output directory
    if path_out is not None:
        directory_out = path_out
    else:
        directory_out = path_in

    # Check if the output file exists, and exit if it is the case
    name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save"
    if not force and len(glob.glob(name_save)) != 0:
        return

    # Compute the bins array
    if binning == "log":
        lrad = np.log10(rad_ext)
        rad = np.logspace(lrad - 2.0, lrad, num=nb_bin)
    elif binning == "lin":
        rad = np.linspace(0.0, rad_ext, num=nb_bin)
    else:
        raise ValueError(
            "Incorrect binning specification, binnning should be 'lin' or 'log'"
        )

    # Get Ramses data
    ro = pymses.RamsesOutput(path_in, num)
    lbox = ro.info["boxlen"]  # boxlen in codeunits (=>pc)

    time = ro.info["time"]  # time in codeunits

    # Get array of cell positions
    amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P", "g", "phi"])
    cell_source = CellsToPoints(amr)
    cells = cell_source.flatten()
    dx = cells.get_sizes() * lbox
    pos = cells.points * lbox
    # Get positions in the frame of the protostar
    pos = pos - pos_star

    # Get cylindrical radius
    rc = np.sqrt(pos[:, 0] ** 2 + pos[:, 1] ** 2)
    # Get velocities
    vel = cells["vel"]
    # Get radial component of velocity
    norm_pos = rc
    norm_pos[rc == 0] = 1.0e-10  # Avoid division per 0
    v_rad = (pos[:, 0] * vel[:, 0] + pos[:, 1] * vel[:, 1]) / norm_pos
    # Get azimuthal component of velocity
    v_az = (pos[:, 0] * vel[:, 1] - pos[:, 1] * vel[:, 0]) / norm_pos
    # Gravitational field
    g = cells["g"]
    g_rad = (pos[:, 0] * g[:, 0] + pos[:, 1] * g[:, 1]) / norm_pos
    g_az = (pos[:, 0] * g[:, 1] - pos[:, 1] * g[:, 0]) / norm_pos

    # Select cells that are actually in the disk, ie within the scale height
    G = 1.0
    cs = np.sqrt(cells["P"] / cells["rho"])  # sound velocity
    height = cs * np.sqrt(rc ** 3 / (G * mass_star))
    mask_pos = np.abs(pos[:, 2]) < height  # condition on position
    mask_dens = cells["rho"] > 1.0e6  # condition on density
    mask = mask_pos  # & mask_dens
    print("Number of selected cells ", np.sum(mask))

    pos_disk = pos[mask]
    rc_disk = rc[mask]
    vel_disk = vel[mask]
    rho_disk = cells["rho"][mask]  # density
    press_disk = cells["P"][mask]  # pressure
    dx_disk = dx[mask]
    dvol_disk = dx_disk ** 3
    v_rad_disk = v_rad[mask]
    v_az_disk = v_az[mask]
    v_kepl = np.sqrt(mass_star * G / rc_disk)
    height_disk = height[mask]
    g_rad_disk = g_rad[mask]
    g_az_disk = g_az[mask]

    total_mass_disk = np.sum(rho_disk * dvol_disk)
    total_mass = np.sum(cells["rho"] * dx ** 3)

    print("Mass disk", total_mass_disk)
    print("Mass box", total_mass)

    # Initialize binned quantities
    cs_rad = np.zeros(nb_bin - 1)
    temp_rad = np.zeros(nb_bin - 1)
    press_rad = np.zeros(nb_bin - 1)
    rho_rad = np.zeros(nb_bin - 1)
    coldens_rad = np.zeros(nb_bin - 1)
    v_az_rad = np.zeros(nb_bin - 1)
    v_kepl_rad = np.zeros(nb_bin - 1)
    v_rad_rad = np.zeros(nb_bin - 1)
    alpha_rey_rad = np.zeros(nb_bin - 1)
    alpha_rey_rad_bis = np.zeros(nb_bin - 1)
    alpha_grav_rad = np.zeros(nb_bin - 1)
    Q_kepl_rad = np.zeros(nb_bin - 1)
    height_rad = np.zeros(nb_bin - 1)
    vol_rad = np.zeros(nb_bin - 1)  # Volume of a bin
    surf_rad = np.zeros(nb_bin - 1)  # Surface of a bin
    mass_rad = np.zeros(nb_bin - 1)  # Mass of a bin

    for i in range(nb_bin - 1):
        mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1])

        # print("Bin #{} : {} cells between {} and {}".format(i, np.sum(mask_bin), rad[i], rad[i + 1]))

        vol_rad[i] = np.sum(dvol_disk[mask_bin])
        mass_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
        press_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i]
        rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i]
        temp_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i]

        # Surface of a bin :  S = dr * 2 * pi * r with
        #  dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2.
        surf_rad[i] = (rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi
        coldens_rad[i] = mass_rad[i] / surf_rad[i]

        v_az_rad[i] = (
            np.sum(v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
            / mass_rad[i]
        )

        v_rad_rad[i] = (
            np.sum(v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
            / mass_rad[i]
        )

        height_rad[i] = (
            np.sum(height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
            / mass_rad[i]
        )

        alpha_rey_rad[i] = (2.0 / 3) * (
            (
                np.sum(
                    v_az_disk[mask_bin]
                    * v_rad_disk[mask_bin]
                    * rho_disk[mask_bin]
                    * dvol_disk[mask_bin]
                )
                / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin])
                - v_az_rad[i] * v_rad_rad[i] * rho_rad[i] / press_rad[i]
            )
            * v_az_rad[i]
            / abs(v_az_rad[i])
        )

        alpha_grav_rad[i] = (2.0 / 3) * (
            np.sum(
                g_az_disk[mask_bin]
                * g_rad_disk[mask_bin]
                * rho_disk[mask_bin]
                * dvol_disk[mask_bin]
            )
            / (4 * np.pi * G)
            / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin])
            * coldens_rad[i]
        )

        v_kepl_rad[i] = (
            np.sum(v_kepl[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin])
            / mass_rad[i]
        )

    # Derived quantities
    cs_rad = np.sqrt(temp_rad)
    Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1])

    # Means
    mask_mean = (0.1 < rad[0 : nb_bin - 1]) & (rad[0 : nb_bin - 1] < 0.2)
    mass_mean = np.sum(mass_rad[mask_mean])
    alpha_rey_mean = np.sum(alpha_rey_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
    alpha_grav_mean = (
        np.sum(alpha_grav_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
    )
    Q_mean = np.sum(Q_kepl_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
    Q_min = np.nanmin(Q_kepl_rad)

    # store the results
    prop_disk = {
        "time": time,
        "mass_disk": total_mass_disk,
        "mass_box": total_mass,
        "rad": rad[0 : nb_bin - 1],
        "center": pos_star,
        "alpha_rey": alpha_rey_rad,
        "alpha_rey_mean": alpha_rey_mean,
        "alpha_grav": alpha_grav_rad,
        "alpha_grav_mean": alpha_grav_mean,
        "v_rad": v_rad_rad,
        "v_az": v_az_rad,
        "v_kepl": v_kepl_rad,
        "coldens": coldens_rad,
        "rho": rho_rad,
        "press": press_rad,
        "temp": temp_rad,
        "cs": cs_rad,
        "Q_kepl": Q_kepl_rad,
        "Q_mean": Q_mean,
        "Q_min": Q_min,
        "height": height_rad,
    }
    f = open(name_save, "w")
    pickle.dump(prop_disk, f)
    f.close()


def plot_disk_prop(path, num, force=False, tag="", interactive=False):
    """
    Plot properties of a disk

    num       id of the ramses output
    path      path to the properties file
    force     if set, redo plots even if already done
    """

    # Load property file
    name_save = path + "/prop_disk_" + str(num).zfill(5) + ".save"

    # Check if the properties file exists
    if len(glob.glob(name_save)) == 0:
        raise ("no pickle file for disk properties. Run disk_prop() first")
    f = open(name_save, "r")
    prop_disk = pickle.load(f)
    f.close()

    # Check if the output file exists, and exit if it is the case
    name_save = path + "/rho_disk_r_" + str(num).zfill(5) + out_ext
    if not force and len(glob.glob(name_save)) != 0:
        return

    time = prop_disk["time"]
    mass = prop_disk["mass_disk"]
    title = "t=" + str(time)[0:5] + " (code)"
    rad = prop_disk["rad"]

    if interactive:
        P.figure()
    else:
        P.close()

    P.xscale("log")
    P.yscale("log")
    P.grid()
    P.plot(rad, prop_disk["rho"], color="k", linewidth=2)
    P.ylabel(r"$n \, (code)$")
    P.xlabel("disk radius")
    P.title(title)
    if interactive:
        P.figure()
    else:
        P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext)
        P.close()

    P.xscale("log")
    P.yscale("log")
    P.grid()
    P.plot(rad, prop_disk["temp"], color="k", linewidth=2)
    P.ylabel(r"$T \, (K)$")
    P.xlabel("disk radius")
    P.title(title)
    if interactive:
        P.figure()
    else:
        P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + out_ext)
        P.close()

    P.xscale("log")
    P.yscale("symlog", linthreshy=0.01)

    P.plot(rad, prop_disk["v_rad"], color="k", linewidth=2)
    P.plot(rad, prop_disk["v_kepl"], color="b", linewidth=2)
    P.plot(rad, abs(prop_disk["v_az"]), color="r", linewidth=2)
    P.plot(rad, prop_disk["cs"], color="c", linewidth=2)
    P.grid()
    P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right")

    P.ylabel(r"$V \, (km s^{-1})$")
    P.xlabel("disk radius")
    if interactive:
        P.figure()
    else:
        P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext)
        P.close()

    P.xscale("log")
    P.yscale("log")
    P.grid()
    P.plot(rad, prop_disk["coldens"], color="k", linewidth=2)
    P.ylabel(r"$N\, (cm^{-2})$")
    P.xlabel("disk radius ")
    P.title(title)
    if interactive:
        P.figure()
    else:
        P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext)
        P.close()

    # Alpha
    alpha_rey_mean, alpha_grav_mean = (
        prop_disk["alpha_rey_mean"],
        prop_disk["alpha_grav_mean"],
    )

    # P.xscale('log')
    P.xlim([1e-2, 0.25])
    P.yscale("log")
    P.ylim([1e-7, 1.0])
    P.grid()
    P.plot(
        rad, abs(prop_disk["alpha_rey"]), "b", linewidth=2, label=r"$\alpha_{Reynolds}$"
    )
    P.plot(rad, abs(alpha_rey_mean * np.ones(len(rad))), "b:", linewidth=1)
    P.plot(
        rad, abs(prop_disk["alpha_grav"]), "r", linewidth=2, label=r"$\alpha_{grav}$"
    )
    P.plot(rad, abs(alpha_grav_mean * np.ones(len(rad))), "r:", linewidth=1)
    P.plot(
        rad,
        abs(prop_disk["alpha_rey"]) + abs(prop_disk["alpha_grav"]),
        "g--",
        linewidth=2,
        label=r"$\alpha_{tot}$",
    )
    P.title(
        title
        + r", $\bar{\alpha}_{Reynolds} = %.1e,  \bar{\alpha}_{grav} = %.1e$"
        % (alpha_rey_mean, alpha_grav_mean)
    )

    P.legend()
    P.ylabel(r"$\alpha$")
    P.xlabel("disk radius ")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext)
        P.close()

    # Q
    P.ylim([0, 10.0])
    P.xlim([0, 0.5])
    P.yticks(np.arange(0.0, 11, 1.0))
    P.grid()
    P.plot(rad, abs(prop_disk["Q_kepl"]), color="b", linewidth=2)
    P.plot(rad, abs(prop_disk["Q_mean"]) * np.ones(len(rad)), "b:", linewidth=1)
    P.ylabel(r"$Q$")
    P.xlabel("disk radius ")
    P.title(title + ", mass of disk = {} (code)".format(mass))

    if interactive:
        pass
    else:
        P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext)
        P.close()

    # height ratio
    P.grid()
    P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2)
    P.ylabel(r"H ratio")
    P.xlabel("disk radius ")
    P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"]))

    if interactive:
        pass
    else:
        P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext)
        P.close()


def disk_pdf(
    path,
    num,
    maps_disk,
    pos_star=[1.0, 1.0],
    force=False,
    interactive=False,
    nb_bin_hist=50,
    tag="",
    rad_min=0.1,
):

    # Load property file
    name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save"

    # Check if the properties file exists
    if len(glob.glob(name_prop)) == 0:
        raise IOError("no pickle file for disk properties. Run disk_prop() first")
    f = open(name_prop, "r")
    prop_disk = pickle.load(f)
    f.close()

    # Load maps file
    print("load maps file")
    name_maps = path + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save"

    # Check if the maps file exists
    if maps_disk is None:
        if len(glob.glob(name_maps)) == 0:
            raise IOError("no pickle file for disk maps. Run make_image_disk() first")
        f = open(name_maps, "r")
        maps_disk = pickle.load(f)
        f.close()
    print("maps file loaded")

    time = prop_disk["time"]
    title = "t=" + str(time)[0:5] + " (code)"

    coldens_map = maps_disk["coldens_z"]
    im_extent = maps_disk["im_extent"]
    rad_bins = prop_disk["rad"]
    coldens_rad = prop_disk["coldens"]

    x = np.linspace(im_extent[0], im_extent[1], coldens_map.shape[0])
    y = np.linspace(im_extent[2], im_extent[3], coldens_map.shape[1])
    xx, yy = np.meshgrid(x, y)
    rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2)

    # Find appropriate bin for each coordinate set
    bins = np.zeros(rr.shape, dtype=int)
    for r in rad_bins[1:]:
        bins = bins + (rr >= r).astype(int)

    coldens_mean = coldens_rad[bins.flatten()]
    coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape)

    #  Kill fluctuations for r < rad_min
    coldens_map[rr < rad_min] = coldens_mean_map[rr < rad_min]

    # Histogramm : density fluctuation distribution function
    dcoldens = np.log10(coldens_map.flatten() / coldens_mean)
    P.grid()
    P.yscale("log")
    P.xlabel(r"$\log(N / \bar{N})$")
    P.ylabel("Probability density")
    P.title(title)
    P.hist(dcoldens, nb_bin_hist, range=(-1, 2))

    if interactive:
        pass
    else:
        P.savefig(path + "/dcol_hist_" + str(num).zfill(5) + out_ext)
        P.close()

    # Fluctuation map
    dcoldens_map = np.reshape(dcoldens, coldens_map.shape)

    # cont = P.contour(coldens_mean_map,
    #                  extent=im_extent,
    #                  origin='lower',
    #                  colors='k',
    #                  linewidths=0.1)

    im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis")

    P.title(title)
    P.xlabel(r"$x$")
    P.ylabel(r"$y$")

    cbar = P.colorbar(im)
    cbar.set_label(r"$log(N/\bar{N})$ (code)")
    name = path + "/dcoldensmap_" + format(num, "05")
    name_im = name + out_ext

    if interactive:
        P.figure()
    else:
        P.savefig(name_im)
        P.close()


def compare(path, runs, num, force=False, interactive=False):
    """
    Compare properties of a disk in several simulations

    num         id of the ramses output
    runs        list of runs to consider
    path        path to the properties file
    force       if set, redo plots even if already done
    interactive interactive mode, to use in a %pylab ipython shell
    """

    # Initialize arrays
    time = np.zeros(len(runs))
    alpha_rey = np.zeros(len(runs))
    alpha_grav = np.zeros(len(runs))
    Q = np.zeros(len(runs))
    Q0 = np.zeros(len(runs))

    for i, run in enumerate(runs):
        path_run = path + "/" + run
        # Load property file
        name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save"

        # Check if the properties file exists
        if len(glob.glob(name_save)) == 0:
            raise ("no pickle file for disk properties. Run disk_prop() first")
        f = open(name_save, "r")
        prop_disk = pickle.load(f)
        f.close()

        time[i] = prop_disk["time"]
        alpha_rey[i] = prop_disk["alpha_rey_mean"]
        alpha_grav[i] = prop_disk["alpha_grav_mean"]
        Q[i] = prop_disk["Q_min"]
        Q0[i] = float(run.split("_")[2][1:])

    # Check if the output file exists, and exit if it is the case
    name_save = path + "/alphaQ_" + str(num).zfill(5) + out_ext
    if not force and len(glob.glob(name_save)) != 0:
        return

    title = "t=" + str(time[0]) + " (code)"

    # alpha = f(Qmin)
    P.yscale("log")
    P.ylim([1e-7, 1.0])
    P.grid()

    P.plot(Q, abs(alpha_rey), "o--", label=r"$\bar{\alpha}_{Reynolds}$")
    P.plot(Q, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")

    P.legend()
    P.ylabel(r"$\bar{\alpha}$")
    P.xlabel(r"$Q_{min}$")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/alphaQ_" + str(num).zfill(5) + out_ext)
        P.close()

    # alpha = f(Q0)
    P.yscale("log")
    P.ylim([1e-7, 1.0])
    P.grid()

    P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
    P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")

    P.legend()
    P.ylabel(r"$\bar{\alpha}$")
    P.xlabel(r"$Q_{0}$")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext)
        P.close()


def evolution(path, nums, force=False, interactive=False):
    """
    Plot properties over time

    path        path to the properties file
    nums        list of id of the ramses output
    force       if set, redo plots even if already done
    interactive interactive mode, to use in a %pylab ipython shell
    """

    # Initialize arrays
    time = np.zeros(len(nums))
    alpha_rey = np.zeros(len(nums))
    alpha_grav = np.zeros(len(nums))
    Qmin = np.zeros(len(nums))
    Qmean = np.zeros(len(nums))
    mass_disk = np.zeros(len(nums))
    mass_box = np.zeros(len(nums))

    for i, num in enumerate(nums):

        # Load property file
        name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save"

        # Check if the properties file exists
        if len(glob.glob(name_prop)) == 0:
            raise ("no pickle file for disk properties. Run disk_prop() first")
        f = open(name_prop, "r")
        prop_disk = pickle.load(f)
        f.close()

        time[i] = prop_disk["time"]
        alpha_rey[i] = prop_disk["alpha_rey_mean"]
        alpha_grav[i] = prop_disk["alpha_grav_mean"]
        Qmin[i] = prop_disk["Q_min"]
        Qmean[i] = prop_disk["Q_mean"]
        mass_disk[i] = prop_disk["mass_disk"]
        mass_box[i] = prop_disk["mass_box"]

    # Check if the output file exists, and exit if it is the case
    name_save = path + "/alpha_time" + out_ext
    if not force and len(glob.glob(name_save)) != 0:
        return

    # Alpha
    P.yscale("log")
    P.ylim([1e-7, 1.0])
    P.grid()

    P.plot(time, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
    P.plot(time, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")

    P.legend()
    P.ylabel(r"$\bar{\alpha}$")
    P.xlabel(r"time (code)")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/alpha_time" + out_ext)
        P.close()

    # Q
    P.grid()
    P.plot(time, Qmin, "o-.", label=r"$Q_{min}$")
    P.plot(time, Qmean, "*--", label=r"$\bar{Q}$")

    P.legend()
    P.ylabel(r"$Q$")
    P.xlabel(r"time (code)")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/Q_time" + out_ext)
        P.close()

    # M
    P.grid()
    P.plot(time, mass_disk, "o-.", label=r"$M_{disk}$")
    P.plot(time, mass_box, "*--", label=r"$M_{box}$")

    P.legend()
    P.ylabel(r"$M / M_{*}$")
    P.xlabel(r"time (code)")

    if interactive:
        P.figure()
    else:
        P.savefig(path + "/mass_time" + out_ext)
        P.close()


def make_clump_hop(
    path_in,
    num,
    name,
    thres_dens,
    thres_level,
    pos_zoom,
    size_zoom,
    path_out=None,
    path_hop="",
    force=False,
    gcomp=True,
):
    """
    This routine use the HOP algorithm to extract clumps defined from their peaks
    as an output it provides a list of cell position ordered by the group to which they belong

    Parameters
    ----------
    path_in      is the path where the data tobe read are located
    path_out     is the path of teh directory where resulting files must be written
    num          output number
    name         a string which is used to write the names of the various files
    thres_dens   density threshold above which cells are considered
    thres_level  level threshold above which cells are considered
    pos_zoom     the center of the zoom coordinates
    size_zoom    the 3 zoom extension (x, y and z)
    """

    if path_out is not None:
        directory_out = path_out
    else:
        directory_out = path_in

    name_txt = name + ".txt"

    # check whether hop entry files have been created (not test is done on .txt only
    if len(glob.glob(name_txt)) == 0 or force:
        ro = pymses.RamsesOutput(path_in, num)
        amr = ro.amr_source(["rho"], grav_compat=gcomp)  # density only is used
        center = np.asarray(pos_zoom)
        radius = size_zoom
        min_coords = np.zeros(3)
        max_coords = np.zeros(3)
        min_coords[:] = center[:] - radius / 2.0
        max_coords[:] = center[:] + radius / 2.0

        region = Box((min_coords, max_coords))

        # region = Sphere(center,radius)
        filt_amr = RegionFilter(region, amr)
        cell_source = CellsToPoints(
            filt_amr,
        )

        # selection of the cells of interest
        def cell_selec_func(dset):
            mask1 = dset["rho"] >= thres_dens
            dx = dset.get_sizes()
            mask2 = dx <= 0.5 ** thres_level
            return mask1 * mask2

        # begin cell_selec
        cells_selec = PointFunctionFilter(cell_selec_func, cell_source).flatten()
        dx = cells_selec.get_sizes()
        ncells = cells_selec.npoints

        # fill the matrice with ID, x,y,z and masses of particles
        val_mat = np.zeros((ncells, 5))

        val_mat[:, 0] = np.arange(ncells)
        val_mat[:, 1:4] = cells_selec.points[:, 0:3]
        val_mat[:, 4] = cells_selec["rho"] * (dx ** 3)

        # write name.txt
        head = str(ncells)

        np.savetxt(
            name_txt,
            val_mat,
            fmt="%10d %.10e %.10e %.10e %.10e",
            header=head,
            delimiter=" ",
            comments=" ",
        )
        # end of creation name.txt

        # creation name.den
        f = open(name + ".den", "wb")
        f.write(pack("I", ncells))
        cells_selec["rho"].astype("f").tofile(f)
        f.close()

        print(name + ".den created")
        # end of creation name.den

    # HOP Algorithm
    print("creation of .hop and .gbound du to hop")
    fname = path_hop + name + ".txt"
    print("look for hop in ", fname)

    h = HOP(fname, path_hop)
    h.process_hop()
    print("hop grouping is finished")
    # end of HOP Algorithm

    idpart = val_mat[:, 0]
    X = val_mat[:, 1]
    Y = val_mat[:, 2]
    Z = val_mat[:, 3]
    mass = val_mat[:, 4]

    # read the gbound file to get list of particle numbers within groups
    f = open(name + ".gbound", "r")
    aline = f.readline()
    ngroups = int(aline)
    npart_v = np.zeros(ngroups, dtype=int)
    for i in range(10):
        aline = f.readline()
    for i in range(ngroups):
        aline = f.readline()
        vec = aline.split()
        igroup = int(vec[0])
        npart_v[igroup] = int(vec[1])
    f.close()

    # get the igroup array
    group_ids = h.get_group_ids()

    # sort it and apply the sorting to the coordinates
    # this means that the particules of group 1 are written first then of group 2 etc...
    ind_sort = np.argsort(group_ids)
    xx_v = X[ind_sort]
    yy_v = Y[ind_sort]
    zz_v = Z[ind_sort]
    vect_id_group = group_ids[ind_sort]  # not so useful

    # write the sorted cells
    name_save_clump = directory_out + name + ".save"

    np.savez(
        name_save_clump,
        ngroups=ngroups,
        npart_v=npart_v,
        xx_v=xx_v,
        yy_v=yy_v,
        zz_v=zz_v,
        vect_id_group=vect_id_group,
        num=num,
        name=name,
        thres_dens=thres_dens,
        thres_level=thres_level,
        pos_zoom=pos_zoom,
        size_zoom=size_zoom,
    )

    return name_save_clump