diff --git a/clump.py b/clump.py
index d56a2dd..9fd5237 100644
--- a/clump.py
+++ b/clump.py
@@ -4,27 +4,4 @@ import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.collections import EllipseCollection
 
-me.make_clump_hop('/home/nbrucy/simus/conv_disk/73_beta3_jeansrefine', 40, 'test_hop', 3, 12,
-                   [0.5, 0.5, 0.5], 1.,  path_out='/home/nbrucy/visus/conv_disk/73_beta3_jeansrefine/', path_hop=
-                   './', force=True)
-me.clump_properties('test_hop', '/home/nbrucy/simus/conv_disk/73_beta3_jeansrefine', 40,  path_out='/home/nbrucy/visus/conv_disk/73_beta3_jeansrefine/', gcomp=False)
 
-f = open('/home/nbrucy/visus/conv_disk/73_beta3_jeansrefine/test_hop_prop_struct.save')
-a = pickle.load(f)
-f.close()
-mask = a['n_max'] > 10**4
-p = a['pos_n_max'][mask]
-c = a['n_max'][mask]
-f = open('/home/nbrucy/visus/conv_disk/73_beta3_jeansrefine/maps_disk_73_beta3_jeansrefine_00040.save')
-m = pickle.load(f)
-plt.imshow(np.log10(m['coldens_z']), extent=m['im_extent'], cmap='plasma')
-plt.scatter(p[:, 0], 2 - p[:, 1], c = c)
-
-fig, ax = plt.subplots(1, 1)
-offsets = list(zip(p[:, 0], 2 - p[:, 1]))
-size = a['size_iner2'][mask]
-
-ax.add_collection(EllipseCollection(widths=size, heights=size, angles=0, units='xy',
-                                    facecolors=plt.cm.hsv(a['n_max'][mask]),
-                                    offsets=offsets, transOffset=ax.transData))
- 
diff --git a/disk_postprocess.py b/disk_postprocess.py
index 05ccaf6..14c6f86 100644
--- a/disk_postprocess.py
+++ b/disk_postprocess.py
@@ -20,6 +20,7 @@ except:
     import pickle
 import tables
 from scipy.stats import linregress
+from numpy.polynomial.polynomial import polyfit
 from pymses.sources.ramses import output
 from pymses.sources.hop.file_formats import *
 from pymses.analysis import Camera, raytracing, slicing, splatting
@@ -27,29 +28,59 @@ from pymses.filters import CellsToPoints
 from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator
 from scipy.stats import gaussian_kde
 
+
+import module_extract as me
+
+
+
 # extension for out files
-out_ext = ".jpeg"
+out_ext = '.jpeg'
 
-P.rcParams["image.cmap"] = "plasma"
-P.rcParams["savefig.dpi"] = 400
+P.rcParams['image.cmap']='plasma'
+P.rcParams['savefig.dpi']=400
 
+params= {'text.latex.preamble' : [r'\usepackage{amsmath}']}
+P.rcParams.update(params)
 
-def compute_image_data(
-    path,
-    num,
-    radius=0.5,
-    path_out=None,
-    order="<",
-    save_data=True,
-    map_size=512,
-    tag="",
-    axes_los=["x", "y", "z"],
-    images=["coldens", "rho", "speed", "Q", "T", "level", "cpu"],
-    pos_star=np.array([1.0, 1.0, 1.0]),
-    put_title=True,
-    force=False,
-    fft=False,
-):
+def get_time(path, num):
+    """
+    Return the time of the output (code units)
+
+    Parameters
+    ----------
+    num       output number
+    path_out  path of the pipeline output
+
+    Returns
+    -------
+    time      the time of the output (code units)
+    """
+    try:
+        f = open(path + '/output_' +  str(num).zfill(5) + '/info_' +  str(num).zfill(5) + '.txt')
+        for line in f:
+            ls = line.split()
+            if len(ls) > 1 and ls[0] == 'time':
+                time = float(ls[2])
+                break
+        # ro = pymses.RamsesOutput(path, num, order='>')
+        # time = ro.info['time'] # time in codeunits
+        f.close()
+        return time
+    except IOError as e:
+        print(e)
+        return np.nan
+
+def compute_image_data(path, num, radius=0.5,
+                       path_out=None, order='<',
+                       save_data=True,
+                       map_size=512,
+                       tag='',
+                       axes_los=['x', 'y', 'z'],
+                       images=['coldens', 'rho', 'speed', 'Q', 'T', 'level', 'cpu'],
+                       pos_star=np.array([1., 1., 1.]),
+                       put_title=True,
+                       force=False,
+                       fft=False):
     """
     Make several useful image of an output of a simulation, auxillary function
 
@@ -86,9 +117,9 @@ def compute_image_data(
         directory = path
 
     # Checking for existing files
-    name_save = directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save"
-    if len(glob.glob(name_save)) == 1 and not force:
-        return
+    name_save = directory + '/maps_disk' + '_' + tag + '_' + format(num,'05') + '.save'
+    if (len(glob.glob(name_save))==1 and not force):
+        return 
 
     # Prepare saving data
     if "T" in images and not "rho" in images:
@@ -231,19 +262,14 @@ def compute_image_data(
         f.close()
     return maps_disk
 
-
-def plot_maps(
-    path,
-    num,
-    force=False,
-    vel_red=40,
-    tag="",
-    images=None,
-    axes_los=None,
-    maps_disk=None,
-    interactive=False,
-    put_title=True,
-):
+def plot_maps(path, num,
+              force=False,
+              vel_red=40, tag='',
+              images=None, axes_los=None,
+              maps_disk=None,
+              interactive=False,
+              put_title=True,
+              set_lim=True) :
     """
     Make several useful image of an output of a simulation, auxillary function
 
@@ -347,30 +373,33 @@ def plot_maps(
             cbar = P.colorbar(im)
 
             if image == 'coldens':
-                cbar.set_label(r'$log(\Sigma)$ (code)')
+                cbar.set_label(r'$\Sigma$ (code)')
+
+                if set_lim:
+                    im.set_clim(0.01, 100)
 
                 if 'levels' in images:
-                      map_level = maps_disk['levels_' + ax_los]
-                      # Computing linewidths
-                      levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1)
-                      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.
+                    map_level = maps_disk['levels_' + ax_los]
+                    # Computing linewidths
+                    levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1)
+                    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.
 
-                      cont = P.contour(map_level,
-                                       extent=im_extent,
-                                       origin='lower',
-                                       colors='k',
-                                       linewidths=lw,
-                                       levels=levels_ar)
-                      cont.levels = cont.levels + 1
+                    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., fmt='%1d')
+                    P.clabel(cont,
+                             levels_ar[levels_ar < 11][1:],
+                             inline=1, fontsize=8., fmt='%1d')
             elif image == 'rho':
-                cbar.set_label(r'$log(n)$ (code)')
+                cbar.set_label(r'$\rho$ (code)')
 
                 if 'speed' in images:
                     dmap_vh =  maps_disk['v' + ax_h + '_' + ax_los]
@@ -380,36 +409,18 @@ def plot_maps(
                     map_vv_red = dmap_vv[::vel_red,::vel_red]
                     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))
+                    vec_h = (np.arange(nh)*2./nh*radius - radius + center[0] + radius/nh) * lbox
+                    vec_v = (np.arange(nv)*2./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",
-                    )
+                    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')
 
-            elif image == "T":
-                cbar.set_label(r"$log(T) \, (K)$")
-            elif image == "Q":
-                cbar.set_label(r"$Q$")
+            elif image == 'T':
+                cbar.set_label(r'$T (code)$')
+            elif image == 'Q':
+                cbar.set_label(r'$Q$')
             else:
                 cbar.set_label(image)
 
@@ -523,7 +534,7 @@ def disk_prop(
     # Select cells that are actually in the disk, ie within the scale height
     G = 1.
     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"] > 0.01 * np.mean(cells["rho"])  # condition on density
@@ -714,12 +725,12 @@ def plot_disk_prop(
     for plot in plots:
         title = "t=" + str(time)[0:5] + " (code)"
         P.grid()
-        P.xlabel('disk radius')
+        P.xlabel(r'$r$')
         if plot == 'rho':
             P.xscale('log')
             P.yscale('log')
             P.plot(rad, prop_disk['rho'], color='k', linewidth=2)
-            P.ylabel(r'$n \, (code)$')
+            P.ylabel(r'$\rho \, (code)$')
         elif plot == 'T':
             P.xscale('log')
             P.yscale('log')
@@ -733,12 +744,12 @@ def plot_disk_prop(
             P.plot(rad, abs(prop_disk['v_az']), color='r', linewidth=2)
             P.plot(rad,prop_disk['cs'], color='c', linewidth=2)
             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.ylabel(r'$V \ (code)$')
         elif plot == 'coldens':
             P.xscale('log')
             P.yscale('log')
             P.plot(rad,prop_disk['coldens'],color='k',linewidth=2)
-            P.ylabel(r'$\Sigma\, (cm^{-2})$')
+            P.ylabel(r'$\Sigma\ (code)$')
         elif plot == 'alpha':
             alpha_rey_mean, alpha_grav_mean = prop_disk['alpha_rey_mean'], prop_disk['alpha_grav_mean']
             P.xlim([1e-2, 0.25])
@@ -768,25 +779,20 @@ def plot_disk_prop(
                 label=r"$\alpha_{tot}$",
             )
             P.legend()
-            P.ylabel(r"$\alpha$")
-            title = (
-                title
-                + r", $\bar{\alpha}_{Reynolds} = %.1e,  \bar{\alpha}_{grav} = %.1e$"
-                % (alpha_rey_mean, alpha_grav_mean)
-            )
-        elif plot == "Q":
-            P.ylim([0, 3.0])
+            P.ylabel(r'$\alpha$')
+            title = title + r', $\bar{\alpha}_{Reynolds} = %.1e,  \bar{\alpha}_{grav} = %.1e$' % (alpha_rey_mean, alpha_grav_mean)
+        elif plot == 'Q':
+            P.ylim([0, 5.])
             P.xlim([0, 0.25])
-            P.yticks(np.arange(0.0, 3, 1.0))
-            P.plot(rad, abs(prop_disk["Q_kepl"]), color="b", linewidth=2)
+            P.yticks(np.arange(0., 5, 1.0))
+            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 ")
-            title = title + ", mass of disk = {} (code)".format(mass)
-        elif plot == "H":
-            P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2)
-            P.ylabel(r"H ratio")
-            title = title + ", mass of box = {} (code)".format(prop_disk["mass_box"])
+            P.ylabel(r'$Q$')
+            title = title 
+        elif plot == 'H':
+            P.plot(rad,abs(prop_disk['height']/ rad),color='b',linewidth=2)
+            P.ylabel(r'H ratio')
+            title = title + ', mass of box = {} (code)'.format( prop_disk['mass_box'])
 
         if put_title:
             P.title(title)
@@ -795,7 +801,7 @@ def plot_disk_prop(
             P.figure()
         else:
             P.tight_layout(pad=1)
-            P.savefig(path + "/" + plot + "_disk_r_" + str(num).zfill(5) + out_ext)
+            P.savefig(path + '/' + plot + '_disk_r_'+ tag + '_' + str(num).zfill(5) + out_ext)
             P.close()
 
 
@@ -867,8 +873,7 @@ def disk_pdf(
     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)
-
+    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:]:
@@ -912,6 +917,9 @@ def disk_pdf(
         for j in range(len(rad_bins) - 1):
             mean_bin[j] = np.mean(map_arr[bins == j])
 
+        # save for later
+        # prop_disk[cur_map + '_rmean'] = mean_bin
+        
         # Add value for borders
         mean_bin = np.concatenate(([mean_bin[0]], mean_bin))
 
@@ -928,8 +936,8 @@ def disk_pdf(
         # Compute fluctuation
         fluct_maps[cur_map] = map_arr / avg_maps[cur_map]
 
-    # Histogramm : density fluctuation distribution function
-    dcoldens = np.log10(fluct_maps["coldens"]).flatten()
+    # Sigma-PDF
+    dcoldens = np.log10(fluct_maps['coldens']).flatten()
 
     nb_cells = np.sum(mask_flat)
     P.grid()
@@ -950,30 +958,18 @@ def disk_pdf(
     # Variance
     var = np.var(dcoldens[mask_flat])
     # Compute the slope of the right part of the histogramm
-    mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0)
-    if np.sum(mask_fit > 0):
-        (a, b, rho, _, stderr) = linregress(
-            centers[mask_fit], np.log10(values[mask_fit])
-        )
-        P.plot(centers, 10 ** (a * centers + b), "--", linewidth=2)
-        print("a=%e, b=%e, rho=%e, var=%e" % (a, b, rho, var))
-        try:
-            beta = int(tag.split("_")[1][4:])
-        except ValueError:
+    mask_fit = (centers > 0.) & (centers < 1.25) & (values > 0)
+    if (np.sum(mask_fit > 0)):
+        (a, b, rho, _, stderr) = linregress(centers[mask_fit], np.log10(values[mask_fit]))
+        P.plot(centers, 10**(a*centers + b), '--', linewidth=2)
+        print("a=%e, b=%e, rho=%e, var=%e"% (a, b, rho, var))
+        try :
+            beta = int(tag.split('_')[1][4:])
+        except ValueError :
             beta = 0
-        fit = {
-            "beta": beta,
-            "slope": a,
-            "origin": b,
-            "correlation": rho,
-            "stderr": stderr,
-            "var": var,
-        }
-        f = open(name_prop, "w")
-        prop_disk["fit"] = fit
-        pickle.dump(prop_disk, f)
-        f.close()
-
+        fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr, 'var':var}
+        prop_disk['fit'] = fit
+        
     if interactive:
         P.figure()
     else:
@@ -981,8 +977,31 @@ def disk_pdf(
         P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext)
         P.close()
 
-    # Derived quantities
+    # rho-PDF
     drho  = fluct_maps['rho'].flatten()
+    nb_cells = np.sum(mask_flat)
+    P.grid()
+    P.yscale('log')
+    P.ylim([0.5 / nb_cells, 1.])
+    P.xlabel(r'$\log(\rho / \bar{\rho})$')
+    P.ylabel(r'$\mathcal{P}_\rho$')
+    if put_title:
+        P.title(title)
+    values, edges, _ = P.hist(drho[mask_flat],
+                              nb_bin_hist, # range=(-1, 6),
+                              weights = np.ones(nb_cells) / nb_cells)
+    centers = 0.5 * (edges[1:] + edges[:-1])
+
+    if interactive:
+        P.figure()
+    else:
+        P.tight_layout(pad=1)
+        P.savefig(path + '/drho_hist_' + tag + '_' + str(num).zfill(5) + out_ext)
+        P.close()
+
+
+
+    # Derived quantities
     dcs   = fluct_maps['cs'].flatten()
     dv    = fluct_maps['v'].flatten()
     dvaz_kepl = abs(maps['vaz'] - v_kepl) / v_kepl
@@ -991,13 +1010,15 @@ def disk_pdf(
     dmach[dmach > 10.] = 10.
     fluct_maps['mach'] = dmach
     dmach_mean = np.mean(dmach[mask_map & (dmach < 5.)].flatten())
-
-    # Fluctuations plots
-    f = open(name_prop, "w")
-    prop_disk["dmach_mean"] = dmach_mean
+    dcoldens_max = np.max(dcoldens[mask_flat])
+    
+    # Save general properties
+    prop_disk['dcoldens_max'] = dcoldens_max
+    prop_disk['dmach_mean'] = dmach_mean
     print("dmach_mean = {}".format(dmach_mean))
-    pickle.dump(prop_disk, f)
-    f.close()
+
+    
+    # Fluctuations plots
 
     # dcs = f(drho)
     P.hist2d(
@@ -1009,6 +1030,18 @@ def disk_pdf(
     P.xlabel(r"$\log(\rho / \bar{\rho})$")
     P.ylabel(r"$\log(c_s / \bar{c_s})$")
     P.ylim([-0.6, 0.6])
+    P.xlim([-1, 1.])
+
+    (a, b, rho, _, stderr) = linregress(np.log10(drho[mask_flat]), np.log10(dcs[mask_flat]))
+    P.plot(np.log10(drho[mask_flat]), a*np.log10(drho[mask_flat]) + b, '--', linewidth=2)
+    print("cs/rho a=%e, b=%e, rho=%e"% (a, b, rho))
+
+    try :
+        beta = int(tag.split('_')[1][4:])
+    except ValueError :
+        beta = 0
+    fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr}
+    prop_disk['fit_cs'] = fit
 
     P.colorbar()
     if interactive:
@@ -1102,18 +1135,17 @@ def disk_pdf(
             P.savefig(name_im)
             P.close()
 
+    # Save on disk
+    f = open(name_prop,'w')
+    pickle.dump(prop_disk, f)
+    f.close()
 
-def compare(
-    path,
-    runs,
-    nums,
-    path_out=None,
-    force=False,
-    interactive=False,
-    Q_in_name=True,
-    pdf=False,
-    gamma=5.0 / 3.0,
-):
+
+
+            
+
+def compare(path, runs, nums, path_out=None, force=False, interactive=False,
+            Q_in_name=True, pdf=False, gamma_ad=5./3., tag=''):
     """
     Compare time averaged properties of a disk in several simulations
 
@@ -1127,36 +1159,49 @@ def compare(
     if type(nums) == int:
         nums = [nums]
 
-    nums_name = "_".join(str(num).zfill(5) for num in [nums[0], nums[-1]])
+    if(type(nums)==dict):
+        nums_name = '_'
+    else:
+        nums_name = "_" +  "_".join(str(num).zfill(5) for num in [nums[0], nums[-1]])
 
     # Initialize arrays
-    alpha_rey = np.zeros(len(runs))
-    alpha_grav = np.zeros(len(runs))
-    Q = np.zeros(len(runs))
+    var_names = ['alpha_rey_mean', 'alpha_grav_mean', 'Q', 'time']
+    std_names = ['alpha_rey_mean', 'alpha_grav_mean', 'Q']
+    plot_names = ['alphaQ', 'alphaQall']
+    
     if Q_in_name:
-        Q0 = np.zeros(len(runs))
+        var_names.append('Q0')
+        plot_names.append('alphaQ0')
     if pdf:
-        beta = np.zeros(len(runs))
-        slope = np.zeros(len(runs))
-        slope_std = np.zeros(len(runs))
-        var = np.zeros(len(runs))
-        var_std = np.zeros(len(runs))
-        dmach = np.zeros(len(runs))
-        dmach_std = np.zeros(len(runs))
+        var_names = var_names + ['beta', 'kappa', 'var', 'dmach', 'gamma']
+        std_names = std_names + ['kappa', 'var', 'dmach', 'gamma']
+        plot_names = plot_names + ['alphabeta', 'kappabeta', 'gammabeta', 'varbeta', 'dmachvarall']
 
-        all_var = []
-        all_dmach = []
-        all_beta  = []
+    var_data = dict()
+    std_data = dict()
+    all_data = dict()
 
+    for name in var_names:
+        var_data[name] = np.zeros(len(runs))
+        all_data[name] = []
+
+    for name in std_names:
+        std_data[name] = np.zeros(len(runs))
+        
     for i, run in enumerate(runs):
         path_run = path + "/" + run
         nb_outputs = 0
-
-        slope_run = []
-        var_run = []
-        dmach_run = []
-
-        for num in nums:
+        
+        if(type(nums)==dict):
+            nums_run = nums[run]
+        else:
+            nums_run = nums
+        
+        run_data = dict()
+        for name in var_names:
+            run_data[name] = []
+              
+        for num in nums_run:
             try:
                 # Load property file
                 name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save"
@@ -1172,16 +1217,35 @@ def compare(
                 prop_disk = pickle.load(f)
                 f.close()
 
-                alpha_rey[i] = alpha_rey[i] + abs(prop_disk["alpha_rey_mean"])
-                alpha_grav[i] = alpha_grav[i] + abs(prop_disk["alpha_grav_mean"])
-                Q[i] = Q[i] + prop_disk["Q_min"]
-                if pdf:
-                    fit = prop_disk["fit"]
-                    beta[i] = fit["beta"]
-                    slope_run.append(fit["slope"])
-                    var_run.append(fit["var"])
-                    dmach_run.append(prop_disk["dmach_mean"])
-
+                for name in var_names:
+                    if pdf:
+                        fit    = prop_disk['fit']
+                        fit_cs = prop_disk['fit_cs']
+                    if name == 'alpha_rey_mean':
+                        data = abs(prop_disk['alpha_rey_mean'])
+                    elif name == 'alpha_grav_mean':
+                        data =  abs(prop_disk['alpha_grav_mean'])
+                    elif name == 'Q':
+                        data =  prop_disk['Q_min']
+                    elif name == 'beta':                  
+                        data = fit['beta']
+                    elif name == 'kappa':                  
+                        data = - fit['slope']
+                    elif name == 'var':                  
+                        data = fit['var']
+                    elif name == 'dmach':
+                        data = prop_disk['dmach_mean']
+                    elif name == 'gamma':
+                        data = 2 * fit_cs['slope'] + 1
+                    elif name == 'Q0':
+                        data = float(run.split('_')[2][1:])
+                        print(data)
+                    else:
+                        data = prop_disk[name]
+                        
+                    run_data[name].append(data)
+                    all_data[name].append(data)
+                    
                 nb_outputs = nb_outputs + 1
                 print(run, num, nb_outputs)
 
@@ -1189,174 +1253,155 @@ def compare(
                 print(run, num, e, nb_outputs)
 
         if nb_outputs > 0:
-            alpha_rey[i] = alpha_rey[i] / nb_outputs
-            alpha_grav[i] = alpha_grav[i] / nb_outputs
-            Q[i] = Q[i] / nb_outputs
-            if pdf:
-                slope[i] = np.mean(slope_run)
-                slope_std[i] = np.std(slope_run)
-                var[i] = np.mean(var_run)
-                var_std[i] = np.std(var_run)
-                dmach[i] = np.mean(dmach_run)
-                dmach_std[i] = np.std(dmach_run)
-                all_var = all_var + var_run
-                all_dmach = all_dmach + dmach_run
-                all_beta = all_beta + [beta[i]]*len(var_run)
+            for name in var_names:
+                var_data[name][i] = np.mean(run_data[name])
+
+            for name in std_names:
+                std_data[name][i] = np.std(run_data[name])
         else:
-            for array in [alpha_rey, alpha_grav, Q]:
-                array[i] = np.nan
-            if pdf:
-                slope[i] = np.nan
-                var[i] = np.nan
-                dmach[i] = np.nan
+            for name in var_names:
+                var_data[name][i] = np.nan
 
-        if Q_in_name:
-            Q0[i] = float(run.split("_")[2][1:])
+            for name in std_names:
+                std_data[name][i] = np.nan
 
-    if pdf:
-        all_dmach = np.array(all_dmach)
-        all_var   = np.array(all_var)
-        all_beta   = np.array(all_beta)
+    for name in var_names:
+        all_data[name] = np.array(all_data[name])
 
-    # alpha = f(Qmin)
-    P.yscale("log")
-    P.ylim([1e-7, 1.0])
-    P.grid()
+    for pname in plot_names:
 
-    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.tight_layout(pad=1)
-        P.savefig(path_out + "/alphaQ_" + nums_name + out_ext)
-        P.close()
-
-    if Q_in_name:
-        # alpha = f(Q0)
-        P.yscale("log")
-        P.ylim([1e-7, 1.0])
         P.grid()
+        print("Ploting "+ pname)
+        
+        if pname == 'alphaQ' or pname == 'alphaQ0':
+            # alpha = f(Q)
 
-        P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
-        P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
+            if pname == 'alphaQ':
+                 P.xlabel(r'$Q_{min}$')
+                 Q = var_data['Q']
+                 xerr = std_data['Q']
+            else:
+                P.xlabel(r'$Q_{0}$')
+                Q = var_data['Q0']
+                xerr = None
+                
+            P.yscale('log')
+            P.ylim([1e-6,1.])
+            P.xlim([np.nanmin(Q)- 0.2, np.nanmax(Q) + 0.2])
+            
+            P.errorbar(Q, abs(var_data['alpha_rey_mean']), xerr=xerr, yerr=std_data['alpha_rey_mean'],
+                       fmt='*--', label=r"$\bar{\alpha}_{Reynolds}$")
+            # P.errorbar(Q, abs(var_data['alpha_grav_mean']), xerr=xerr, yerr=std_data['alpha_grav_mean'],
+            #           fmt='o--', label=r"$\bar{\alpha}_{grav}$")
 
-        P.legend()
-        P.ylabel(r"$\bar{\alpha}$")
-        P.xlabel(r"$Q_{0}$")
+            P.legend()
+            P.ylabel(r'$\bar{\alpha}$')
+        elif pname == 'kappabeta':
+            P.errorbar(var_data['beta'], var_data['kappa'], yerr=std_data['kappa'], fmt='*')
+            
+            P.ylabel(r'$\kappa = - dlog\mathcal{P}_\Sigma / d\log\left(\Sigma/\overline{\Sigma}\right)$')
+            P.xlabel(r'$\beta$')
+            P.xlim([None, np.max(var_data['beta']) + 0.3])
+            (a, b, rho, _, stderr) = linregress(var_data['beta'], var_data['kappa'])
+            #P.plot(var_data['beta'], a*var_data['beta'] + b, '--', linewidth=1.5)
+            c = polyfit(var_data['beta'], var_data['kappa'], 1, w = [1.0 / ty for ty in std_data['kappa']], full=False)
+            P.plot(var_data['beta'], c[1]*var_data['beta'] + c[0], '--', linewidth=1.5)
+            print("slope #1 : a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
+            rho_p = np.sqrt(np.sum((var_data['kappa'] - c[1]*var_data['beta'] - c[0])**2))
+            print("slope #2 : a=%e, b=%e, err^2=%e"% (c[1],c[0],rho_p**2))
+        elif pname =='gammabeta':
+            # gamma = f(beta)
+            P.errorbar(var_data['beta'], var_data['gamma'], yerr=std_data['gamma'], fmt='*')
+            P.xlim([None, np.max(var_data['beta']) + 0.3])
+            P.ylabel(r'$\Gamma_{eff}$')
+            P.xlabel(r'$\beta$')
 
+            (a, b, rho, _, stderr) = linregress(var_data['beta'], var_data['gamma'])
+            # P.plot(var_data['beta'], a*var_data['beta'] + b, '--', linewidth=1.5)
+            print("gamma : a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
+        elif pname == 'varbeta':
+            P.errorbar(var_data['beta'], var_data['var'], yerr=std_data['var'], fmt='*')
+            P.ylabel(r'$Var(\log(\Sigma / \bar{\Sigma}))$')
+            P.xlabel(r'$\beta$')
+
+            (a, b, rho, _, stderr) = linregress(var_data['beta'], var_data['var'])
+            P.plot(var_data['beta'], a*var_data['beta'] + b, '--', linewidth=1.5)
+            print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
+        elif pname == 'dmachvarall':
+            mask = all_data['dmach'] < 1
+            all_var = all_data['var'][mask]
+            all_beta = all_data['beta'][mask]
+            all_dmach = all_data['dmach'][mask]
+            beta_max = np.max(var_data['beta'])
+            beta_min = np.min(var_data['beta'])
+            P.grid()
+            cmap = mpl.cm.get_cmap('RdYlBu', beta_max - beta_min + 1)
+            P.scatter(all_dmach, np.exp(all_var), c=all_beta,
+                      vmin=beta_min-0.5,
+                      vmax=beta_max+0.5,
+                      cmap=cmap)
+            cbar = P.colorbar()
+            cbar.set_ticks(var_data['beta'])
+            cbar.set_label(r'$\beta$')
+
+            P.xlabel(r'$<(v - \bar{v}) / c_s>$')
+            P.ylabel(r'$\exp(Var(\log(\Sigma / \bar{\Sigma}))$')
+        if pname == 'alphaQall' :            
+            Q         = all_data['Q'] 
+            alpha_rey = abs(all_data['alpha_rey_mean'])
+            time      = all_data['time']
+            P.yscale('log')
+            P.ylim([1e-7,1.])
+            P.scatter(Q, alpha_rey, c=time)
+            P.xlabel(r'$Q_{min}$')
+            P.ylabel(r'$\bar{\alpha}_{Reynolds}$')
+            cbar = P.colorbar()
+            cbar.set_label(r'time (code)')
+        if pname == 'alphabeta' :
+            alpha_th = (4./9.) / (gamma_ad * (gamma_ad - 1.) * var_data['beta'])
+
+            mask = var_data['beta'] >= 10
+            P.errorbar(var_data['beta'][mask], abs(var_data['alpha_rey_mean'][mask]),  yerr=std_data['alpha_rey_mean'][mask],
+                       fmt='*--', label=r"$\bar{\alpha}_{Reynolds}$")
+            # P.plot(var_data['beta'], abs(var_data['alpha_grav_mean']), '*--', label=r"$\bar{\alpha}_{grav}$")
+            P.plot(var_data['beta'][mask], abs(alpha_th[mask]), ':', label=r"$\bar{\alpha}_{th}$")
+            P.legend()
+            P.ylabel(r'$\bar{\alpha}$')
+            P.xlabel(r'$\beta$')
+
+            (a, b, rho, _, stderr) = linregress(np.log(var_data['beta']), np.log(var_data['alpha_rey_mean']))
+            print("log(alpha)= a*log(beta) + b, a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
         if interactive:
             P.figure()
         else:
             P.tight_layout(pad=1)
-            P.savefig(path_out + "/alphaQ0_" + nums_name + out_ext)
-            P.close()
-
-    if pdf:
-        # slope of the pdf = f(beta)
-        P.grid()
-        P.errorbar(beta, slope, yerr=slope_std, fmt="o")
-
-        P.legend()
-        P.ylabel(r'$d\log\mathcal{P}_\Sigma / d\log\Sigma$')
-        P.xlabel(r'$\beta$')
-
-        (a, b, rho, _, stderr) = linregress(beta, slope)
-        P.plot(beta, a * beta + b, "--", linewidth=1.5)
-        print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2))
-
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(path_out + "/dcolslopebeta_" + nums_name + out_ext)
-            P.close()
-
-        # var of the pdf = f(beta)
-        P.grid()
-        P.errorbar(beta, var, yerr=var_std, fmt="o")
-
-        P.legend()
-        P.ylabel(r'$Var(\log(\Sigma / \bar{\Sigma})$')
-        P.xlabel(r'$\beta$')
-
-        (a, b, rho, _, stderr) = linregress(beta, var)
-        P.plot(beta, a * beta + b, "--", linewidth=1.5)
-        print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2))
-
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(path_out + "/varbeta_" + nums_name + out_ext)
-            P.close()
-
-        # var = f(log(dmach)
-        mask = all_dmach < 1
-        all_var = all_var[mask]
-        all_beta = all_beta[mask]
-        all_dmach = all_dmach[mask]
-
-        P.grid()
-        cmap = mpl.cm.get_cmap('RdYlBu', np.max(beta) - np.min(beta) + 1)
-        P.scatter(all_dmach, np.exp(all_var), c=all_beta,
-                  vmin=np.min(beta)-0.5,
-                  vmax=np.max(beta)+0.5,
-                  cmap=cmap)
-        cbar = P.colorbar()
-        cbar.set_ticks(beta)
-        cbar.set_label(r'$\beta$')
-        # P.errorbar(dmach, np.exp(var), xerr=dmach_std, yerr=np.exp(var) * var_std, fmt='+')
-
-        P.xlabel(r'$<(v - \bar{v}) / c_s>$')
-        P.ylabel(r'$\exp(Var(\log(\Sigma / \bar{\Sigma}))$')
-
-        (a, b, rho, _, stderr) = linregress(all_dmach, np.exp(all_var))
-        #P.plot(all_dmach, a*all_dmach + b, '--', linewidth=1.5)
-        print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
-
-        (a, b, rho, _, stderr) = linregress(dmach, np.exp(var))
-       # P.plot(all_dmach, a*all_dmach + b, '-.', linewidth=1.5)
-        print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
-
-
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(path_out + "/dmachvar_" + nums_name + out_ext)
+            P.savefig(path_out + '/'+  pname + '_' + tag + nums_name  + out_ext)
             P.close()
 
         # alpha = f(beta)
-        P.yscale("log")
-        P.ylim([1e-7, 1.0])
-        P.grid()
+        # # P.yscale('log')
+        # P.ylim([None,0.3])
+        # P.grid()
 
-        # theoritical alpha (Gammie 2001)
-        alpha_th = (4.0 / 9.0) / (gamma * (gamma - 1.0) * beta)
+        # # var of the pdf = f(var_data['beta'])
+        # P.grid()
+        # # theoritical alpha (Gammie 2001)
+        # alpha_th = (4./9.) / (gamma_ad * (gamma_ad - 1.) * var_data['beta'])
 
-        P.plot(beta, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
-        P.plot(beta, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
-        P.plot(beta, abs(alpha_th), ":", label=r"$\bar{\alpha}_{th}$")
+        # P.plot(var_data['beta'], abs(alpha_rey_mean), 'o-.', label=r"$\bar{\alpha}_{Reynolds}$")
+        # P.plot(var_data['beta'], abs(alpha_grav), '*--', label=r"$\bar{\alpha}_{grav}$")
+        # P.plot(var_data['beta'], abs(alpha_th), ':', label=r"$\bar{\alpha}_{th}$")
 
-        P.legend()
-        P.ylabel(r"$\bar{\alpha}$")
-        P.xlabel(r"$\beta_{cool}$")
+        # P.legend()
+        # P.ylabel(r'$\bar{\alpha}$')
+        # P.xlabel(r'$\beta_{cool}$')
 
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(path_out + "/alphabeta_" + nums_name + out_ext)
-            P.close()
+        # (a, b, rho, _, stderr) = linregress(np.log(var_data['beta']), np.log(alpha_rey_mean))
+        # print("log(alpha)= a*log(beta) + b, a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
 
 
-def evolution(path, nums, force=False, interactive=False, pdf=False):
+
+def evolution(path, nums, force=False, interactive=False, pdf=False, tag=''):
     """
     Plot properties over time
 
@@ -1367,15 +1412,19 @@ def evolution(path, nums, force=False, interactive=False, pdf=False):
     """
 
     # 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))
+    var_names = ['time', 'alpha_rey_mean', 'alpha_grav_mean', 'Q_min', 'Q_mean', 'mass_disk', 'mass_box']
+    plot_names = ['alpha', 'Q', 'mass']
+    plot_labels = [r'$\alpha$', r'$Q$', 'mass (code)']
     if pdf:
-        slope = np.zeros(len(nums))
+        var_names = var_names + ['kappa', 'var', 'dmach_mean', 'dcoldens_max', 'gamma']
+        plot_names = plot_names + ['kappa', 'var', 'dmach_mean', 'dcoldens_max', 'gamma']
+        plot_labels = plot_labels + [r'$\kappa$', r'$Var(\log(\Sigma/\overline{\Sigma})$', 'dmach_mean', r'$\max\left(\log(\Sigma/\overline{\Sigma})\right)$', r'$\Gamma_{eff}$']
+
+    plot_labels = dict(zip(plot_names, plot_labels))
+    var_data = dict()
+
+    for name in var_names:
+        var_data[name] = np.zeros(len(nums))
 
     for i, num in enumerate(nums):
 
@@ -1389,90 +1438,106 @@ def evolution(path, nums, force=False, interactive=False, pdf=False):
         prop_disk = pickle.load(f)
         f.close()
 
-        time[i] = prop_disk["time"]
+        for name in var_names:
+            try:
+                if name == 'kappa':
+                    var_data[name][i] = - prop_disk['fit']['slope']
+                elif name == 'var':
+                    var_data[name][i] = prop_disk['fit']['var']
+                elif name == 'gamma':
+                    var_data[name][i] = 2 * prop_disk['fit_cs']['slope'] + 1
+                else:
+                    var_data[name][i] = prop_disk[name]
+            except KeyError as e:
+                print('[{}, {}] Key error {}'.format(path, num, e))
+                var_data[name][i] = np.nan
+                
 
-        try:
-            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"]
-        except:
-            for array in [alpha_rey, alpha_grav, Qmin, Qmean, mass_disk, mass_box]:
-                array[i] = np.nan
+    time = var_data['time']
+    for pname in plot_names:
+        
+        if pname == 'alpha':    
+            P.yscale('log')
+            P.plot(time, abs(var_data['alpha_rey_mean']), 'o-.', label=r"$\bar{\alpha}_{Reynolds}$")
+            P.plot(time, abs(var_data['alpha_grav_mean']), '*--', label=r"$\bar{\alpha}_{grav}$")
+            P.legend()
+        elif pname == 'Q':
+            P.plot(time, var_data['Q_min'], 'o-.', label=r"$Q_{min}$")
+            P.plot(time, var_data['Q_mean'], '*--', label=r"$\overline{Q}$")
+            P.legend()
+        elif pname == 'mass':        
+            P.plot(time, var_data['mass_disk'], 'o-.', label=r"$M_{disk}$")
+            P.plot(time, var_data['mass_box'], '*--', label=r"$M_{box}$")
+            P.legend()
+        else:
+            P.plot(time, var_data[pname], '*--')
 
-        if pdf:
-            slope[i] = prop_disk["fit"]["slope"]
-
-    # 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.tight_layout(pad=1)
-        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.tight_layout(pad=1)
-        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.tight_layout(pad=1)
-        P.savefig(path + "/mass_time" + out_ext)
-        P.close()
-
-    # PDF
-    if pdf:
-        P.grid()
-        print(time, slope)
-        P.plot(time, slope, "o-.")
-        P.legend()
-        P.ylabel(r'$d\log\mathcal{P}_{\Sigma} / d\log\Sigma$')
+        if pname =='dcoldens_max':
+            P.plot(time, 1.5 * np.ones(len(time)), 'r-', lw=4)
+        
+        P.ylabel(plot_labels[pname])
         P.xlabel(r'time (code)')
+        P.grid()
         if interactive:
-            P.figure()
+	    P.figure()
         else:
             P.tight_layout(pad=1)
-            P.savefig(path + "/dcolslope_time" + out_ext)
+            P.savefig(path + '/' + pname + 'time_'  + tag + out_ext)
             P.close()
+ 
+def clump_study(path, num, path_out, tag, force=False, rho_thres=1e3, lvl_thres=10):
+    name =  tag + '_' +  str(num).zfill(5)
+    hop_save =  name  + '_hop' + '_prop_struct.save'
 
+    if not (len(glob.glob(path_out + '/' + hop_save))==1) or  force:
+        me.make_clump_hop(path, num, name  + '_hop',  rho_thres, lvl_thres,
+                          [0.5, 0.5, 0.5], 1,  path_out=path_out + '/', path_hop='./',
+                          force=True, gcomp=False)
+        hop_save = me.clump_properties(name  + '_hop', path, num,
+                                       path_out=path_out + '/', gcomp=False)
+
+    f = open(path_out + '/' + hop_save)
+    cp = pickle.load(f)
+    f.close()
+
+
+    ener = cp['grav_vir'] + cp['therm_vir'] * 3./2
+    rad =  np.sqrt((cp['pos_n_max'][:, 0] - 1.)**2 + (cp['pos_n_max'][:, 1] - 1.)**2)
+    mask = (ener < 0) & (rad > 0.07)
+    pos = cp['pos_n_max'][mask]
+    color = np.log10(cp['n_max'][mask])
+
+    print("Load maps files")
+    f = open(path_out + '/maps_disk_' + tag + '_' + str(num).zfill(5) + '.save')
+    m = pickle.load(f)
+    f.close()
+    print("Mpas files loaded")
+
+    P.scatter(pos[:, 0],pos[:, 1], c = color)
+    
+    im = P.imshow(m['coldens_z'],
+                  extent=m['im_extent'],
+                  origin='lower',
+                  norm=mpl.colors.LogNorm(),
+                  vmin=0.01,
+                  vmax=100)
+    P.xlabel(r'$x$')
+    P.ylabel(r'$y$')
+    cbar = P.colorbar(im)
+    cbar.set_label(r"$\Sigma$")
+    P.tight_layout(pad=1)
+    P.title("time = {} (code)".format(m['time']))
+    P.savefig(path_out + '/clumps_' +  tag + '_' +  str(num).zfill(5) + out_ext)
+    P.close()
+
+    # fig, ax = P.subplots(1, 1)
+    # offsets = list(zip(pos[:, 0], 2 - pos[:, 1]))
+    # size = cp['size_iner2'][mask]
+
+    # ax.add_collection(EllipseCollection(widths=size, heights=size, angles=0, units='xy',
+    #                                 facecolors=P.cm.hsv(color),
+    #                                 offsets=offsets, transOffset=ax.transData))
+
+    # P.imshow(np.log10(m['coldens_z']), extent=m['im_extent'], cmap='plasma')
+
+    # P.savefig(path_out + '/clumps_size_' + tag + '_' +  str(num).zfill(5) + out_ext)
diff --git a/pipeline_disk.py b/pipeline_disk.py
index 0d8d457..a2e682f 100644
--- a/pipeline_disk.py
+++ b/pipeline_disk.py
@@ -1,11 +1,13 @@
 # coding: utf-8
 
 import os
+import glob
 from shutil import copy
 import argparse
 import time
 import numpy as np
 import disk_postprocess as dp
+import module_extract as me
 from functools import reduce
 
 storage_in = "/home/nbrucy/simus/"
@@ -13,57 +15,99 @@ storage_out = "/home/nbrucy/visus/"
 
 parser = argparse.ArgumentParser()
 
-parser.add_argument("runs", help="name of runs", nargs="*", default=["015_iso"])
-parser.add_argument("-b", "--begin", help="id of first output", type=int, default=1)
-parser.add_argument("-e", "--end", help="id of last output", type=int, default=100)
-parser.add_argument("-s", "--step", help="step between two output", type=int, default=1)
-parser.add_argument("-p", "--project", help="specify project name", default="disk")
-parser.add_argument(
-    "-fr",
-    "--force_redo",
-    help="redo plots even if the files already exist",
-    action="store_true",
-)
-parser.add_argument(
-    "-w", "--watch", help="wait and watch for missing outputs", action="store_true"
-)
-parser.add_argument("--skip", help="skip failed loadings", action="store_true")
-parser.add_argument(
-    "-wt",
-    "--waiting_time",
-    help="time between to successive try when watching new outputs (in second)",
-    type=int,
-    default=120,
-)
-parser.add_argument(
-    "-af",
-    "--allowed_failures",
-    help="number of allowed failures when waiting",
-    default=30,
-)
-parser.add_argument("-i", "--interactive", help="Interactive mode", action="store_true")
+parser.add_argument("runs", help='name of runs', nargs='*',
+                    default=[])
 
 
-parser.add_argument(
-    "-d", "--disk", help="compute specific disk radial analysis", action="store_true"
-)
-parser.add_argument(
-    "-pd", "--plot_disk", help="plot specific disk radial analysis", action="store_true"
-)
-parser.add_argument("-m", "--maps", help="compute generic maps", action="store_true")
-parser.add_argument("-pm", "--plot_maps", help="plot generic maps", action="store_true")
-parser.add_argument(
-    "-c", "--compare", help="compare different runs", action="store_true"
-)
-parser.add_argument(
-    "-ev",
-    "--evolution",
-    help="plot evolution of quantities over time",
-    action="store_true",
-)
-parser.add_argument(
-    "--pdf", help="plot pdf of fluctuations of column density", action="store_true"
-)
+parser.add_argument("-wo", "--which_outputs",
+                    choices=['all', 'id', 'time'],
+                    help="Select outputs by time range, id range or all of them",
+                    default='all')
+parser.add_argument("-b", "--begin",
+                    help="id of first output",
+                    type=int,
+                    default=1)
+parser.add_argument("-e", "--end",
+                    help="id of last output",
+                    type=int,
+                    default=100)
+parser.add_argument("-s", "--step",
+                    help="step between two output",
+                    type=int,
+                    default=1)
+parser.add_argument("-tb", "--time_begin",
+                    help="time of first output",
+                    type=float,
+                    default=0.)
+parser.add_argument("-te", "--time_end",
+                    help="time of last output",
+                    type=float,
+                    default=6.)
+
+
+parser.add_argument("-p", "--project",
+                    help="specify project name",
+                    default="disk")
+parser.add_argument("-fr", "--force_redo",
+                    help="redo plots even if the files already exist",
+                    action='store_true')
+parser.add_argument("-w", "--watch",
+                    help="wait and watch for missing outputs",
+                    action='store_true')
+parser.add_argument("--skip",
+                    help="skip failed loadings",
+                    action='store_true')
+parser.add_argument("-wt", "--waiting_time",
+                    help="time between to successive try when watching new outputs (in second)",
+                    type=int,
+                    default=120)
+parser.add_argument("-af", "--allowed_failures",
+                    help="number of allowed failures when waiting",
+                    default=30)
+parser.add_argument("-i", "--interactive",
+                    help="Interactive mode",
+                    action='store_true')
+
+
+parser.add_argument("--tag",
+                    help="Add a special tag on output filemanes",
+                    default='')
+
+
+
+parser.add_argument("-d", "--disk",
+                    help="compute specific disk radial analysis",
+                    action='store_true')
+parser.add_argument("-pd", "--plot_disk",
+                    help="plot specific disk radial analysis",
+                    action='store_true')
+parser.add_argument("-m", "--maps",
+                    help="compute generic maps",
+                    action='store_true')
+parser.add_argument("-pm", "--plot_maps",
+                    help="plot generic maps",
+                    action='store_true')
+parser.add_argument("-c", "--compare",
+                    help="compare different runs",
+                    action='store_true')
+parser.add_argument("-ev", "--evolution",
+                    help="plot evolution of quantities over time",
+                    action='store_true')
+parser.add_argument("--pdf",
+                    help="plot pdf of fluctuations of column density",
+                    action='store_true')
+parser.add_argument("--pdf2",
+                    help="plot  patrick's PDF",
+                    action='store_true')
+parser.add_argument("--cpdf",
+                    help="compare pdf",
+                    action='store_true')
+parser.add_argument("--clump",
+                    help="Do clump study",
+                    action='store_true')
+parser.add_argument("--print_outputs",
+                    help="print names of outputs",
+                    action='store_true')
 
 parser.add_argument(
     "--fft", help="use quick and dirty fft rendering", action="store_true"
@@ -87,44 +131,59 @@ parser.add_argument(
 )
 
 
-parser.add_argument(
-    "--nb_bin", help="Number of bins for azimuthal averages", type=int, default=50
-)
-parser.add_argument(
-    "--binning",
-    help="Kind of binning (logarithmic or linear)",
-    choices=["log", "lin"],
-    default="log",
-)
-parser.add_argument(
-    "--rad_ext", help="Value of the highest bin", type=float, default=1.0
-)
-parser.add_argument(
-    "-x", help="x position of the central point", type=float, default=1.0
-)
-parser.add_argument(
-    "-y", help="y position of the central point", type=float, default=1.0
-)
-parser.add_argument(
-    "-z", help="z position of the central point", type=float, default=1.0
-)
+parser.add_argument("--nb_bin",
+                    help="Number of bins for azimuthal averages",
+                    type=int,
+                    default=50)
+parser.add_argument("--pdf_nb_bin",
+                    help="Number of bins for PDF",
+                    type=int,
+                    default=50)
+parser.add_argument("--binning",
+                    help="Kind of binning (logarithmic or linear)",
+                    choices=['log', 'lin'],
+                    default='log')
+parser.add_argument("--rad_ext",
+                    help="Value of the highest bin",
+                    type=float,
+                    default=1.)
+parser.add_argument("-x",
+                    help="x position of the central point",
+                    type=float,
+                    default=1.)
+parser.add_argument("-y",
+                    help="y position of the central point",
+                    type=float,
+                    default=1.)
+parser.add_argument("-z",
+                    help="z position of the central point",
+                    type=float,
+                    default=1.)
 
 parser.add_argument(
     "-ta", "--time_avg", help="Plot time averaged comparaison", action="store_true"
 )
 
 
-parser.add_argument(
-    "--colormap", help="Colormap used", choices=dp.P.colormaps(), default="plasma"
-)
-parser.add_argument(
-    "--format",
-    help="Format of the output",
-    choices=["png", "jpeg", "pdf", "svg", "ps"],
-    default="jpeg",
-)
-parser.add_argument("--dpi", help="Resolution of the output", type=int, default=400)
-parser.add_argument("--beamer", help="Beamer mode", action="store_true")
+parser.add_argument("--colormap",
+                    help="Colormap used",
+                    choices=dp.P.colormaps(),
+                    default='plasma')
+parser.add_argument("--format",
+                    help="Format of the output",
+                    choices=['png', 'jpeg', 'pdf', 'svg', 'ps'],
+                    default='jpeg')
+parser.add_argument("--dpi",
+                    help="Resolution of the output",
+                    type=int,
+                    default=400)
+parser.add_argument("--beamer",
+                    help="Beamer mode",
+                    action='store_true')
+parser.add_argument("--Q_in_name",
+                    help="Whether to use the name to have Q",
+                    action='store_true')
+
 
 
 args = parser.parse_args()
@@ -139,17 +198,25 @@ step = args.step
 rad = 0.5 / args.zoom
 
 # extension for out files
-dp.out_ext = "." + args.format
-if format == "pdf":
-    dp.P.style.use("pdf")
+dp.out_ext = '.' + args.format
+dp.P.style.use("seaborn-deep")
+if args.format == 'pdf':
+     dp.P.style.use("~/.config/matplotlib/pdf.mplstyle")
+
+
 
 if args.beamer:
-    dp.P.rcParams["font.family"] = "sans-serif"
-    dp.P.rcParams["figure.figsize"] = (5, 3.5)
+    dp.P.rcParams['font.family'] = 'sans-serif'
+    dp.P.rcParams['figure.figsize'] = (7, 4.5)
 
 # Plot properties
-dp.P.rcParams["image.cmap"] = args.colormap
-dp.P.rcParams["savefig.dpi"] = args.dpi
+dp.P.rcParams['image.cmap'] = args.colormap
+dp.P.rcParams['savefig.dpi'] = args.dpi
+dp.P.rcParams['lines.linewidth'] = 2
+dp.P.rcParams['lines.markersize'] = 10
+dp.P.rcParams["errorbar.capsize"] = 4
+
+me.P = dp.P
 
 # List of id that were successfully computed
 run_succeded = {}
@@ -159,6 +226,11 @@ for run in runs:
     path_in = storage_in + path_suffix
     path_out = storage_out + path_suffix
 
+    if args.tag == '':
+        tag_run = run
+    else:
+        tag_run = run + '_' + args.tag
+
     if not os.path.exists(path_out):
         os.makedirs(path_out)
         try:
@@ -169,92 +241,95 @@ for run in runs:
 
     run_succeded[run] = []
 
-    for i in range(first, last + 1, step):
+    if args.which_outputs in ['all', 'time'] :
+        names    = glob.glob(path_in + "/output_[0-9][0-9][0-9][0-9][0-9]")
+        nums_all = [int(n.split('/')[-1].split('_')[1]) for n in names]
+        nums_all = np.sort(nums_all)
+        if  args.which_outputs == 'all':
+            nums = nums_all
+        else:
+            time = [dp.get_time(path_in, n) for n in nums_all]
+            nums = [n for i,n in enumerate(nums_all) if time[i] >= args.time_begin and time[i] < args.time_end]
+    else:
+        nums = range(first, last + 1, step)
+
+    for num in nums:
         failures = 0
         success = False
 
         while not success:
             try:
                 maps_disk = None
+                if args.print_outputs:
+                    print("[{}, {}]".format(run, num))
                 if args.maps:
-                    print("[{}, {}] computing maps".format(run, i))
-                    maps_disk = dp.compute_image_data(
-                        path_in,
-                        i,
-                        radius=rad,
-                        path_out=path_out,
-                        tag=run,
-                        map_size=args.mapsize,
-                        force=args.force_redo,
-                        axes_los=args.axes,
-                        images=args.images,
-                        pos_star=np.array([args.x, args.y, args.z]),
-                        fft=args.fft,
-                    )
-                    print("[{}, {}] maps computed".format(run, i))
+                    print("[{}, {}] computing maps".format(run, num))
+                    maps_disk = dp.compute_image_data(path_in, num,
+                                                   radius=rad,
+                                                   path_out=path_out,
+                                                   tag=tag_run,
+                                                   map_size=args.mapsize,
+                                                   force=args.force_redo,
+                                                   axes_los=args.axes,
+                                                   images=args.images,
+                                                   pos_star=np.array([args.x, args.y, args.z]),
+                                                   fft=args.fft)
+                    print("[{}, {}] maps computed".format(run, num))
                 if args.plot_maps:
-                    print("[{}, {}] plotting maps".format(run, i))
-                    maps_disk = dp.plot_maps(
-                        path_out,
-                        i,
-                        maps_disk=maps_disk,
-                        axes_los=args.axes,
-                        images=args.images,
-                        tag=run,
-                        force=args.force_redo,
-                        interactive=args.interactive,
-                        put_title=(not args.beamer),
-                    )
-                    print("[{}, {}] maps plotted".format(run, i))
+                    print("[{}, {}] plotting maps".format(run, num))
+                    maps_disk = dp.plot_maps(path_out, num,
+                                             maps_disk=maps_disk,
+                                             axes_los=args.axes,
+                                             images=args.images,
+                                             tag=tag_run,
+                                             force=args.force_redo,
+                                             interactive=args.interactive,
+                                             put_title=(not args.beamer))
+                    print("[{}, {}] maps plotted".format(run, num))
                 if args.disk:
-                    print("[{}, {}] computing disk props".format(run, i))
-                    dp.disk_prop(
-                        path_in,
-                        i,
-                        path_out=path_out,
-                        nb_bin=args.nb_bin,
-                        binning=args.binning,
-                        rad_ext=args.rad_ext,
-                        force=args.force_redo,
-                        pos_star=np.array([args.x, args.y, args.z]),
-                    )
-                    print("[{}, {}] disk_props computed".format(run, i))
+                    print("[{}, {}] computing disk props".format(run, num))
+                    dp.disk_prop(path_in, num, path_out=path_out,
+                                 nb_bin=args.nb_bin,
+                                 binning=args.binning,
+                                 rad_ext=args.rad_ext,
+                                 force=args.force_redo,
+                                 pos_star=np.array([args.x, args.y, args.z]))
+                    print("[{}, {}] disk_props computed".format(run, num))
                 if args.plot_disk:
-                    print("[{}, {}] plotting disk props".format(run, i))
-                    dp.plot_disk_prop(
-                        path_out,
-                        i,
-                        tag=run,
-                        force=args.force_redo,
-                        interactive=args.interactive,
-                        put_title=(not args.beamer),
-                    )
-                    print("[{}, {}] disk props plotted".format(run, i))
+                    print("[{}, {}] plotting disk props".format(run, num))
+                    dp.plot_disk_prop(path_out, num, tag=tag_run,
+                                      force=args.force_redo,
+                                      interactive=args.interactive,
+                                      put_title=(not args.beamer))
+                    print("[{}, {}] disk props plotted".format(run, num))
                 if args.pdf:
-                    print("[{}, {}] computing pdf".format(run, i))
-                    dp.disk_pdf(
-                        path_out,
-                        i,
-                        maps_disk,
-                        pos_star=np.array([args.x, args.y, args.z]),
-                        force=args.force_redo,
-                        tag=run,
-                        interactive=args.interactive,
-                        put_title=(not args.beamer),
-                    )
-                    print("[{}, {}] pdf computed".format(run, i))
+                    print("[{}, {}] computing pdf #1".format(run, num))
+                    dp.disk_pdf(path_out, num, maps_disk,
+                                pos_star=np.array([args.x, args.y, args.z]),
+                                force=args.force_redo,
+                                tag=tag_run,
+                                nb_bin_hist=args.pdf_nb_bin,
+                                interactive=args.interactive,
+                                put_title=(not args.beamer))
+                    print("[{}, {}] pdf #1 computed".format(run, num))
+                if args.pdf2:
+                    print("[{}, {}] computing pdf #2".format(run, num))
+                    me.get_pdf(path_in, num, path_out=path_out,
+                               force=args.force_redo,
+                               tag=tag_run)
+                    print("[{}, {}] pdf #2 computed".format(run, num))
+                if args.clump:
+                    print("[{}, {}] computing clumps".format(run, num))
+                    dp.clump_study(path_in, num, path_out, tag_run)
+                    print("[{}, {}] clumps computed".format(run, num))
                 # If we are here, success !
                 success = True
-                run_succeded[run].append(i)
+                run_succeded[run].append(num)
             except (ValueError, IOError, KeyError) as e:
                 print(e)
                 if args.watch and failures < args.allowed_failures:
                     failures = failures + 1
-                    print(
-                        "Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(
-                            run, i, args.waiting_time, args.allowed_failures - failures
-                        )
-                    )
+                    print("Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(run, num, args.waiting_time, args.allowed_failures - failures))
                     time.sleep(args.waiting_time)
                 elif args.skip:
                     break
@@ -262,34 +337,25 @@ for run in runs:
                     raise
     if args.evolution:
         print("[{}] plotting evolution".format(run))
-        dp.evolution(
-            path_out,
-            run_succeded[run],
-            force=args.force_redo,
-            interactive=args.interactive,
-            pdf=args.pdf,
-        )
+        dp.evolution(path_out, run_succeded[run],
+                     force=args.force_redo,
+                     tag=args.tag,
+                     interactive=args.interactive,
+                     pdf=args.pdf or args.cpdf)
         print("[{}] evolution plotted".format(run))
 
 if args.compare:
     path_suffix = project
     path = storage_out + path_suffix
 
-    if args.time_avg:
-        # Select output availables for at least on run
-        output_ok = reduce(np.union1d, [run_succeded[run] for run in runs]).astype(int)
-        print(output_ok)
-        if len(output_ok) >= 1:
-            dp.compare(
-                path,
-                runs,
-                output_ok,
-                path_out=path + "/comp",
-                force=args.force_redo,
-                interactive=args.interactive,
-                Q_in_name=(not args.pdf),
-                pdf=args.pdf,
-            )
+    if (args.time_avg):
+        dp.compare(path, runs, run_succeded,
+                   path_out=path + '/comp',
+                   force=args.force_redo,
+                   interactive=args.interactive,
+                   Q_in_name=args.Q_in_name,
+                   tag=args.tag,
+                   pdf=args.pdf or args.cpdf)
     else:
         for i in range(first, last + 1, step):
             # Select usable runs
@@ -297,16 +363,13 @@ if args.compare:
 
             try:
                 if len(runs_ok) > 1:
-                    dp.compare(
-                        path,
-                        runs_ok,
-                        i,
-                        path_out=path + "/comp",
-                        force=args.force_redo,
-                        interactive=args.interactive,
-                        Q_in_name=(not args.pdf),
-                        pdf=args.pdf,
-                    )
+                    dp.compare(path, runs_ok, i,
+                               path_out=path + '/comp',
+                               force=args.force_redo,
+                               interactive=args.interactive,
+                               Q_in_name=args.Q_in_name,
+                               tag=args.tag,
+                               pdf=args.pdf or args.cpdf)
             except (KeyError, IOError) as e:
                 print(e)
                 if args.skip: