diff --git a/.gitignore b/.gitignore
index c240fe6..3fc7b58 100644
--- a/.gitignore
+++ b/.gitignore
@@ -17,3 +17,11 @@ scgenerator_log*
 sc-*.log
 
 .vscode
+
+
+# latex
+*.aux
+*.fdb_latexmk
+*.fls
+*.log
+*.synctex.gz
diff --git a/formulas/main.tex b/formulas/main.tex
new file mode 100644
index 0000000..c8e3f1f
--- /dev/null
+++ b/formulas/main.tex
@@ -0,0 +1,34 @@
+\documentclass[preview]{standalone}
+\usepackage{siunitx}
+\usepackage{tikz}
+\usepackage[siunitx]{circuitikz}
+\usepackage{babel}
+\usepackage{acronym}
+\newcommand{\At}{A(z, t)}
+\newcommand{\Et}{E(z, t)}
+\newcommand{\Ew}{\tilde{E}(z, \omega)}
+\newcommand{\PNL}{\tilde{\mathrm{P}}^\mathrm{NL}}
+
+\begin{document}
+Complex envelope such that $|A|^2$ is in W with E in V/m~:
+\begin{equation}
+    \At = \sqrt{\frac12 \epsilon_0 c n A_\mathrm{eff}} \Et
+\end{equation}
+
+Equations with E in V/m
+\begin{equation}
+    \frac{\partial \Ew}{\partial z} = i\left(\beta(\omega) - \frac{\omega}{\nu} \right)\Ew + i \frac{\omega^2}{2c\epsilon_0 \beta(\omega)}\PNL(z, \omega)
+\end{equation}
+
+\begin{equation}
+    \tau = t - \frac{z}{\nu}
+    \,, \quad \xi = z
+    \,, \quad \frac{\partial}{\partial z} = -\frac1{\nu}\frac{\partial}{\partial \tau} + \frac{\partial}{\partial \xi} = -\frac{i \omega}{\nu} + \frac{\partial}{\partial \xi}
+    \,, \quad \frac{\partial}{\partial t} = \frac{\partial}{\partial \tau}
+\end{equation}
+$\Rightarrow$
+\begin{equation}
+    \frac{\partial \Ew}{\partial z} = i\beta(\omega)\Ew + i \frac{\omega^2}{2c\epsilon_0 \beta(\omega)}\PNL(z, \omega)
+\end{equation}
+
+\end{document}
\ No newline at end of file
diff --git a/src/scgenerator/cli/cli.py b/src/scgenerator/cli/cli.py
index fa7c374..75e06de 100644
--- a/src/scgenerator/cli/cli.py
+++ b/src/scgenerator/cli/cli.py
@@ -2,6 +2,7 @@ import argparse
 import os
 from collections import ChainMap
 from pathlib import Path
+import re
 
 import numpy as np
 
@@ -72,6 +73,7 @@ def create_parser():
         help="comma-separated list of left limit, right limit and unit. "
         "One plot is made for each limit set provided. Example : 600,1200,nm or -2,2,ps",
     )
+    plot_parser.add_argument("--options", "-o", default=None)
     plot_parser.set_defaults(func=plot_all)
 
     dispersion_parser = subparsers.add_parser(
@@ -180,8 +182,11 @@ def resume_sim(args):
 
 
 def plot_all(args):
+    opts = {}
+    if args.options is not None:
+        opts |= dict([o.split("=")[:2] for o in re.split("[, ]", args.options)])
     root = Path(args.sim_dir).resolve()
-    scripts.plot_all(root, args.spectrum_limits)
+    scripts.plot_all(root, args.spectrum_limits, **opts)
 
 
 def plot_init_field_spec(args):
diff --git a/src/scgenerator/physics/materials.py b/src/scgenerator/physics/materials.py
index da091bc..7114626 100644
--- a/src/scgenerator/physics/materials.py
+++ b/src/scgenerator/physics/materials.py
@@ -196,11 +196,11 @@ def ionization_rate_ADK(
     nstar = Z * np.sqrt(2.1787e-18 / ionization_energy)
     omega_t = lambda field: e * np.abs(field) / np.sqrt(2 * me * ionization_energy)
     Cnstar = 2 ** (2 * nstar) / (scipy.special.gamma(nstar + 1) ** 2)
-    omega_C = omega_p / 4 * math.abs2(Cnstar)
+    omega_pC = omega_p * Cnstar
 
     def rate(field: np.ndarray) -> np.ndarray:
-        opt4 = 4 * omega_t(field) / om
-        return omega_C * (4 * omega_p / ot) ** (2 * nstar - 1) * np.exp(-4 * omega_p / (3 * ot))
+        opt4 = 4 * omega_p / omega_t(field)
+        return omega_pC * opt4 ** (2 * nstar - 1) * np.exp(-opt4 / 3)
 
     return rate
 
@@ -229,7 +229,8 @@ class Plasma:
         """
         Ne = free_electron_density(self.t, field, N0, self.rate)
         return cumulative_trapezoid(
-            np.gradient(Ne, self.t) * self.Ip / field, self.t, initial=0
-        ) + e ** 2 / me * cumulative_trapezoid(
-            cumulative_trapezoid(Ne * field, self.t, initial=0), self.t, initial=0
+            np.gradient(Ne, self.t) * self.Ip / field
+            + e ** 2 / me * cumulative_trapezoid(Ne * field, self.t, initial=0),
+            self.t,
+            initial=0,
         )
diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py
index aa37c29..2d6195a 100644
--- a/src/scgenerator/physics/pulse.py
+++ b/src/scgenerator/physics/pulse.py
@@ -125,6 +125,26 @@ def modify_field_ratio(
     return ratio
 
 
+def convert_field_units(envelope: np.ndarray, n: np.ndarray, A_eff: float) -> np.ndarray:
+    """[summary]
+
+    Parameters
+    ----------
+    envelope : np.ndarray, shape (n,)
+        complex envelope in units such that |envelope|^2 is in W
+    n : np.ndarray, shape (n,)
+        refractive index
+    A_eff : float
+        effective mode field area in m^2
+
+    Returns
+    -------
+    np.ndarray, shape (n,)
+        real field in V/m
+    """
+    return 2 * envelope.real / np.sqrt(2 * units.epsilon0 * units.c * n * A_eff)
+
+
 def conform_pulse_params(
     shape: Literal["gaussian", "sech"],
     width: float = None,
diff --git a/src/scgenerator/plotting.py b/src/scgenerator/plotting.py
index 103ac96..f3e3dec 100644
--- a/src/scgenerator/plotting.py
+++ b/src/scgenerator/plotting.py
@@ -8,6 +8,8 @@ import numpy as np
 from matplotlib.colors import ListedColormap
 from scipy.interpolate import UnivariateSpline
 
+from .logger import get_logger
+
 from . import io, math
 from .defaults import default_plotting as defaults
 from .math import abs2, make_uniform_1D, span
@@ -250,7 +252,7 @@ def _finish_plot_2D(
     file_name,
     file_type,
 ):
-
+    logger = get_logger(__name__)
     # apply log transform if required
     if log is not False:
         vmax = defaults["vmax"] if vmax is None else vmax
@@ -272,7 +274,7 @@ def _finish_plot_2D(
         elif log == "unique 1D":
             try:
                 ref = _finish_plot_2D.ref
-                print(f"recovered reference value {ref} for log plot")
+                logger.info(f"recovered reference value {ref} for log plot")
             except AttributeError:
                 ref = np.max(values, axis=1)
                 ind = np.argmax((ref[:-1] - ref[1:]) < 0)
@@ -335,7 +337,7 @@ def _finish_plot_2D(
 
     if is_new_plot:
         fig.savefig(out_path, bbox_inches="tight", dpi=200)
-        print(f"plot saved in {out_path}")
+        logger.info(f"plot saved in {out_path}")
     if cbar_label is not None:
         return fig, (ax, cbar.ax)
     else:
diff --git a/src/scgenerator/scripts/__init__.py b/src/scgenerator/scripts/__init__.py
index 41e8456..047d0ab 100644
--- a/src/scgenerator/scripts/__init__.py
+++ b/src/scgenerator/scripts/__init__.py
@@ -6,6 +6,7 @@ from cycler import cycler
 
 import matplotlib.pyplot as plt
 import numpy as np
+from tqdm import tqdm
 
 from ..utils.parameter import BareParams
 
@@ -22,24 +23,26 @@ def fingerprint(params: BareParams):
     return h1, h2
 
 
-def plot_all(sim_dir: Path, limits: list[str]):
-    for p in sim_dir.glob("*"):
-        if not p.is_dir():
-            continue
-
-        pulse = Pulse(p)
-        for lim in limits:
-            left, right, unit = lim.split(",")
-            left = float(left)
-            right = float(right)
-            pulse.plot_2D(
-                left,
-                right,
-                unit,
-                file_name=p.parent
-                / f"{pretty_format_from_file_name(p.name)} {left} {right} {unit}",
-            )
-        plt.close("all")
+def plot_all(sim_dir: Path, limits: list[str], **opts):
+    dir_list = list(p for p in sim_dir.glob("*") if p.is_dir())
+    limits = [
+        tuple(func(el) for func, el in zip([float, float, str], lim.split(","))) for lim in limits
+    ]
+    print(limits)
+    with tqdm(total=len(dir_list) * len(limits)) as bar:
+        for p in dir_list:
+            pulse = Pulse(p)
+            for left, right, unit in limits:
+                pulse.plot_2D(
+                    left,
+                    right,
+                    unit,
+                    file_name=p.parent
+                    / f"{pretty_format_from_file_name(p.name)} {left} {right} {unit}",
+                    **opts,
+                )
+                bar.update()
+            plt.close("all")
 
 
 def plot_init_field_spec(