use MathML for the HTML manual

All the major browser now supports it pretty well.

README.md

@@ -76,18 +76,16 @@ To build the documentation (GRAY user manual and man pages) you need:
 - [pandoc] (≥2.19 tested)
 - [librsvg]
 - a basic [Tex Live] installation (≥2022 tested)
-- a [KaTeX] release archive
 - [help2man]
 
 [pandoc]: https://pandoc.org/installing.html
 [TeX Live]: https://tug.org/texlive/quickinstall.html
 [librsvg]: https://wiki.gnome.org/Projects/LibRsvg
-[KaTeX]: https://github.com/KaTeX/KaTeX/releases
 [help2man]: https://www.gnu.org/software/help2man/
 
 Then run:
 ```
-./configure --with-katex=/path/to/katex/archive/unpacked
+./configure
 make docs
 ```
 The results will be placed in `build/share`, to install, run:
@@ -95,6 +93,9 @@ The results will be placed in `build/share`, to install, run:
 make install-doc
 ```
 
+Note: you may have to specify a math font to be bundled with the manual using
+`./configure --with-math-font=/path/to/font.otf`.
+
 ## Running the tests
 
 GRAY includes a number of integration tests that are useful to check whether

configure

@@ -17,7 +17,7 @@ Options:
   --enable-static         statically link programs and libraries [no]
   --disable-static        dynamically link programs and libraries [yes]
   --enable-deterministic  try to make a bit-for-bit deterministic build [no]
-  --with-katex=PATH       specify path to the KaTeX library [automatically downloaded]
+  --with-math-font=PATH   specify path to OpenType math font for docs [Libertinus Math]
 EOF
 }
 
@@ -44,7 +44,7 @@ for arg in "$@"; do
     --enable-static) printf 'STATIC=1\n' >> configure.mk ;;
     --disable-static) ;;
     --enable-deterministic) printf 'DETERMINISTIC=1\n' >> configure.mk ;;
-    --with-katex=*) printf 'KATEX_URL==file://%s/\n' "${arg#*=}" >> configure.mk ;;
+    --with-math-font=*) printf 'MATHFONT=%s\n' "${arg#*=}" >> configure.mk ;;
     *=*)
       printf '%s\n' "${arg?}" >> configure.mk
       printf 'set %s\n' "$arg"

default.nix

@@ -11,12 +11,6 @@ let
       sha256 = "166yxg4ff2jxvl9mbngd90dr1k3rdj6xmiccga41xynhh2zr0vmb";
     };
 
-  # Needed for HTML manual
-  katex = builtins.fetchTarball
-    { url = "https://github.com/KaTeX/KaTeX/releases/download/v0.15.1/katex.tar.gz";
-      sha256 = "007nv11r0z9fz593iwzn55nc0p0wj5lpgf0k2brhs1ynmikq9gjr";
-    };
-
   # Exclude this file, .git and more from source
   sieve = path: type:
     !builtins.elem (baseNameOf path) (
@@ -66,7 +60,7 @@ in
       makefile2graph graphviz
 
       # documentation
-      help2man pandoc librsvg
+      help2man pandoc librsvg fontconfig
       haskellPackages.pandoc-crossref
       (texlive.combine {
         inherit (texlive) scheme-small xetex fontspec;
@@ -95,7 +89,6 @@ in
     hardeningDisable = [ "format" ];
     configureFlags = [
       (lib.enableFeature static "static")
-      "--with-katex=${katex}"
       "--enable-deterministic"
     ];
 

doc/2.physics.md

@@ -12,7 +12,7 @@ viewed from above.
 
 In addition to the right handed Cartesian orthogonal system specified above, we
 introduce also a right-handed cylindrical system $(R,φ,Z)$ with transformation
-from the cylindrical to the Cartesian system given by $x= R\cosφ$, $y=R\sinφ$,
+from the cylindrical to the Cartesian system given by $x= R\cos φ$, $y=R\sin φ$,
 $z=Z$.
 
 
@@ -20,28 +20,25 @@ $z=Z$.
 
 In the complex eikonal framework, a solution of the wave equation for the
 electric field is looked for in the form
-
 $$
-  {\bf E}({\bf x},t) =
+  {\mathbf E}({\mathbf x},t) =
-    {\bf e}({\bf x}) E_0({\bf x})
+    {\mathbf e}({\mathbf x}) E_0({\mathbf x})
-    e^{-i k_0 S({\bf x}) + iωt}
+    e^{-i k_0 S({\mathbf 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 normalised polarisation (Jones)
+amplitude in vacuum, ${\mathbf e}({\mathbf x})$ the normalised polarisation (Jones)
-vector and $E_0({\bf x})$ the slowly varying wave amplitude.
+vector and $E_0({\mathbf x})$ the slowly varying wave amplitude.
 
-The function $S({\bf x})$ is the complex eikonal, $S = S_R({\bf x}) + i S_I
+The function $S({\mathbf x})$ is the complex eikonal, $S = S_R({\mathbf x}) + i S_I
-({\bf x})$, in which the real part $S_R({\bf x})$ is related to the beam
+({\mathbf x})$, in which the real part $S_R({\mathbf x})$ is related to the beam
-propagation as in the geometric optics (GO), and the imaginary part $S_I({\bf
+propagation as in the geometric optics (GO), and the imaginary part $S_I({\mathbf
 x}) (<0)$ to the beam intensity profile shape, as it is apparent writing
 [@eq:eikonal-ansatz] as
-
 $$
-  {\bf E}({\bf x},t) =
+  {\mathbf E}({\mathbf x},t) =
-    {\bf e}({\bf x}) E_0({\bf x})
+    {\mathbf e}({\mathbf x}) E_0({\mathbf x})
-    e^{k_0 S_I({\bf x})} e^{-i k_0 S_R({\bf x})+i ωt}
+    e^{k_0 S_I({\mathbf x})} e^{-i k_0 S_R({\mathbf x})+i ωt}
 $$ {#eq:efri}
 
 
@@ -53,7 +50,6 @@ 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
@@ -67,9 +63,8 @@ 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
+  E ({\mathbf x}) \propto
     \exp{\left[- \left(\frac{{\xi}_w^2}{w_\xi^2}
           +\frac{{\eta}_w^2}{w_\eta^2}\right)
          \right]}
@@ -83,19 +78,17 @@ parameters: the beam widths $w_{\xi,\eta}$, the phase front curvature radii
 $R_{c\xi,\eta}$ and the intensity and phase ellipses rotation angles $φ_{w,R}$.
 
 Simple astigmatic beams can be described in terms of 5 parameters only, because
-the phase and intensity ellipses are aligned, i.e., $φ_w=φ_R\equivφ$:
+the phase and intensity ellipses are aligned, i.e., $φ_w=φ_R\equiv φ$:
 $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}) \\
     w_j  &= w_{0j} \sqrt{1+({\bar z}- d_{0j})^2/z_{Rj}^2},
   \end{aligned}
 $$ {#eq:rciw}
-
 and $z_{Rj}= k_0 w_{0j}^2/2$ is the Raylegh length. According to
 [@eq:rciw], a convergent beam (${\bar z} < d_{0j}$) has $R_{cj}<0$, while a
 divergent beam has $R_{cj}>0$.
@@ -109,28 +102,23 @@ 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σ} &=
+  \frac{d {\mathbf x}}{dσ} &= +\frac{∂Λ}{∂\mathbf N} \biggr |_{Λ=0} \\
-    +{∂ Λ \over ∂ {\bf N}} \biggr |_{Λ=0} \\
+  \frac{d {\mathbf N}}{dσ} &= -\frac{∂Λ}{∂\mathbf x} \biggr |_{Λ=0} \\
-  \frac{d {\bf N}}{dσ} &=
+  \frac{∂ Λ}{∂ {\mathbf N}} &\cdot ∇ S_I = 0
-    -{∂ Λ \over ∂ {\bf x}} \biggr |_{Λ=0} \\
-  \frac{∂ Λ}{∂ {\bf N}} &\cdot ∇ S_I = 0
 \end{aligned}
 $$
-
+where the function $Λ ({\mathbf x},{\mathbf k},ω)$ is the QO dispersion
-where the function $Λ ({\bf x},{\bf k},ω)$ is the QO dispersion
 relation, which reads
-
 $$
-  Λ  = N² -  N_c²({\bf x}, N_\parallel, ω)
+  Λ  = N² -  N_c²({\mathbf 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
+and $N_c({\mathbf 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.:
@@ -139,16 +127,15 @@ In GRAY three choices for the integration variable $σ$ are available, i.e.:
 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} &=
+  \frac{d{\mathbf x}}{ds} &=
-    +\frac{∂ Λ /∂ {\bf N}}
+    +\frac{∂ Λ /∂ {\mathbf N}}
-    {|∂ Λ /∂ {\bf N}|} \biggr |_{Λ=0}  \\
+    {|∂ Λ /∂ {\mathbf N}|} \biggr |_{Λ=0}  \\
-  \frac{d{\bf N}}{ds} &=
+  \frac{d{\mathbf N}}{ds} &=
-    -\frac{∂ Λ /∂ {\bf x}}
+    -\frac{∂ Λ /∂ {\mathbf x}}
-    {|∂ Λ /∂ {\bf N}|} \biggr |_{Λ=0} \\
+    {|∂ Λ /∂ {\mathbf N}|} \biggr |_{Λ=0} \\
-  \frac{∂ Λ}{∂ {\bf N}} &\cdot ∇ S_I = 0
+  \frac{∂ Λ}{∂ {\mathbf N}} &\cdot ∇ S_I = 0
 \end{aligned}
 $$ {#eq:qort}
 
@@ -182,46 +169,40 @@ 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 \\
+    x_0 &= R_0\cos φ_0 \\
-    y_0 &= R_0\sinφ_0 \\
+    y_0 &= R_0\sin φ_0 \\
     z_0 &= Z_0.
   \end{aligned}
 $$
-
+and the launched wavevector $\mathbf N$ will have components $(N_{x0}, N_{y0},
-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, \\
+    N_{x0} &= N_{R0} \cos φ_0 - N_{φ 0} \sin φ_0, \\
-    N_{y0} &= N_{R0} \sinφ_0 + N_{φ 0} \cosφ_0, \\
+    N_{y0} &= N_{R0} \sin φ_0 + N_{φ 0} \cos φ_0, \\
     N_{z0} &= N_{Z0}
   \end{aligned}
 $$
 
-
 ## EC Launching angles ($α,β$)
 
 The poloidal and toroidal angles $α, β$ are defined in terms of the
 cylindrical components of the wavevector
-
 $$
   \begin{aligned}
-    N_{R0}  &= -\cosβ \cosα, \\
+    N_{R0}  &= -\cos β \cos α, \\
-    N_{φ0}  &= +\sinβ,       \\
+    N_{φ0}  &= +\sin β,        \\
-    N_{Z0}  &= -\cosβ \sinα
+    N_{Z0}  &= -\cos β \sin α
   \end{aligned}
 $$ {#eq:ncyl}
 
 with $-180° ≤ α ≤ 180°$, and $-90° ≤ β ≤ 90°$, so that
-
 $$
   \begin{aligned}
-    \tanα &= N_{Z0}/N_{R0}, \\
+    \tan α &= N_{Z0}/N_{R0}, \\
-    \sinβ &= N_{φ 0}
+    \sin β &= N_{φ 0}
   \end{aligned}
 $$ {#eq:albt}
 
@@ -239,24 +220,19 @@ 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}
+        {|∂ Λ /∂ {\mathbf{N}}|} \biggr|_{Λ=0}
     ≈ 4 \frac{ω}{c} {{\text {Im}}(N_{\perp w})}
-        \frac {N_{\perp}} {|{∂ Λ}/{∂ {\bf N}|}}  \biggr|_{Λ=0}
+        \frac {N_{\perp}} {|{∂ Λ}/{∂ {\mathbf 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
 $$
@@ -267,13 +243,11 @@ 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
@@ -281,11 +255,9 @@ 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.
 
@@ -294,63 +266,54 @@ 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
-where
+a ratio between two integrals in momentum space
-${\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}) \,
+  \frac{\int{d{\mathbf u} {\mathcal P}({\mathbf u}) \,
-  \eta_{\bf u}({\bf u})}}{\int{d{\bf u} {\mathcal P}({\bf u}) }}
+  \eta_{\mathbf u}({\mathbf u})}}{\int{d{\mathbf u} {\mathcal P}({\mathbf 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
+$Λ_c$ the Coulomb logarithm, and $B_m$, $\langle B \rangle$ are the minimum
-minimum value and the flux surface averaged value of the magnetic field on the
+value and the flux surface averaged value of the magnetic field on the given
-given magnetic surface, respectively.
+magnetic surface, respectively.
-The functions ${\mathcal P}({\bf u})$ and $\eta_{\bf u}({\bf u})$ are the
+
-normalized absorbed power density and current drive efficiency per unit
+The functions ${\mathcal P}({\mathbf u})$ and $\eta_{\mathbf u}({\mathbf u})$
-momentum ${\bf u}={\bf p}/mc$ [@gray].
+are the normalized absorbed power density and current drive efficiency per unit
-Note that the warm wave polarisation is used to compute ${\mathcal P}({\bf u})$.
+momentum ${\mathbf u}={\mathbf p}/mc$ [@gray].
-In the adjoint formulation adopted here, the function $\eta_{\bf u}({\bf u})$
+Note that the warm wave polarisation is used to compute ${\mathcal P}({\mathbf
-is written in terms of the response function for the current, and its explicit
+u})$. In the adjoint formulation adopted here, the function $\eta_{\mathbf
-expression is related to the chosen ECCD model.
+u}({\mathbf u})$ is written in terms of the response function for the current,
+and its explicit 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}
+    = \left\langle\frac{1}{R^2}\right\rangle \frac{f(\psi)}{ \langle B\rangle}
     = \frac{ \langle {B_φ}/{R} \rangle}{ \langle B\rangle}
 $$
-
 being $f(\psi) =B_φ R$ the poloidal flux function.
 
 
 ## ECCD Models
 
-Two models for $\eta_{\bf u}({\bf u})$ efficiency in [@eq:effr] are implemented
+Two models for $\eta_{\mathbf u}({\mathbf u})$ efficiency in [@eq:effr] are implemented
 for ECCD calculations, a Cohen-like module in the high-velocity limit and the
 momentum conserving model developed by Marushenko.
 The used Cohen-like module, developed explicitly for GRAY, is described in
@@ -364,37 +327,31 @@ 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
+    = \left \langle\frac{{\mathbf J}_{cd} \cdot {\mathbf B}}{B} \right \rangle
-    = \frac{\langle {{\bf J}_{cd} \cdot {\bf B}}\rangle}
+    = \frac{\langle {{\mathbf J}_{cd} \cdot {\mathbf B}}\rangle}
            {{\langle B^2 \rangle/}{\langle B \rangle}}.
 $$
-
 a *toroidal* driven current density $J_φ$ defined as
-
 $$
   J_φ =\frac{δ I_{cd}} {δ A}
 $$ {#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}}
+  J_{cd} = \frac{\langle {\mathbf J} \cdot {\mathbf 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}
@@ -426,17 +383,14 @@ 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
+  \qquad  \text{with } \qquad
   \langle δ ρ \rangle_a^2 = \langle ρ^2 \rangle_a-(\langle ρ \rangle_a)^2
 $$ {#eq:drav}
 
@@ -446,7 +400,6 @@ 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}},
@@ -472,8 +425,8 @@ The beam vector at launch is computed from the polarisation ellipse parameters
 using the formula:
 $$
   \begin{aligned}
-    \hat{e}₁ &= \cosχ\cosψ + i\sinχ\sinψ \\
+    \hat{e}₁ &= \cos χ\cos ψ + i\sin χ\sin ψ \\
-    \hat{e}₂ &= \cosχ\sinψ - i\sinχ\cosψ
+    \hat{e}₂ &= \cos χ\sin ψ - i\sin χ\cos ψ
   \end{aligned}
 $$ {#eq:ellipse2field}
 
@@ -481,7 +434,7 @@ The following convention is assumed (illustrated in [@fig:ellipse]):
 
 - $ψ$ is the angle between the $x$ axis and the major axis.
 
-- $χ = \tan(b/a)$ where $a,b$ are ellipse major and minor semi-axes, respectively.
+- $χ{=}\tan(b/a)$ where $a,b$ are ellipse major and minor semi-axes, respectively.
 
 - A positive $ψ$ corresponds to an ellipse rotated counterclockwise in the
   $x∧y$ plane.
@@ -501,18 +454,18 @@ used to describe the beam reflection on the inner wall of the tokamak in the
 cases where only partial absorption occurs at the first pass in the plasma.
 An ideal 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 Jones vector $\hat {\bf
+refractive index ${\mathbf N}_\text{refl}$ and the Jones vector $\hat {\mathbf
-e}_{\rm{refl}}$ of the reflected beam are
+e}_\text{refl}$ of the reflected beam are
 $$
-  {\bf N}_{\rm{refl}} =
+  {\mathbf N}_\text{refl} =
-    {\bf N}_{\rm{in}} - 2 ({\bf N}_{\rm{in}}
+    {\mathbf N}_\text{in} - 2 ({\mathbf N}_\text{in}
-    \cdot \hat {\bf n}) \hat {\bf n}, \qquad
+    \cdot \hat {\mathbf n}) \hat {\mathbf n}, \qquad
-  \hat {\bf e}_{\rm{refl}} =
+  \hat {\mathbf e}_\text{refl} =
-    -\hat {\bf e}_{\rm{in}}
+    -\hat {\mathbf e}_\text{in}
-    + 2 (\hat {\bf e}_{\rm{in}} \cdot \hat {\bf n}) \hat {\bf n},
+    + 2 (\hat {\mathbf e}_\text{in} \cdot \hat {\mathbf n}) \hat {\mathbf n},
 $$
-being ${\bf N}_{\rm{in}}$ and $\hat {\bf e}_{\rm{in}}$ the vector refractive
+being ${\mathbf N}_\text{in}$ and $\hat {\mathbf e}_\text{in}$ the vector refractive
-index and the Jones vector of the incoming wave, and $\hat {\bf n}$ the
+index and the Jones vector of the incoming wave, and $\hat {\mathbf n}$ the
 normal unit vector to the wall at the beam incidence point.
 
 The reflected beam Jones vector is again used to compute the coupling
@@ -520,14 +473,14 @@ to the plasma modes at the second and successive pass, with potentially
 $2^{n-1}$ independent modes being traced after n reflections.
 
 Note that the Jones vectors of the ordinary and extraordinary modes are
-orthogonal w.r.t. the standard Hermitian product: $\hat{\mathbf e}_{\rm
+orthogonal w.r.t. the standard Hermitian product:
-O}⋅\hat{\mathbf e}_{\rm X}^* = 0$.
+$\hat{\mathbf{e}}_\mathrm{O}⋅\hat{\mathbf{e}}_\mathrm{X}^* = 0$.
 From [@eq:ellipse2field] it then follows that these relations hold:
 $$
   \begin{aligned}
-    ψ_{\rm O} &= ψ_{\rm X} + \frac{π}{2} \\
+    ψ_\textrm{O} &= ψ_\textrm{X} + \frac{π}{2} \\
-    χ_{\rm O} &= -χ_{\rm X} \\
+    χ_\textrm{O} &= -χ_\textrm{X} \\
-    1 &= c_{\rm O} + c_{\rm X}
+    1 &= c_\textrm{O} + c_\textrm{X}
   \end{aligned}
 $$
 with the latter meaning all the incoming power is coupled to the plasma.

doc/3.io-files.md

@@ -197,7 +197,6 @@ 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
 $$
@@ -252,12 +251,12 @@ Variable         Type          Units          Valid range       Definition
 `P0`             Real        MW         $x>0$                   EC injected power.
 
 `Nrayr`          Integer     1          $1 ≤ n ≤ 31$            Number of rays $N_r$ in radial direction + 1 in
-                                                                the center, ${\tt Nrayr} = N_r +1$.
+                                                                the center, $\texttt{Nrayr} = N_r +1$.
 
 `Nrayth`         Integer     1          $1 ≤ n ≤ 36$            Number of rays $N_{\theta}$ in angular direction.
 
 `rwmax`          Real                   $x > 0$                 "cut-off" size $\tilde ρ_{max}$ of the gaussian beam.
-                                                                (typical $1 ≤ {\tt rwmax} ≤ 1.5$), defined at page.
+                                                                (typical $1 ≤ \texttt{rwmax} ≤ 1.5$), defined at page.
 
 `x0`             Real        cm                                 $X$ coordinate of the launching point.
 
@@ -291,7 +290,7 @@ Variable         Type          Units          Valid range       Definition
                                                                   2=$ψ$.
 
 `filenmbm`       String                 len$(s)≤ 24$            Name of file (extension `.txt` appended) with beam.
-                                                                parameters, used if ${\tt ibeam} >0$.
+                                                                parameters, used if $\texttt{ibeam} >0$.
 
 `iox`            Integer                1, 2                    1=Ordinary mode (OM), 2=Extraordinary (XM) mode.
 
@@ -331,19 +330,19 @@ Variable         Type          Units          Valid range       Definition
                                                                 1=normalized (obsolete, for some old files
                                                                 psi$(R,z)$ in `filenmeqq`.
 
-`sspl`           Real        1              $x > 0$             Tension of spline fit for psi (${\tt sspl} \ll 1$ typical),
+`sspl`           Real        1              $x > 0$             Tension of spline fit for psi ($\texttt{sspl} \ll 1$ typical),
                                                                 0=interpolation.
 
 `factb`          Real        1              $x > 0$             Numerical factor to rescale the magnetic field
-                                                                $B \rightarrow B \cdot {\tt factb}$.
+                                                                $B \rightarrow B \cdot \texttt{factb}$.
 
 `factt`          Real        1              $x > 0$             Numerical factor to rescale the electron temperature
-                                                                $T_e \rightarrow T_e \cdot {\tt factt}
+                                                                $T_e \rightarrow T_e \cdot \texttt{factt}
-                                                                \cdot {\tt factb}^{\tt ab}$.
+                                                                \cdot \texttt{factb}^\texttt{ab}$.
 
 `factn`          Real        1              $x > 0$             Numerical factor to rescale the electron density
-                                                                $n_e \rightarrow n_e \cdot {\tt factn}
+                                                                $n_e \rightarrow n_e \cdot \texttt{factn}
-                                                                \cdot {\tt factb}^{\tt ab}$.
+                                                                \cdot \texttt{factb}^\text{ab}$.
 
 `iscal`          Int                        1, 2                Model for $n_e$, $T_e$ scaling with $B$:
                                                                   1=constant $n_{Greenwald}$ (`ab`=1),
@@ -352,7 +351,7 @@ Variable         Type          Units          Valid range       Definition
 `filenmprf`      String                     len$(s) ≤ 24$       Name of file for kinetic profiles
 
 `psdbnd`         Real        1              $x > 0$             Normalized psi value at the plasma boundary where
-                                                                $n_e$ is set to zero (typ. $1 ≤ {\tt psdbnd} ≤ 1.1$).
+                                                                $n_e$ is set to zero (typ. $1 ≤ \texttt{psdbnd} ≤ 1.1$).
 
 `sgnbphi`        Real        1              -1, +1, 0           Force signum of toroidal B, used if nonzero
 
@@ -374,7 +373,7 @@ Variable         Type          Units          Valid range       Definition
                                                                   3=fully relativistic (slow).
 
 `ilarm`          Integer     1          $n ≥ 1$                 Order of Larmor radius expansion for absorption
-                                                                computation ${\tt ilarm}>$ local EC harmonic number.
+                                                                computation $\texttt{ilarm}>$ local EC harmonic number.
 
 `ieccd`          Integer                0, 1, 11                Current drive model:
                                                                   0=none,
@@ -384,7 +383,7 @@ Variable         Type          Units          Valid range       Definition
 
 `igrad`          Integer                0, 1                    Ray-tracing model:
                                                                   0=optical,
-                                                                  1=quasi-optical (requires ${\tt nrayr} \ge 5$).
+                                                                  1=quasi-optical (requires $\texttt{nrayr} \ge 5$).
 
 `idst`           Integer                0, 1, 2                 Ray-tracing integration variable:
                                                                   0=$s$,
@@ -396,10 +395,10 @@ Variable         Type          Units          Valid range       Definition
 
 `nstep`          Integer     1          $0 ≤ n ≤ 8000$          Maximum number of integration steps.
 
-`istpr`          Integer     1          $1 ≤ n ≤ {\tt nstep}$   Subsampling factor for beam cross section
+`istpr`          Integer     1          $1 ≤ n ≤ \texttt{nstep}$   Subsampling factor for beam cross section
                                                                 data output (units 8, 12).
 
-`istpl`          Integer     1          $1 ≤ n ≤ {\tt nstep}$   Subsampling factor for outmost rays
+`istpl`          Integer     1          $1 ≤ n ≤ \texttt{nstep}$   Subsampling factor for outmost rays
                                                                 data output (unit 33).
 
 `ipec`           Integer                1, 2                    Grid spacing for ECRH&CD profiles:
@@ -408,7 +407,7 @@ Variable         Type          Units          Valid range       Definition
 
 `nnd`            Integer     1          $2 ≤ n ≤ 5001$          Number of points in the ECRH&CD profile grid.
 
-`ipass`          Integer                $-2,1, 2$               $\vert{\tt ipass}\vert$=number of passes into plasma:
+`ipass`          Integer                $-2,1, 2$               $\vert\texttt{ipass}\vert$=number of passes into plasma:
                                                                   -2=reflection at `rwallm`,
                                                                   +2=reflection at limiter.
                                                                 Surface given in EQDSK.
@@ -431,7 +430,7 @@ Variable         Type          Units          Valid range       Definition
 
 `ne`             real        10¹⁹m⁻³    $x ≥ 0$                 Electron density
 
-`Zeff`           real        1          $x ≥ 0$                 Effective charge, $Z_{\rm eff}$
+`Zeff`           real        1          $x ≥ 0$                 Effective charge, $Z_\text{eff}$
 ----------------------------------------------------------------------------------------------------------------------------
 
 Table: GRAY input data `filenmprf.prf` - plasma profiles {#tbl:profiles1}
@@ -525,11 +524,11 @@ Variable         Type          Units   Definition
                                        defined in [@eq:drav].
 
 `ratjbmx`        Real        1         Ratio $J_{cd}/J_{φ}$ at $ρ=$`rhotjava`,
-                                       with $J_{cd}=\langle \bold{J}\cdot \bold{B} \rangle
+                                       with $J_{cd}=\langle \mathbf{J}\cdot \mathbf{B} \rangle
                                                     /\langle \mathbf{B} \rangle$.
 
 `ratjamx`        Real        1         Ratio $J_{cd}/J_{φ}$ at $ρ=$`rhotjava`,
-                                       with $J_{cd}=\langle \bold{J}\cdot \bold{B} \rangle/B_0$
+                                       with $J_{cd}=\langle \mathbf{J}\cdot \mathbf{B} \rangle/B_0$
 
 `stmx`           Real        cm        Path from the launching point and peak $dP/dV$
                                        for the central ray.
@@ -555,7 +554,7 @@ Variable         Type          Units   Definition
                                        with $A$ the area of the poloidal section labelled by $ρ$.
 
 `Jcdb`           Real        MA⋅m⁻²    EC current density
-                                        $J_{cd}=\langle \bold{J}\cdot\bold{B} \rangle
+                                        $J_{cd}=\langle \mathbf{J}\cdot\mathbf{B} \rangle
                                                 /\langle \mathbf{B} \rangle$.
 
 `dPdV`           Real        MW⋅m⁻³    EC power density $p(ρ)=dP/dV$.

doc/4.implementation.md

@@ -14,12 +14,11 @@ 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)
+indices and the number of passes is given by $\lfloor \log₂(1 + \texttt{index\_rt})
 \rfloor$. For example, an `index_rt`=19 denotes the following chain:
-
 $$
   \begin{aligned}
     \text{mode:}         && O &→ X → O → O  \\
-    \text{\tt index\_rt:}  && 1 &→ 4 → 9 → 19
+    \texttt{index\_rt:}  && 1 &→ 4 → 9 → 19
   \end{aligned}
 $$

doc/Makefile

@@ -8,7 +8,8 @@
 # Pandoc flags
 FLAGS      = --filter pandoc-crossref --citeproc --toc --number-sections
 PDF_FLAGS  = --pdf-engine=xelatex
-HTML_FLAGS = --standalone --katex$(KATEX_URL) --lua-filter filters/add-anchors.lua
+HTML_FLAGS = --standalone --mathml --lua-filter filters/add-anchors.lua
+MATHFONT ?= $(shell fc-match --format=%{file} 'Libertinus Math')
 
 # Rebuild everything if the makefile changed
 .EXTRA_PREREQS += Makefile
@@ -34,6 +35,7 @@ $(DOCDIR)/manual.html: $(SECTIONS) res | $(DOCDIR)
 
 $(DOCDIR)/res: res | $(DOCDIR)
 	cp -Tfr $^ '$@'
+	cp $(MATHFONT) '$@/mathfont.otf'
 
 $(DOCDIR)/manual.pdf: $(SECTIONS) | $(DOCDIR)
 	pandoc $(FLAGS) $(PDF_FLAGS) $^ -o '$@'

doc/res/style.css

@@ -3,15 +3,29 @@
   box-sizing: border-box;
 }
 
+/* Set a compatible math font */
+@font-face {
+  font-family: "Libertinus Math";
+  src: local("Libertinus Math"),
+       url("res/mathfont.otf") format("opentype");
+}
+body { font-family: "Libertinus Math", sans-serif; }
+math { font-family: "Libertinus Math", math; }
+
 /* Justify all text */
 body {
   text-align: justify;
-  hypens: auto;
+  hyphens: auto;
+  word-wrap: break-word;
+  overflow-wrap: break-word;
 }
 
+/* Slighlty increase size */
+body { font-size: 1.1em }
 
 /* Make headings smaller */
-h1 { font-size: 1.8em }
+header > h1 { font-size: 2em }
+h1 { font-size: 1.7em }
 h2 { font-size: 1.4em }
 h3 { font-size: 1.1em }
 h1, h2, h3 { margin: 1em 0 }
@@ -37,16 +51,35 @@ h1:hover > .header-section-number,
 h2:hover > .header-section-number,
 h3:hover > .header-section-number { opacity: 0; }
 
+/* Fix equation alignment in chromium */
+mtd  {
+  margin-left: 0.1rem;
+  margin-right: 0.1rem;
+  padding-left: 0;
+  padding-right: 0;
+}
+mtr > mtd[columnalign="right"] { text-align: right; }
+mtr > mtd[columnalign="left"] { text-align: left; }
 
-/* Match math and normal mode fonts */
+/* Fix in-paragraph math spacing */
-  body { font-size: 1.1em }
+p > math[display="block"] {
-.katex { font-size: 1.0em !important }
+  margin-top: 1em;
+  margin-bottom: 1em;
+}
 
-/* Fix spacing of numbered equations */
+/* Fix numbered equations alignment */
-td .katex-display { margin: 0 }
+div[id^="eq"] table {
-
+  width: 100%;
-/* Fix equations width */
+  position: relative;
-div[id^="eq"] table { width: 100%; }
+}
+div[id^="eq"] td[style="text-align: center;"] {
+  width: 100%;
+}
+div[id^="eq"] td[style="text-align: right;"] {
+  position: absolute;
+  right: 0;
+  top: calc(50% - 0.6em);
+}
 
 /* Center the title */
 header { text-align: center }