diff --git a/Makefile b/Makefile
index 339a773..e3e90a0 100644
--- a/Makefile
+++ b/Makefile
@@ -88,7 +88,7 @@ endif
 
 .PHONY: all clean install docs
 
-all: $(BINARIES) $(LIBRARIES) docs
+all: $(BINARIES) $(LIBRARIES)
 
 # Remove all generated files
 clean:
@@ -102,6 +102,9 @@ install: $(BINARIES) $(LIBRARIES) $(SHAREDIR)/doc $(SHAREDIR)/gray.1
 	install -Dm644 -t $(PREFIX)/share/doc/res $(SHAREDIR)/doc/res/*
 	install -Dm644 -t $(PREFIX)/share/man/man1 $(SHAREDIR)/gray.1
 
+# dependencies
+$(OBJDIR)/%.o: $(OBJDIR)/%.d
+
 # GRAY binary
 $(GRAY): $(shell ./depend $(OBJDIR)/main.o) | $(BINDIR)
 	$(LD) $(LDFLAGS) -o '$@' $^
diff --git a/doc/1.intro.md b/doc/1.intro.md
index 6ec45c0..053c852 100644
--- a/doc/1.intro.md
+++ b/doc/1.intro.md
@@ -106,9 +106,9 @@ references:
     issued: 2003-3-5, 2003-4-5
 
 # Font
-mainfont: Libertinus Serif
-mathfont: Libertinus Math
-monofont: Fira Mono
+# mainfont: Libertinus Serif
+# mathfont: Libertinus Math
+# monofont: Fira Mono
 
 # PDF output options
 classoptions:
diff --git a/doc/2.physics.md b/doc/2.physics.md
index 121f5d0..b5fa18c 100644
--- a/doc/2.physics.md
+++ b/doc/2.physics.md
@@ -20,11 +20,13 @@ $z=Z$.
 
 In the complex eikonal framework, the a solution of the wave equation for the
 electric field is looked for in the form
+
 $$
   {\bf E}({\bf x},t) =
     {\bf e}({\bf x}) E_0({\bf x})
     e^{-i k_0 S({\bf x}) + iωt}
 $$ {#eq:eikonal-ansatz}
+
 such that it allows for Gaussian beam descriptions.
 In [@eq:eikonal-ansatz], $ω$ is the real frequency, $k_0 = ω/c$ the
 wavevector amplitude in vacuum, ${\bf e}({\bf x})$ the polarisation versor and
@@ -35,6 +37,7 @@ The function $S({\bf x})$ is the complex eikonal, $S = S_R({\bf x}) + i S_I
 propagation as in the geometric optics (GO), and the imaginary part $S_I({\bf
 x}) (<0)$ to the beam intensity profile shape, as it is apparent writing
 [@eq:eikonal-ansatz] as
+
 $$
   {\bf E}({\bf x},t) =
     {\bf e}({\bf x}) E_0({\bf x})
@@ -50,6 +53,7 @@ the $\bar x$ axis lies in the horizontal plane (i.e., $z=\text{const}$), and
 two additional coordinate systems, $(\xi_w,\eta_w)$ and $(\xi_R,\eta_R)$ in the
 $(\bar x, \bar y)$ plane, rotated by the angles $φ_w$ and $φ_R$,
 respectively,
+
 $$
 \begin{aligned}
   \bar x &= \xi_w \cos φ_w - \eta_w \sin φ_w
@@ -58,10 +62,12 @@ $$
           = \xi_R \sin φ_R + \eta_R \cos φ_R
 \end{aligned}
 $$ {#eq:phiwr}
+
 In the $(\xi_w,\eta_w)$ and $(\xi_R,\eta_R)$ systems, the axes are aligned
 with the major and minor axes of the intensity and phase ellipses respectively,
 and the general astigmatic Gaussian beam in vacuum takes the simple form
 [@gaussian-beam]
+
 $$
   E ({\bf x}) \propto
     \exp{\left[- \left(\frac{{\xi}_w^2}{w_\xi^2}
@@ -82,6 +88,7 @@ $w_{\xi,\eta}$, $R_{c\xi,\eta}$, $φ$ or alternatively by the beam waists
 $w_{0\xi,\eta}$, the waists $\bar z$ coordinates $d_{0\xi,\eta}$, and $φ$,
 where $R_{c\xi,\eta}$, $w_{\xi,\eta}$ are related to $d_{0\xi,\eta}$,
 $w_{0\xi,\eta}$ by the following equations:
+
 $$
   \begin{aligned}
     R_{cj} &= [({\bar z}- d_{0j})^2+z_{Rj}^2]/({\bar z}- d_{0j}) \\
@@ -102,6 +109,7 @@ in detail in [@gray].
 The "extended" rays obey to the following quasi-optical ray-tracing equations
 that are coupled together through an additional constraint in the form of a
 partial differential equation:
+
 $$
 \begin{aligned}
   \frac{d {\bf x}}{dσ} &=
@@ -111,24 +119,27 @@ $$
   \frac{∂ Λ}{∂ {\bf N}} &\cdot ∇ S_I = 0
 \end{aligned}
 $$
+
 where the function $Λ ({\bf x},{\bf k},ω)$ is the QO dispersion
 relation, which reads
+
 $$
   Λ  = N² -  N_c²({\bf x}, N_\parallel, ω)
           - |∇ S_I|² + \frac{1}{2}(\mathbf{b} ⋅ ∇ S_I )^2
                        \frac{∂² N_s²}{∂{N_\parallel}²} = 0
 $$ {#eq:eqlam}
+
 being $\mathbf{b}=\mathbf{B}/B$, $N_\parallel = {\mathbf N} \cdot \mathbf{b}$,
 and $N_c({\bf x}, N_\parallel, ω)$ the solution of the cold dispersion relation
 for the considered mode.
 
 In GRAY three choices for the integration variable $σ$ are available, i.e.:
-
 1. the arclength along the trajectory $s$,
 2. the time $τ=ct$, and
 3. the real part of the eikonal function $S_R$.
 
 The default option is the variable $s$ and the QO ray equations become:
+
 $$
 \begin{aligned}
   \frac{d{\bf x}}{ds} &=
@@ -141,6 +152,7 @@ $$
 \end{aligned}
 $$ {#eq:qort}
 
+
 ## Ray initial conditions
 
 The QO ray equations [@eq:qort] are solved for $N_T= N_r \times N_\vartheta
@@ -170,6 +182,7 @@ z_{Rj}$, the ray distribution used for the QO ray-tracing.
 The launching coordinates of the central ray of the EC beam will be denoted
 either as $(x_0, y_0, z_0)$, or $(R_0, φ_0, Z_0)$, depending on the
 coordinate system used (cartesian or cylindrical)
+
 $$
   \begin{aligned}
     x_0 &= R_0\cosφ_0 \\
@@ -177,8 +190,10 @@ $$
     z_0 &= Z_0.
   \end{aligned}
 $$
+
 and the launched wavevector $\bf N$ will have components $(N_{x0}, N_{y0},
 N_{z0})$, and $(N_{R0}, N_{φ 0}, N_{Z0})$, related by
+
 $$
   \begin{aligned}
     N_{x0} &= N_{R0} \cosφ_0 - N_{φ 0} \sinφ_0, \\
@@ -192,6 +207,7 @@ $$
 
 The poloidal and toroidal angles $α, β$ are defined in terms of the
 cylindrical components of the wavevector
+
 $$
   \begin{aligned}
     N_{R0}  &= -\cosβ \cosα, \\
@@ -199,7 +215,9 @@ $$
     N_{Z0}  &= -\cosβ \sinα
   \end{aligned}
 $$ {#eq:ncyl}
+
 with $-180° ≤ α ≤ 180°$, and $-90° ≤ β ≤ 90°$, so that
+
 $$
   \begin{aligned}
     \tanα &= N_{Z0}/N_{R0}, \\
@@ -221,10 +239,13 @@ This convention is the same used for the EC injection angles in ITER
 
 The EC power $P$ is assumed to evolve along the ray trajectory obeying to the
 following equation
+
 $$
   \frac{dP}{ds} = -α P,
 $$ {#eq:pincta}
+
 where here $α$ is the absorption coefficient
+
 $$
   α = 2 \frac{ω}{c} \frac {{\text{Im}}(Λ_w)}
         {|∂ Λ /∂ {\bf{N}}|} \biggr|_{Λ=0}
@@ -232,8 +253,10 @@ $$
         \frac {N_{\perp}} {|{∂ Λ}/{∂ {\bf N}|}}  \biggr|_{Λ=0}
     = 2{{\text{Im}}(k_{\perp w})} \frac{v_{g\perp}} v_{g}.
 $$ {#eq:alpha}
+
 being $N_{\perp w}$ (and $k_{\perp w}$) the perpendicular refractive index (and
 wave vector) solution of the relativistic dispersion relation for EC waves
+
 $$
   Λ_w = N^2-N_{\parallel}^2-N_{\perp w}^2=0
 $$
@@ -244,11 +267,13 @@ described in [@dispersion].
 
 Integration of [@eq:pincta] yields the local transmitted and deposited
 power in terms of the optical depth $τ= \int_0^{s}{α(s') d s'}$ as
+
 $$
   P(s)=P_0 e^{-τ(s)},
   \quad \mathrm{and} \quad
   P_{abs} (s)=P_0 [1-e^{-τ}] ,
 $$
+
 respectively, being $P_0$ the injected power.
 
 The flux surface averaged absorbed power density $p(ρ)=dP_{abs}/dV$ is
@@ -256,9 +281,11 @@ computed as the the ratio between the power deposited within the volume $dV$
 between two adjacent flux surfaces and the volume itself. At each position
 along the ray trajectory (parametrized by $s$), the absorbed power density can
 be written in terms of the absorption coefficient as
+
 $$
   p = P₀ α(s) e^{-τ(s)} \frac{δs}{δV}
 $$ {#eq:pav}
+
 $δs$ being the ray length between two adjacent magnetic surfaces, and $δV$ the
 associated volume.
 
@@ -267,17 +294,21 @@ associated volume.
 
 Within the framework of the linear adjoint formulation, the flux surface
 averaged EC driven current density is given by
+
 $$
   \langle J_{\parallel}\rangle = {\mathcal R}^* \, p
 $$ {#eq:jav}
+
 where
 ${\mathcal R}^*$ is a current drive efficiency, which can be expressed as a ratio
 between two integrals in momentum space
+
 $$
   {\mathcal R}^*= \frac{e}{m c \nu_c} \frac{\langle B \rangle}{B_m}
   \frac{\int{d{\bf u} {\mathcal P}({\bf u}) \,
   \eta_{\bf u}({\bf u})}}{\int{d{\bf u} {\mathcal P}({\bf u}) }}
 $$ {#eq:effr}
+
 where $\nu_c=4 \pi n e^4 Λ_c/(m^2 c^3)$ is the collision frequency, with
 $Λ_c$ the Coulomb logarithm, and $B_m$, $\langle B \rangle$ are the
 minimum value and the flux surface averaged value of the magnetic field on the
@@ -292,22 +323,28 @@ expression is related to the chosen ECCD model.
 
 The flux surface average driven current density [@eq:jav] can be written as
 [@gray]
+
 $$
   \langle J_{\parallel}\rangle =
     P_0 α(s) e^{-τ(s)} {\mathcal R}^*(s) \frac{δs}{δV}
 $$ {#eq:jrtav}
+
 and the equation for the current evolution $I_{cd}$ along the ray trajectory as
+
 $$
   \frac{dI_{cd}}{ds} =
     -{\mathcal R}^*(s)\frac{1}{2 \pi R_J } \frac{dP}{ds},
 $$
+
 where $R_J(\psi)$ is an effective radius for the computation of the driven
 current
+
 $$
   \frac{1}{R_J}
     = \langle \frac{1}{R^2} \rangle \frac{f(\psi)}{ \langle B\rangle}
     = \frac{ \langle {B_φ}/{R} \rangle}{ \langle B\rangle}
 $$
+
 being $f(\psi) =B_φ R$ the poloidal flux function.
 
 
@@ -327,31 +364,37 @@ trapping is based on a local development.
 In GRAY, three outputs for the EC driven current density are given.
 The EC flux surface averaged driven *parallel* current density $\langle
 J_{\parallel}\rangle$, that is the output of the ECCD theory, defined as
+
 $$
   \langle J_{\parallel}\rangle
     = \left \langle\frac{{\bf J}_{cd} \cdot {\bf B}}{B} \right \rangle
     = \frac{\langle {{\bf J}_{cd} \cdot {\bf B}}\rangle}
            {{\langle B^2 \rangle/}{\langle B \rangle}}.
 $$
+
 a *toroidal* driven current density $J_φ$ defined as
-\begin{equation}
+
+$$
   J_φ =\frac{δ I_{cd}} {δ A}
-  \label{eq:jphia}
-\end{equation}
+$$ {#eq:jphia}
+
 being $δ I_{cd}$ the current driven within the volume $δ V$ between
 two adjacent flux surfaces, and $δ A$ the poloidal area between the two
 adjacent flux surfaces, such that the total driven current is computed as
 $I_{cd}= \int J_φ dA$.
 Finally, an EC flux surface averaged driven current density $J_{cd}$ to be
 compared with transport code outputs
+
 $$
   J_{cd} = \frac{\langle {\bf J} \cdot {\bf B} \rangle} {B_{ref}}
 $$ {#eq:jcd}
+
 with the $B_{ref}$ value dependent on the transport code, i.e, $B_{ref}=B_0$
 for ASTRA and CRONOS, and $B_{ref}={\langle B \rangle}$ for JINTRAC.
 
 The above definitions are related to each other in terms of flux surface
 averaged quantities, dependent on the equilibrium, i.e.,
+
 $$
   \begin{aligned}
     J_φ &= \frac{f(\psi)}{\langle B \rangle}
@@ -383,22 +426,27 @@ via [@eq:ratj] for the two $B_{ref}$ choices.
 The second approach applies also to non monotonic profiles. Two average
 quantities are computed for both power and current density profiles, namely,
 the average radius $\langle ρ \rangle_a$ $(a=p,j)$
+
 $$
   \langle ρ \rangle_p = \frac{\int dV ρ p(ρ)}{\int dV  p(ρ)} , \qquad
   \langle ρ \rangle_j = \frac{\int dA ρ | J_{φ}(ρ)|}  {\int dA |J_{φ}(ρ)|}
 $$ {#eq:rav}
+
 and average profile width ${δρ}_a$ defined in terms of the variance as
+
 $$
   δ ρ_a = 2 \sqrt{2} \langle δ ρ \rangle_a
   \qquad  \mathrm {with } \qquad
   \langle δ ρ \rangle_a^2 = \langle ρ^2 \rangle_a-(\langle ρ \rangle_a)^2
 $$ {#eq:drav}
+
 Factor $\sqrt{8}$ is introduced to match with the definition of the full
 profile width in case of Gaussian profiles.
 Consistently with the above average definitions, we introduce suitable peak
 values $p_{0}$ and $J_{φ 0}$, corresponding to those of a Gaussian profile
 characterized by [@eq:rav;@eq:drav] and same total absorbed power $P_{abs}$ and
 driven current $I_{cd}$
+
 $$
   p_0 = \frac{2}{\sqrt{\pi}} \frac{P_{abs}}{{ δ ρ}_p
         \left ({dV}/{d ρ}\right)_{\langle ρ \rangle_p}},
@@ -417,20 +465,23 @@ conductor is assumed for the reflecting surface, so that the full power of the
 incident beam is transferred to the reflected one. The vector refractive index
 ${\bf N}_{\rm{refl}}$ and the unit electric field $\hat {\bf e}_{\rm{refl}}$ of
 the reflected wave are
-\begin{equation}
+
+$$
   {\bf N}_{\rm{refl}} =
     {\bf N}_{\rm{in}} - 2 ({\bf N}_{\rm{in}}
     \cdot \hat {\bf n}) \hat {\bf n}, \qquad
   \hat {\bf e}_{\rm{refl}} =
     -\hat {\bf e}_{\rm{in}}
     + 2 (\hat {\bf e}_{\rm{in}} \cdot \hat {\bf n}) \hat {\bf n},
-\end{equation}
+$$
+
 being ${\bf N}_{\rm{in}}$ and $\hat {\bf e}_{\rm{in}}$ the vector refractive
 index and the unit electric field of the incoming wave, and $\hat {\bf n}$ the
 normal unit vector to the wall at the beam incidence point.
 
 The Stokes parameter for the unit electric vector $\hat {\bf e}$ in vacuum are
 defined in the beam reference system $({\bar x},{\bar y},{\bar z})$ as
+
 $$
   \begin{aligned}
     I &= \vert \hat e_{\bar x} \vert^2 + \vert \hat e_{\bar y} \vert^2 = 1 \\
@@ -439,7 +490,9 @@ $$
     V &= 2 \cdot {\rm Im} (\hat e_{\bar x} \hat e_{\bar y}^*).
   \end{aligned}
 $$ {#eq:stokes}
+
 Alternatively, the two angles $\psi_p$ and $\chi_p$ can be used:
+
 $$
   \begin{aligned}
     Q &= \cos {2 \psi_p} \cos {2 \chi_p} \\
@@ -447,12 +500,14 @@ $$
     V &= \sin {2 \chi_p}
   \end{aligned}
 $$
+
 which define respectively the major axis orientation and the ellipticity of the
 polarisation ellipse. The polarisation parameters of the reflected wave are
 used to compute the coupling with the Ordinary (OM) and Extraordinary (XM)
 modes at the vacuum-plasma interface before the calculation of the second pass
 in the plasma. At the second pass both modes are traced, taking into account
 that the power fraction coupled to each mode is
+
 $$
   P_{\rm O,X} =
     \frac{P_{\rm in}}{2}
diff --git a/doc/3.io-files.md b/doc/3.io-files.md
index c74fc8d..006c57d 100644
--- a/doc/3.io-files.md
+++ b/doc/3.io-files.md
@@ -197,6 +197,7 @@ Table: **Plasma boundary** {#tbl:eqdisk-bound}
 ### Toroidal Current Density
 
 The toroidal current $J_T$ (A/m²) is related to $P'(ψ)$ and $FF'(ψ)$ through
+
 $$
   J_T = R P'(ψ) + FF'(ψ) / R
 $$
diff --git a/doc/4.implementation.md b/doc/4.implementation.md
index 0a52b83..6c2ea8d 100644
--- a/doc/4.implementation.md
+++ b/doc/4.implementation.md
@@ -16,6 +16,7 @@ modes and the `index_rt` is updated as:
 It follows that ordinary(extraordinary) modes respectively have odd(even)
 indices and the number of passes is given by $\lfloor \log₂(1 + \tt index\_rt)
 \rfloor$. For example, an `index_rt`=19 denotes the following chain:
+
 $$
   \begin{aligned}
     \text{mode:}           && O &→ X → O → O  \\