diff --git a/fragdisk_galactica.py b/fragdisk_galactica.py
index f5dbc2d..a548cb4 100644
--- a/fragdisk_galactica.py
+++ b/fragdisk_galactica.py
@@ -28,15 +28,17 @@
 from __future__ import print_function, unicode_literals
 
 import os
-
+import argparse
 import h5py
 import numpy as np
+import matplotlib as mpl
 from astrophysix import units as U
 from astrophysix.simdm import Project, ProjectCategory, SimulationStudy
 from astrophysix.simdm.datafiles import PlotInfo, PlotType
 from astrophysix.utils.file import FileType
 
-from plotter import Plotter, default_params
+from plotter import Plotter
+from pp_params import default_params
 from ramses_astrophysix import ramses
 
 pp_params = default_params()
@@ -53,6 +55,9 @@ pp_params.pdf.nb_bin = 100
 pp_params.process.num_process = 10
 in_dir = "/drf/projets/alfven-data/nbrucy/simus/fragdisk"
 out_dir = "/dsm/anais/storageA/nbrucy/visus/fragdisk/mnras"
+in_dir_conv = "/drf/projets/alfven-data/nbrucy/simus/conv_disk"
+out_dir_conv = "/dsm/anais/storageA/nbrucy/visus/conv_disk"
+
 nml_key = "cloud_params/beta_cool"
 
 
@@ -177,32 +182,36 @@ pl_units = Plotter(
 info = pl_units.comp.info
 rd = U.Unit.create_unit(
     "r_d",
+    latex="r_\\mathrm{d}",
     base_unit=0.25 * info["unit_length"] / info["boxlen"],
 )
 rhod = U.Unit.create_unit(
     "rho_d",
-    latex="\rho_d",
+    latex="\\rho_\\mathrm{d}",
     base_unit=6.36 * info["unit_density"],
 )
 rho_u = U.Unit.create_unit(
-    "M_*/r_d^3",
-    latex="M_\\star.r_d^{-3}",
+    "rho_u",
+    latex="\mathrm{M}_\\star.r_\\mathrm{d}^{-3}",
     base_unit=info["unit_mass"] * rd ** (-3),
 )
+orp = U.Unit.create_unit(
+    "ORP",
+    base_unit=2 * np.pi * np.sqrt(0.25 ** 3) * info["unit_time"],
+)
+vkd = U.Unit.create_unit(
+    "vkd",
+    latex="v_\\mathrm{k,d}",
+    base_unit=2 * np.pi * rd / orp,
+)
 Pd = U.Unit.create_unit(
     "P_d",
     base_unit=(0.2 * info["unit_velocity"]) ** 2 * rhod,
 )
-
-orp = U.Unit.create_unit(
-    "ORP",
-    base_unit=2 * np.pi * 0.25 ** -3 * info["unit_time"],
-)
-vkd = rd * orp
 P_u = U.Unit.create_unit(
-    "(k_kd * M_*)/r_d^3",
-    latex="v_{k,d}^2.M_\\star.r_d^{-3}",
-    base_unit=vkd * rho_u,
+    "P_u",
+    latex="v_\\mathrm{k,d}^2.\\mathrm{M}_\\star.r_\\mathrm{d}^{-3}",
+    base_unit=vkd ** 2 * rho_u,
 )
 
 Sigmad = U.Unit.create_unit(
@@ -211,7 +220,7 @@ Sigmad = U.Unit.create_unit(
     base_unit=0.25 * info["unit_density"] * info["unit_length"] / info["boxlen"],
 )
 Sigma_u = U.Unit.create_unit(
-    "M*/r_d^2",
+    "Sigma_u",
     latex="M_\\star.r_d^{-2}",
     base_unit=rd ** (-2) * info["unit_mass"],
 )
@@ -219,82 +228,95 @@ Sigma_u = U.Unit.create_unit(
 
 # ---- Runs  -----
 
-# JR13_TIC
+groups = ["jr11", "jr12", "jr12_tic", "jr13_tic"]
+# groups = ["jr13_tic"]
+pls = []
 
-pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
-pp_params.astrophysix.descr_fmt = """
-    <p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
-    """
+for group in groups:
+    if group == "jr13_tic":
+        # JR13_TIC
 
-runs = "*_jr13"
-pl_jr13 = Plotter(
-    in_dir,
-    filter_name=runs,
-    in_nums="all",
-    sort_run_by=nml_key,
-    path_out=out_dir,
-    tag="jr13_tic_mnras",
-    pp_params=pp_params,
-    unit_time=orp,
-)
-print("JR13_TIC defined")
+        pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
+        pp_params.astrophysix.descr_fmt = """
+        <p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
+        """
 
-# JR12_TIC
-runs_12 = "0[0-9][0-9]_beta*_jr12"
-pl_jr12_tic = Plotter(
-    in_dir,
-    filter_name=runs_12,
-    in_nums="all",
-    sort_run_by=nml_key,
-    filter_nml=("cloud_params", "!=", 7),
-    path_out=out_dir,
-    tag="jr12_tic_mnras",
-    pp_params=pp_params,
-    unit_time=orp,
-)
+        runs = "*_jr13"
+        pl = Plotter(
+            in_dir,
+            filter_name=runs,
+            in_nums="all",
+            sort_run_by=nml_key,
+            path_out=out_dir,
+            tag="jr13_tic_mnras",
+            pp_params=pp_params,
+            unit_time=orp,
+        )
 
-pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
-pp_params.astrophysix.descr_fmt = """
-    <p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
-    """
-print("JR12_TIC defined")
+    if group == "jr12_tic":
+        # JR12_TIC
+        pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}"
+        pp_params.astrophysix.descr_fmt = """
+        <p>Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.</p>
+        """
+        runs_12 = "0[0-9][0-9]_beta*_jr12"
+        pl = Plotter(
+            in_dir,
+            filter_name=runs_12,
+            in_nums="all",
+            sort_run_by=nml_key,
+            filter_nml=("cloud_params", "!=", 7),
+            path_out=out_dir,
+            tag="jr12_tic_mnras",
+            pp_params=pp_params,
+            unit_time=orp,
+        )
 
-# JR12
-in_dir_conv = "/drf/projets/alfven-data/nbrucy/simus/conv_disk"
-out_dir_conv = "/dsm/anais/storageA/nbrucy/visus/conv_disk"
-runs = "[7-8][0-9]_beta*_j*"
-pl_jr12 = Plotter(
-    in_dir_conv,
-    filter_name=runs,
-    in_nums="all",
-    sort_run_by=nml_key,
-    path_out=out_dir,
-    tag="jr12_mnras",
-    pp_params=pp_params,
-    unit_time=orp,
-)
-print("JR12 defined")
+    if group == "jr12":
 
-# JR11
-runs = "*beta*_jr11"
-pl_l11 = Plotter(
-    in_dir_conv,
-    filter_name=runs,
-    sort_run_by=nml_key,
-    filter_nml=("cloud_params/beta_cool", ">", 3),
-    path_out=out_dir_conv,
-    tag="jr11_mnras",
-    pp_params=pp_params,
-    unit_time=orp,
-)
+        pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
+        pp_params.astrophysix.descr_fmt = """
+        <p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
+        """
 
-print("JR11 defined")
+        # JR12
+        runs = "[7-8][0-9]_beta*_j*"
+        pl = Plotter(
+            in_dir_conv,
+            filter_name=runs,
+            in_nums="all",
+            sort_run_by=nml_key,
+            path_out=out_dir,
+            tag="jr12_mnras",
+            pp_params=pp_params,
+            unit_time=orp,
+        )
 
+    if group == "jr11":
 
-pls = [pl_l11, pl_jr12, pl_jr12_tic, pl_jr13]
+        pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}"
+        pp_params.astrophysix.descr_fmt = """
+        <p>Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}</p>
+        """
 
+        # JR11
+        runs = "*beta*_jr11"
+        pl = Plotter(
+            in_dir_conv,
+            filter_name=runs,
+            sort_run_by=nml_key,
+            filter_nml=("cloud_params/beta_cool", ">", 3),
+            path_out=out_dir_conv,
+            tag="jr11_mnras",
+            pp_params=pp_params,
+            unit_time=orp,
+        )
+
+    pls.append(pl)
+    print("{} defined".format(group.upper()))
 
 # -------------------------------------------------------------------------------------------------------------------- #
+
 redo = True
 for pl in pls:
     pl.pp_params.process.verbose = True
@@ -307,13 +329,17 @@ for pl in pls:
 # -------------------------------------------------------------------------------------------------------------------- #
 
 
+select = {
+    "time": 4.5,
+    #    "filter_nml" : (nml_key, "=", 8),
+}
+
 for pl in pls:
-    select = {"time": 4.5}
-
     # Edit descriptions
-
     pl.rules["slice_rho"].description = "Density slice"
     pl.rules["slice_P"].description = "Pressure slice"
+    pl.rules["slice_velr"].description = "Radial velocity slice"
+    pl.rules["slice_velphi"].description = "Orthoradial velocity slice"
     pl.rules[
         "pdf_coldens"
     ].description = """
@@ -321,42 +347,49 @@ for pl in pls:
     $\\sigma = \\Sigma/\\overline{\\Sigma}$ with respect to its azimuthal average
     """
 
-    coldens_kwargs = {
+    # define plot parameters
+    map_kwargs = {
         "overwrite": redo,
         "overwrite_dep": False,
         "unit_space": rd,
         "center_space": True,
         "unit_time": orp,
-        "unit": Sigma_u,
         "nml_key": "cloud_params/beta_cool",
         "select": select,
-        "label": r"$\Sigma$",
+    }
+
+    coldens_kwargs = {
+        **map_kwargs,
+        **{
+            "unit": Sigma_u,
+            "label": r"$\Sigma$",
+        },
     }
 
     rho_kwargs = {
-        "overwrite": redo,
-        "overwrite_dep": False,
-        "unit_space": rd,
-        "center_space": True,
-        "unit_time": orp,
-        "unit": rho_u,
-        "nml_key": "cloud_params/beta_cool",
-        "select": select,
-        "label": r"$\rho$",
+        **map_kwargs,
+        **{
+            "unit": rho_u,
+            "label": r"$\rho$",
+        },
     }
 
     P_kwargs = {
-        "overwrite": redo,
-        "overwrite_dep": False,
-        "unit_space": rd,
-        "center_space": True,
-        "unit_time": orp,
-        "unit": P_u,
-        "nml_key": "cloud_params/beta_cool",
-        "select": select,
-        "label": r"$P$",
+        **map_kwargs,
+        **{
+            "unit": P_u,
+            "label": r"$P$",
+        },
     }
 
+    vel_kwargs = {
+        **map_kwargs,
+        **{
+            "unit": vkd,
+        },
+    }
+
+    # do plots
     pl.coldens("z", **coldens_kwargs)
     pl.coldens("y", **coldens_kwargs)
 
@@ -365,6 +398,15 @@ for pl in pls:
 
     pl.slice_P("z", **P_kwargs)
 
+    pl.slice_velr(
+        "z",
+        label=r"$v_r$",
+        norm=mpl.colors.SymLogNorm(0.1, vmin=-2, vmax=2, base=10),
+        autoscale=False,
+        **vel_kwargs,
+    )
+    pl.slice_velphi("z", label=r"$v_\varphi$", **vel_kwargs)
+
     pl.pdf_coldens(
         "z",
         overwrite=redo,
@@ -400,15 +442,17 @@ for simu in proj.simulations:
         # Change unit to things Galactica understands
         snap.time = (snap.time[0], U.none)
         snap.description = """
-        <p>Snapshot at {time} ORPs.</p>
+        <p>Snapshot at {time:.3g} ORPs.</p>
+
+        <h5 class="sn_panel">Notations:</h5>
 
-        <h5> Notations </h5>
         <p>
         $r_d$ : radius of the disk <br />
         ORP : outer rotation period, that is period of the gas at $r_d$ <br />
         $M_\\star$ : Mass of the central object <br />
         $v_{kepl}$ : Keplerian speed at $r_d$
         </p>
+        <p>All slices are taken at $z=0$ or $y=0$.
         """.format(
             time=snap.time[0], kepl="{k,d}"
         )