diff --git a/disk_postprocess.py b/disk_postprocess.py
index 3ef29df..07a801e 100644
--- a/disk_postprocess.py
+++ b/disk_postprocess.py
@@ -825,7 +825,7 @@ def disk_pdf(
     interactive=False,
     nb_bin_hist=50,
     tag="",
-    rad_min=0.05,
+    rad_min=0.075,
     rad_max=0.3,
     put_title=True,
     do_speed=True,
@@ -863,10 +863,14 @@ def disk_pdf(
         f.close()
     print("maps file loaded")
 
+    # Properties
     time = prop_disk["time"]
     im_extent = maps_disk["im_extent"]
     title = tag.split("_")[1] + " t=" + str(time)[0:5] + " (code)"
 
+    # Load coldens
+    coldens_map = maps_disk["coldens_z"]
+
     # radius of the corner of the box plus a margin
     rad_box = (
         np.sqrt((im_extent[1] - pos_star[0]) ** 2 + (im_extent[3] - pos_star[1]) ** 2)
@@ -877,13 +881,23 @@ def disk_pdf(
     rad_bins = 0.5 * (rad_bins[0:-1] + rad_bins[1:])
     rad_bins = np.concatenate(([0.0], rad_bins, [rad_box]))
 
-    coldens_map = maps_disk["coldens_z"]
-
     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)
+    bins_flat = bins.flatten()
+    rr_flat = rr.flatten()
+
+    # Mask selecting the zone of interest
+    mask_map = (rr > rad_min) & (rr < rad_max)
+    mask_flat = mask_map.flatten()
+
+    # Additionnal maps
     rho_map = maps_disk["rho_z"]
     cs_map = np.sqrt(maps_disk["T_z"])
 
@@ -896,49 +910,45 @@ def disk_pdf(
     vrad_map = (vx_map * xx_star + vy_map * yy_star) / rr
     vaz_map = (vy_map * xx_star - vx_map * yy_star) / rr
 
-    maps = [coldens_map, rho_map, cs_map, v_map, vaz_map]
+    maps = {
+        "coldens": coldens_map,
+        "rho": rho_map,
+        "cs": cs_map,
+        "v": v_map,
+        "vaz": vaz_map,
+    }
 
-    # 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)
-    bins_flat = bins.flatten()
-    rr_flat = rr.flatten()
+    avg_maps = {}
+    fluct_maps = {}
 
-    coldens_rad = np.zeros(len(rad_bins) - 1)
-    rho_rad = np.zeros(len(rad_bins) - 1)
-    cs_rad = np.zeros(len(rad_bins) - 1)
-    v_rad = np.zeros(len(rad_bins) - 1)
-    vaz_rad = np.zeros(len(rad_bins) - 1)
-
-    rads = [coldens_rad, rho_rad, cs_rad, v_rad, vaz_rad]
-    means = [None] * len(rads)
-
-    for i, rad_arr in enumerate(rads):
-        map_arr = maps[i]
+    for cur_map in maps:
+        map_arr = maps[cur_map]
 
+        # mean of all the cells in the bin
+        mean_bin = np.zeros(len(rad_bins) - 1)
         for j in range(len(rad_bins) - 1):
-            rad_arr[j] = np.mean(map_arr[bins == j])
+            mean_bin[j] = np.mean(map_arr[bins == j])
 
         # Add value for borders
-        rad_arr = np.concatenate(([rad_arr[0]], rad_arr))
+        mean_bin = np.concatenate(([mean_bin[0]], mean_bin))
 
-        # Retrieve the mean of the bins
+        # Compute the map azimuthally averaged
         # use linear interpolation to improve accuracy
-        means[i] = (rad_bins[bins_flat + 1] - rr_flat) * rad_arr[bins_flat]
-        means[i] = means[i] + (rr_flat - rad_bins[bins_flat]) * rad_arr[bins_flat + 1]
-        means[i] = means[i] / (rad_bins[bins_flat + 1] - rad_bins[bins_flat])
+        avg_flat = (rad_bins[bins_flat + 1] - rr_flat) * mean_bin[bins_flat]
+        avg_flat = avg_flat + (rr_flat - rad_bins[bins_flat]) * mean_bin[bins_flat + 1]
+        avg_flat = avg_flat / (rad_bins[bins_flat + 1] - rad_bins[bins_flat])
+        avg_maps[cur_map] = np.reshape(avg_flat, rr.shape)
 
-    coldens_mean = means[0]
-    coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape)
-    #  Kill fluctuations for r < rad_min and r > rad_max
-    mask_map = (rr > rad_min) & (rr < rad_max)
-    coldens_map[np.logical_not(mask_map)] = np.nan
+        #  Select zone of interest
+        avg_maps[cur_map][np.logical_not(mask_map)] = np.nan
+
+        # Compute fluctuation
+        fluct_maps[cur_map] = map_arr / avg_maps[cur_map]
 
     # Histogramm : density fluctuation distribution function
-    dcoldens = np.log10(coldens_map.flatten() / coldens_mean)
-    mask_hist = (rr_flat > rad_min) & (rr_flat < rad_max)
-    nb_cells = np.sum(mask_hist)
+    dcoldens = np.log10(fluct_maps["coldens"]).flatten()
+
+    nb_cells = np.sum(mask_flat)
     P.grid()
     P.yscale("log")
     P.ylim([0.5 / nb_cells, 1.0])
@@ -947,7 +957,7 @@ def disk_pdf(
     if put_title:
         P.title(title)
     values, edges, _ = P.hist(
-        dcoldens[mask_hist],
+        dcoldens[mask_flat],
         nb_bin_hist,
         range=(-1, 3),
         weights=np.ones(nb_cells) / nb_cells,
@@ -955,7 +965,7 @@ def disk_pdf(
     centers = 0.5 * (edges[1:] + edges[:-1])
 
     # Variance
-    var = np.var(dcoldens[mask_hist])
+    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):
@@ -988,26 +998,27 @@ def disk_pdf(
         P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext)
         P.close()
 
+    # Derived quantities
+    drho = fluct_maps["rho"].flatten()
+    dcs = fluct_maps["cs"].flatten()
+    dv = fluct_maps["v"].flatten()
+    dvaz_kepl = abs(maps["vaz"] - v_kepl) / v_kepl
+    fluct_maps["vaz_kepl"] = dvaz_kepl
+    dmach = abs(maps["v"] - avg_maps["v"]) / maps["cs"]
+    fluct_maps["mach"] = dmach
+    dmach_mean = np.mean(dmach[mask_map].flatten())
+
     # Fluctuations plots
-    rho_mean = means[1]
-    cs_mean = means[2]
-    v_mean = means[3]
-    vaz_mean = means[4]
-    drho = abs(rho_map.flatten() / rho_mean)
-    dcs = abs(cs_map.flatten() / cs_mean)
-    dvaz = abs(vaz_map.flatten() / vaz_mean)
-    dv = abs(v_map.flatten() / v_mean)
-    dmach = abs((v_map.flatten() - v_mean) / cs_map.flatten())
-    dmach_mean = np.mean(dmach[mask_hist])
     f = open(name_prop, "w")
     prop_disk["dmach_mean"] = dmach_mean
+    print("dmach_mean = {}".format(dmach_mean))
     pickle.dump(prop_disk, f)
     f.close()
 
-    # P.scatter(drho[mask_hist], dcs[mask_hist], s=0.1)
+    # dcs = f(drho)
     P.hist2d(
-        np.log10(drho[mask_hist]),
-        np.log10(dcs[mask_hist]),
+        np.log10(drho[mask_flat]),
+        np.log10(dcs[mask_flat]),
         (1000, 1000),
         norm=mpl.colors.LogNorm(),
     )
@@ -1023,9 +1034,10 @@ def disk_pdf(
         P.savefig(path + "/dcs_" + tag + "_" + str(num).zfill(5) + out_ext)
         P.close()
 
+    # dv = f(drho)
     P.hist2d(
-        np.log10(drho[mask_hist]),
-        dvaz[mask_hist],
+        np.log10(drho[mask_flat]),
+        dv[mask_flat],
         (1000, 1000),
         norm=mpl.colors.LogNorm(),
     )
@@ -1037,160 +1049,67 @@ def disk_pdf(
         P.figure()
     else:
         P.tight_layout(pad=1)
-        P.savefig(path + "/dvaz_" + tag + "_" + str(num).zfill(5) + out_ext)
+        P.savefig(path + "/dv_" + tag + "_" + str(num).zfill(5) + out_ext)
         P.close()
 
     # Fluctuations maps
-    dcoldens_map = np.reshape(dcoldens, coldens_map.shape)
 
-    im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis")
-
-    if put_title:
-        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_" + tag + "_" + format(num, "05")
-    name_im = name + out_ext
-
-    if interactive:
-        P.figure()
-    else:
-        P.tight_layout(pad=1)
-        P.savefig(name_im)
-        P.close()
-
-    dcs_map = np.reshape(dcs, coldens_map.shape)
-    dcs_map[np.logical_not(mask_map)] = np.nan
-
-    im = P.imshow(
-        np.log10(dcs_map),
-        extent=im_extent,
-        origin="lower",
-        cmap="RdBu_r",
-        vmin=-0.6,
-        vmax=0.6,
-    )
-
-    if put_title:
-        P.title(title)
-    P.xlabel(r"$x$")
-    P.ylabel(r"$y$")
-
-    cbar = P.colorbar(im)
-    cbar.set_label(r"$log(c_s/\bar{c_s})$ (code)")
-    name = path + "/dcsmap_" + tag + "_" + format(num, "05")
-    name_im = name + out_ext
-
-    if interactive:
-        P.figure()
-    else:
-        P.tight_layout(pad=1)
-        P.savefig(name_im)
-        P.close()
-
-    drho_map = np.reshape(drho, coldens_map.shape)
-    drho_map[np.logical_not(mask_map)] = np.nan
-
-    im = P.imshow(
-        np.log10(drho_map),
-        extent=im_extent,
-        origin="lower",
-        cmap="RdBu_r",
-        vmin=-2.0,
-        vmax=2.0,
-    )
-
-    if put_title:
-        P.title(title)
-    P.xlabel(r"$x$")
-    P.ylabel(r"$y$")
-
-    cbar = P.colorbar(im)
-    cbar.set_label(r"$log(\rho/\bar{\rho})$ (code)")
-    name = path + "/drhomap_" + tag + "_" + format(num, "05")
-    name_im = name + out_ext
-
-    if interactive:
-        P.figure()
-    else:
-        P.tight_layout(pad=1)
-        P.savefig(name_im)
-        P.close()
-
-    if do_speed:
-        dvel_kepl = abs(v_map - v_kepl) / v_kepl
-        dvaz_kepl = abs(vaz_map - v_kepl) / v_kepl
-        dvrad_vaz = abs(vrad_map / vaz_map)
-        dvaz_map = np.reshape(dvaz, coldens_map.shape)
-        dvaz_map[np.logical_not(mask_map)] = np.nan
-
-        im = P.imshow(
-            dvaz_kepl,
-            extent=im_extent,
-            origin="lower",
-            cmap="RdBu_r",
-            norm=mpl.colors.LogNorm(),
-            vmax=1.0,
-            vmin=0.01,
-        )
-        cbar = P.colorbar(im)
-        cbar.set_label(r"$|v_\varphi - v_{kepl}|/v_{kepl}$")
-        P.xlabel(r"$x$")
-        P.ylabel(r"$y$")
+    for cur_map in fluct_maps:
+        fluct_map = fluct_maps[cur_map]
+        if cur_map == "coldens":
+            im = P.imshow(
+                fluct_map,
+                norm=mpl.colors.LogNorm(),
+                extent=im_extent,
+                origin="lower",
+                cmap="viridis",
+            )
+            label = r"$log(N/\bar{N})$"
+        elif cur_map == "cs":
+            im = P.imshow(
+                np.log10(fluct_map),
+                extent=im_extent,
+                origin="lower",
+                cmap="RdBu_r",
+                vmin=-0.6,
+                vmax=0.6,
+            )
+            label = r"$log(c_s/\bar{c_s})$"
+        elif cur_map == "rho":
+            im = P.imshow(
+                np.log10(fluct_map),
+                extent=im_extent,
+                origin="lower",
+                cmap="RdBu_r",
+                vmin=-2.0,
+                vmax=2.0,
+            )
+            label = r"$log(\rho/\bar{\rho})$"
+        elif cur_map == "vaz_kepl":
+            im = P.imshow(
+                fluct_map,
+                extent=im_extent,
+                origin="lower",
+                cmap="RdBu_r",
+                norm=mpl.colors.LogNorm(),
+                vmax=1.0,
+                vmin=0.01,
+            )
+            label = r"$|v_\varphi - v_{kepl}|/v_{kepl}$"
+        elif cur_map == "vaz":
+            im = P.imshow(fluct_map, extent=im_extent, origin="lower", cmap="RdBu_r")
+            label = r"$v_\varphi / \bar{v_\varphi}$"
+        elif cur_map == "mach":
+            im = P.imshow(fluct_map, extent=im_extent, origin="lower", cmap="RdBu_r")
+            label = r"$|v - \bar{v}| / c_s$"
 
         if put_title:
             P.title(title)
-
-        name = path + "/dvaz_kepl_" + tag + "_" + format(num, "05")
-        name_im = name + out_ext
-
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(name_im)
-            P.close()
-
-        im = P.imshow(
-            vrad_map,
-            extent=im_extent,
-            origin="lower",
-            cmap="RdBu_r",
-            norm=mpl.colors.SymLogNorm(linthresh=0.03),
-            vmax=8.0,
-            vmin=-8,
-        )
-        cbar = P.colorbar(im)
-        cbar.set_label(r"$v_r$")
         P.xlabel(r"$x$")
         P.ylabel(r"$y$")
-
-        if put_title:
-            P.title(title)
-
-        name = path + "/vradmap_" + tag + "_" + format(num, "05")
-        name_im = name + out_ext
-
-        if interactive:
-            P.figure()
-        else:
-            P.tight_layout(pad=1)
-            P.savefig(name_im)
-            P.close()
-
-        im = P.imshow(dvaz_map, extent=im_extent, origin="lower", cmap="RdBu_r")
         cbar = P.colorbar(im)
-        cbar.set_label(r"$v_\varphi / \bar{v_\varphi}$")
-        P.xlabel(r"$x$")
-        P.ylabel(r"$y$")
-
-        if put_title:
-            P.title(title)
-
-        name = path + "/dvazmap_" + tag + "_" + format(num, "05")
+        cbar.set_label(label)
+        name = path + "/d" + cur_map + "map_" + tag + "_" + format(num, "05")
         name_im = name + out_ext
 
         if interactive:
@@ -1242,6 +1161,9 @@ def compare(
         dmach = np.zeros(len(runs))
         dmach_std = np.zeros(len(runs))
 
+        all_var = []
+        all_dmach = []
+
     for i, run in enumerate(runs):
         path_run = path + "/" + run
         nb_outputs = 0
@@ -1276,6 +1198,9 @@ def compare(
                     var_run.append(fit["var"])
                     dmach_run.append(prop_disk["dmach_mean"])
 
+                all_var = all_var + var_run
+                all_dmach = all_dmach + dmach_run
+
                 nb_outputs = nb_outputs + 1
                 print(run, num, nb_outputs)
 
@@ -1387,17 +1312,18 @@ def compare(
             P.savefig(path_out + "/varbeta_" + nums_name + out_ext)
             P.close()
 
-        # dmach = f(var)
+        # var = f(log(dmach)
         P.grid()
-        P.errorbar(var, dmach, xerr=var_std, yerr=dmach_std, fmt="o")
+        P.plot(all_dmach, all_var, "o")
 
         P.legend()
-        P.ylabel(r"$<(v - \bar{v}) / c_s>$")
-        P.xlabel(r"$Var(\log(N / \bar(N))$")
+        P.xlabel(r"$<(v - \bar{v}) / c_s>$")
+        P.ylabel(r"$Var(\log(N / \bar(N))$")
+        P.yscale("log")
 
-        (a, b, rho, _, stderr) = linregress(var, dmach)
-        P.plot(var, a * var + b, "--", linewidth=1.5)
-        print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2))
+        #  (a, b, rho, _, stderr) = linregress(var, log(dmach)
+        #  P.plot(var, a*var + b, '--', linewidth=1.5)
+        #  print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
 
         if interactive:
             P.figure()
diff --git a/pipeline_disk.py b/pipeline_disk.py
index f26bf1b..b231186 100644
--- a/pipeline_disk.py
+++ b/pipeline_disk.py
@@ -72,7 +72,7 @@ parser.add_argument(
     "--images",
     nargs="*",
     default=["coldens", "rho", "speed", "Q", "T"],
-    choices=["coldens", "rho", "speed", "Q", "T", "level", "cpu"],
+    choices=["coldens", "rho", "speed", "Q", "T", "levels", "cpu"],
 )
 parser.add_argument(
     "--axes", nargs="*", default=["x", "y", "z"], choices=["x", "y", "z"]