diff --git a/doc/conf.py b/doc/conf.py
new file mode 100644
index 0000000000000000000000000000000000000000..e918f68fff53c12bfaedcd566e3b1e04a81f85c6
--- /dev/null
+++ b/doc/conf.py
@@ -0,0 +1,11 @@
+project = 'usadelndsoc'
+copyright = '2024, Pauli Virtanen'
+author = 'Pauli Virtanen'
+release = '0.1'
+
+extensions = []
+#templates_path = ['_templates']
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+html_theme = 'alabaster'
+#html_static_path = ['_static']
diff --git a/doc/discretization.md b/doc/discretization.md
deleted file mode 100644
index cee8b0d4c3e47002a893e3fbd423186aa0dae6bc..0000000000000000000000000000000000000000
--- a/doc/discretization.md
+++ /dev/null
@@ -1,249 +0,0 @@
-# Gauge-invariant discretization
-
-We do the discretization of the Usadel equation as follows. We first
-discretize the action $`S`$. The discrete saddle-point equation or expressions
-for currents are not derived analytically, but instead obtained from the computer
-implementation of the action via automatic differentiation (AD).
-This simplifies things somewhat, because proper handling of the gauge fields is
-simpler to do in the action formulation.
-
-The action is
-```math
-  S[Q] = \frac{i\pi\nu}{8} \mathrm{Tr}[D(\hat{\nabla}Q)^2 + \eta F_{ij}Q\hat{\nabla}_iQ\hat{\nabla}_jQ + 4 i \Omega Q]
-  \,,
-```
-where $`\hat{\nabla}_j Q = \partial_j Q - i[\mathcal{A}_j, Q]`$
-and $`\Omega = \epsilon\tau_3 + i\Delta`$, $`\Delta^\dagger = \Delta`$.
-We'd like the discretized functional to be also gauge invariant.
-
-We subdivide the space to cells, centered on a rectangular lattice
-$`\vec{r}_j=(h j_x, h j_y, h j_z)`$ with $`j_{x,y,z}\in\mathbb{Z}`$
-and *h* are the lattice spacings. See the figure below:
-
-![FIG 1. Lattice discretization](img/discretization.svg "FIG 1. Lattice discretization.")
-
-We choose $`Q(j) = Q(\vec{r}_j)`$ to be values of $`Q`$ at the lattice
-sites. We define the Wilson link matrices $`U_{ij} =U_{ji}^{-1}
-=U(\vec{r}_{i},\vec{r}_{j}) =
-\mathrm{Pexp}[i\int_{L(\vec{r}_i,\vec{r}_j)}d\vec{r}'\,\cdot\vec{\mathcal{A}}(\vec{r}')]`$
-where $`L(\vec{r}_i,\vec{r}_j)`$ is straight line from $`\vec{r}_j`$
-to $`\vec{r}_i`$ and $`\mathrm{Pexp}`$ is the path-ordered integral.
-Since the gauge field $`\vec{\mathcal{A}}`$ is fixed, the link
-matrices can be computed ahead of time.
-
-We want a discretized theory that is invariant under the discrete gauge transformation
-```math
-Q(j) \mapsto u(j) Q u(j)^{-1} \,,
-```
-and
-```math
-  U_{ij}
-  \mapsto
-  u(i)U_{ij}u(j)^{-1}
-  \,,
-```
-and in the continuum limit reduces to the original theory.  Such
-discretization has the advantage of being insensitive to the gauge
-choice, and having exact discrete conservation laws associated with
-the symmetry.
-
-## Derivative term
-
-We discretize the derivatives a
-```math
-  D_{ij}
-  :=
-  \frac{1}{h}[Q(i) U_{ij} - U_{ij} Q(j)]
-  \,,
-  \qquad
-  \Rightarrow
-  Q(i) D_{ij} = - D_{ij} Q(j)
-  \,,
-```
-which satisfies an anticommutation relation analogous to continuum derivative.
-
-We then discretize
-```math
-S_2
-=
-\mathrm{Tr} D(\hat{\nabla} Q)^2
-\mapsto
--D \frac{2d}{|\mathrm{neigh}|} h^d \sum_{\mathrm{neigh}(i,j)}\mathrm{Tr} D_{ij} D_{ji}
-```
-where $`d`$ is the space dimension and `|neigh|` the number of neighbors on the grid.
-
-Expanding in small $`h`$ we have
-```math
-S_2
-\simeq
-\frac{h^{d}d}{|\mathrm{neigh}|}
-\sum_{i}\sum_{j\in{}\mathrm{neigh}(i)}\mathrm{Tr}\Bigl(\frac{Q(i) - Q(j)}{h} - [i\mathcal{A}(\frac{\vec{r}_i+\vec{r}_j}{2}), Q(j)]\Bigr)^2
-\,,
-```
-which verifies the reduction to the continuum limit.
-
-## Field strength
-
-The link matrices have the usual continuum limit expansion
-```math
-U_{ji}
-=
-\sum_{n=0}^\infty i^n |\vec{r}_i-\vec{r}_j|^n \int_{-1/2}^{1/2}ds_1\int_{-1/2}^{s_1}ds_2\ldots\int_{-1/2}^{s_{n-1}}ds_n\, \mathcal{A}(s_1)\cdots\mathcal{A}(s_n)
-```
-where $`\mathcal{A}(s)=\mathcal{A}[(\frac{1}{2} - s)\vec{r}_i + (\frac{1}{2} + s)\vec{r}_j]`$.
-
-As well known in lattice QCD, the field strength can be expressed in
-terms of the link matrices $`U`$ via a plaquette loop.  Consider then
-the plaquette (see Fig. 1), with edges $`j\leftarrow{}i`$ and $`k\leftarrow{}l`$ in
-direction $`\mu`$ and $`l\leftarrow{}i`$, $`k\leftarrow{}j`$ in direction $`\nu`$, and expand around
-$`\vec{r}_0=\frac{\vec{r}_i+\vec{r}_j+\vec{r}_k+\vec{r}_l}{4}`$:
-```math
-  P_{lkji}
-  =
-  U_{il}U_{lk}U_{kj}U_{ji}
-  \,,
-```
-```math
-  \Rightarrow
-  P_{lkji}
-  \simeq
-  1
-  - i h^2
-  \bigl(
-  \partial_\mu \mathcal{A}_\nu - \partial_\nu \mathcal{A}_\mu - i[\mathcal{A}_\mu,\mathcal{A}_\nu]
-  \bigr)
-  +
-  \mathcal{O}(h^3)
-  =
-  1
-  + i h^2 F_{\mu \nu}(\vec{r}_0)
-  +
-  \mathcal{O}(h^3)
-```
-which then allows expressing $`F_{ij}`$ in terms of the link matrices.
-
-## Hall term
-
-In the plaquette of Fig. 1 we have
-```math
-  \mathrm{tr}
-  F_{\mu\nu}Q \hat{\nabla}_\mu Q \hat{\nabla}_\nu Q
-  \rvert_{\mu=ji, \nu=li}
-  \simeq
-  \frac{i}{h^4}
-  \mathrm{tr}
-  (1 - P_{lkji}) U_{ij} D_{ji} Q(i) D_{il} U_{li}
-```
-```math
-  =
-  \frac{-i}{h^4}
-  \mathrm{tr}
-  (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
-```
-So the structure is $`\frac{-i}{h^4}`$ $`\times`$ (gauge factors for
-counterclockwise loop $`-`$ return back clockwise) $`\times`$ $`DQD`$ for
-three sites in counterclockwise order.  Since only $`Q(i)`$, $`Q(j)`$,
-$`Q(l)`$ appear here, one can also think of this as a plaquette corner
-with $`Q(i)`$ being the corner.
-
-To do the sum over the $`\mu`$, $`\nu`$ indices, one can average over
-corners of the plaquettes as follows, which gives the final discretization
-of the Hall term:
-```math
-  S_H
-  =
-  \eta
-  \mathrm{Tr} F_{\mu\nu}Q\hat{\nabla}_\mu Q\hat{\nabla}_\nu Q
-```
-```math
-  S_H
-  \mapsto
-  \eta
-  \frac{h^d d}{4ih^4}
-  \sum_{\mathrm{plaqc}(ijkl; i)}
-  \mathrm{Tr}
-  \Bigl(
-  (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
-  \Bigr)
-```
-```math
-  S_H
-  =
-  \frac{\eta h^{d-4} d}{4i}
-  \sum_{\mathrm{plaqc}(i; jkl)}
-  \mathrm{Tr}
-  \Bigl(
-  (P_{lkji} - 1)(Q(i) - U_{ij} Q(j) U_{ji} Q(i) U_{il} Q(l) U_{li})
-  \Bigr)
-  \,,
-```
-where $`\mathrm{plaqc}(i;jkl)`$ are the plaquette corners. The last line
-follows with direct calculation, and considering the cyclic
-permutations.
-
-## Matrix current
-
-We can also calculate the matrix current $`j\mapsto{}i`$ between neighboring cells:
-```math
-  J_{ij}^a
-  =
-  \frac{\partial}{\partial\xi}
-  S[
-    U_{ij}
-    \mapsto
-    U_{ij}
-    +
-    T^a U_{ij}\xi
-    \,,
-    \;
-    U_{ji}
-    \mapsto
-    U_{ji}
-    -
-    U_{ji} T^a \xi
-  ]\rvert_{\xi=0}
-```
-where $`T^a`$ is an appropriate matrix generator, and the transformation
-is inserted to only one of the links.
-
-The local gauge invariance now implies that there is an exact discrete
-continuity equation
-```math
-  \sum_{j\in{}N_i} J_{ij}^a
-  =
-  \frac{\partial}{\partial\xi}
-  S[Q(i)\mapsto{}e^{\xi T^a}Q(i)e^{-\xi T^a}]\rvert_{\xi=0}
-  =
-  \frac{\partial}{\partial\xi}
-  S[\Omega(i)\mapsto{}e^{-\xi T^a}\Omega(i)e^{\xi T^a}]\rvert_{\xi=0}
-  =
-  R_{i}^a
-  \,.
-```
-This is one main advantage of the procedure, in addition to its
-invariance vs. gauge fixing.
-
-However, the definition of the matrix current is *asymmetric*, so
-that $`J_{ij}\ne{}-J_{ji}`$ as they differ by a quantity of order
-$`\mathcal{O}(h^2)`$. The asymmetry arises when $`T^a`$ does not commute
-with $`U_{ij}`$, so it is an issue only for the SU(2) component. The
-problem is only in the interpretation.
-
-The currents measured at sites *i* and *j* are different: since
-the parallel transport of spin between these two locations can imply
-rotation due to the gauge field, there's no reason why we should have
-$`J_{ij}=-J_{ji}`$, except in the continuum limit where the two points
-become close to each other.
-
-
-## Implementation
-
-The above is sufficient for numerical implementation: it is not
-necessary to find out the saddle point equation symbolically. It can
-be deduced by parametrizing first so that $`Q^2=1`$ and then using
-automatic differentiation (AD). Only a routine evaluating the value of
-the action is then needed, in terms of the AD dual variables.
-
-When looking for the saddle point in equilibrium, the problem becomes
-even simpler; *Q* has then additional restrictions, and we are looking
-for the minimum of the free energy *S*.
diff --git a/doc/discretization.rst b/doc/discretization.rst
new file mode 100644
index 0000000000000000000000000000000000000000..2b5a5ef68b926ba2217b9c4fe20492a8d81a940c
--- /dev/null
+++ b/doc/discretization.rst
@@ -0,0 +1,293 @@
+Gauge-invariant discretization
+==============================
+
+We do the discretization of the Usadel equation as follows. We first
+discretize the action :math:`S`. The discrete saddle-point equation or
+expressions for currents are not derived analytically, but instead
+obtained from the computer implementation of the action via automatic
+differentiation (AD). This simplifies things somewhat, because proper
+handling of the gauge fields is simpler to do in the action formulation.
+
+The action is
+
+.. math::
+
+     S[Q] = \frac{i\pi\nu}{8} \mathrm{Tr}[D(\hat{\nabla}Q)^2 + \eta F_{ij}Q\hat{\nabla}_iQ\hat{\nabla}_jQ + 4 i \Omega Q]
+     \,,
+
+where :math:`\hat{\nabla}_j Q = \partial_j Q - i[\mathcal{A}_j, Q]` and
+:math:`\Omega = \epsilon\tau_3 + i\Delta`,
+:math:`\Delta^\dagger = \Delta`. We’d like the discretized functional to
+be also gauge invariant.
+
+We subdivide the space to cells, centered on a rectangular lattice
+:math:`\vec{r}_j=(h j_x, h j_y, h j_z)` with
+:math:`j_{x,y,z}\in\mathbb{Z}` and *h* are the lattice spacings. See the
+figure below:
+
+.. figure:: img/discretization.svg
+   :alt: FIG 1. Lattice discretization.
+
+   FIG 1. Lattice discretization
+
+We choose :math:`Q(j) = Q(\vec{r}_j)` to be values of :math:`Q` at the
+lattice sites. We define the Wilson link matrices
+:math:`U_{ij} =U_{ji}^{-1} =U(\vec{r}_{i},\vec{r}_{j}) = \mathrm{Pexp}[i\int_{L(\vec{r}_i,\vec{r}_j)}d\vec{r}'\,\cdot\vec{\mathcal{A}}(\vec{r}')]`
+where :math:`L(\vec{r}_i,\vec{r}_j)` is straight line from
+:math:`\vec{r}_j` to :math:`\vec{r}_i` and :math:`\mathrm{Pexp}` is the
+path-ordered integral. Since the gauge field :math:`\vec{\mathcal{A}}`
+is fixed, the link matrices can be computed ahead of time.
+
+We want a discretized theory that is invariant under the discrete gauge
+transformation
+
+.. math:: Q(j) \mapsto u(j) Q u(j)^{-1} \,,
+
+and
+
+.. math::
+
+     U_{ij}
+     \mapsto
+     u(i)U_{ij}u(j)^{-1}
+     \,,
+
+and in the continuum limit reduces to the original theory. Such
+discretization has the advantage of being insensitive to the gauge
+choice, and having exact discrete conservation laws associated with the
+symmetry.
+
+Derivative term
+---------------
+
+We discretize the derivatives a
+
+.. math::
+
+     D_{ij}
+     :=
+     \frac{1}{h}[Q(i) U_{ij} - U_{ij} Q(j)]
+     \,,
+     \qquad
+     \Rightarrow
+     Q(i) D_{ij} = - D_{ij} Q(j)
+     \,,
+
+which satisfies an anticommutation relation analogous to continuum
+derivative.
+
+We then discretize
+
+.. math::
+
+   S_2
+   =
+   \mathrm{Tr} D(\hat{\nabla} Q)^2
+   \mapsto
+   -D \frac{2d}{|\mathrm{neigh}|} h^d \sum_{\mathrm{neigh}(i,j)}\mathrm{Tr} D_{ij} D_{ji}
+
+where :math:`d` is the space dimension and ``|neigh|`` the number of
+neighbors on the grid.
+
+Expanding in small :math:`h` we have
+
+.. math::
+
+   S_2
+   \simeq
+   \frac{h^{d}d}{|\mathrm{neigh}|}
+   \sum_{i}\sum_{j\in{}\mathrm{neigh}(i)}\mathrm{Tr}\Bigl(\frac{Q(i) - Q(j)}{h} - [i\mathcal{A}(\frac{\vec{r}_i+\vec{r}_j}{2}), Q(j)]\Bigr)^2
+   \,,
+
+which verifies the reduction to the continuum limit.
+
+Field strength
+--------------
+
+The link matrices have the usual continuum limit expansion
+
+.. math::
+
+   U_{ji}
+   =
+   \sum_{n=0}^\infty i^n |\vec{r}_i-\vec{r}_j|^n \int_{-1/2}^{1/2}ds_1\int_{-1/2}^{s_1}ds_2\ldots\int_{-1/2}^{s_{n-1}}ds_n\, \mathcal{A}(s_1)\cdots\mathcal{A}(s_n)
+
+where
+:math:`\mathcal{A}(s)=\mathcal{A}[(\frac{1}{2} - s)\vec{r}_i + (\frac{1}{2} + s)\vec{r}_j]`.
+
+As well known in lattice QCD, the field strength can be expressed in
+terms of the link matrices :math:`U` via a plaquette loop. Consider then
+the plaquette (see Fig. 1), with edges :math:`j\leftarrow{}i` and
+:math:`k\leftarrow{}l` in direction :math:`\mu` and
+:math:`l\leftarrow{}i`, :math:`k\leftarrow{}j` in direction :math:`\nu`,
+and expand around
+:math:`\vec{r}_0=\frac{\vec{r}_i+\vec{r}_j+\vec{r}_k+\vec{r}_l}{4}`:
+
+.. math::
+
+     P_{lkji}
+     =
+     U_{il}U_{lk}U_{kj}U_{ji}
+     \,,
+
+.. math::
+
+     \Rightarrow
+     P_{lkji}
+     \simeq
+     1
+     - i h^2
+     \bigl(
+     \partial_\mu \mathcal{A}_\nu - \partial_\nu \mathcal{A}_\mu - i[\mathcal{A}_\mu,\mathcal{A}_\nu]
+     \bigr)
+     +
+     \mathcal{O}(h^3)
+     =
+     1
+     + i h^2 F_{\mu \nu}(\vec{r}_0)
+     +
+     \mathcal{O}(h^3)
+
+which then allows expressing :math:`F_{ij}` in terms of the link
+matrices.
+
+Hall term
+---------
+
+In the plaquette of Fig. 1 we have
+
+.. math::
+
+     \mathrm{tr}
+     F_{\mu\nu}Q \hat{\nabla}_\mu Q \hat{\nabla}_\nu Q
+     \rvert_{\mu=ji, \nu=li}
+     \simeq
+     \frac{i}{h^4}
+     \mathrm{tr}
+     (1 - P_{lkji}) U_{ij} D_{ji} Q(i) D_{il} U_{li}
+
+.. math::
+
+     =
+     \frac{-i}{h^4}
+     \mathrm{tr}
+     (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
+
+So the structure is :math:`\frac{-i}{h^4}` :math:`\times` (gauge factors
+for counterclockwise loop :math:`-` return back clockwise)
+:math:`\times` :math:`DQD` for three sites in counterclockwise order.
+Since only :math:`Q(i)`, :math:`Q(j)`, :math:`Q(l)` appear here, one can
+also think of this as a plaquette corner with :math:`Q(i)` being the
+corner.
+
+To do the sum over the :math:`\mu`, :math:`\nu` indices, one can average
+over corners of the plaquettes as follows, which gives the final
+discretization of the Hall term:
+
+.. math::
+
+     S_H
+     =
+     \eta
+     \mathrm{Tr} F_{\mu\nu}Q\hat{\nabla}_\mu Q\hat{\nabla}_\nu Q
+
+.. math::
+
+     S_H
+     \mapsto
+     \eta
+     \frac{h^d d}{4ih^4}
+     \sum_{\mathrm{plaqc}(ijkl; i)}
+     \mathrm{Tr}
+     \Bigl(
+     (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
+     \Bigr)
+
+.. math::
+
+     S_H
+     =
+     \frac{\eta h^{d-4} d}{4i}
+     \sum_{\mathrm{plaqc}(i; jkl)}
+     \mathrm{Tr}
+     \Bigl(
+     (P_{lkji} - 1)(Q(i) - U_{ij} Q(j) U_{ji} Q(i) U_{il} Q(l) U_{li})
+     \Bigr)
+     \,,
+
+where :math:`\mathrm{plaqc}(i;jkl)` are the plaquette corners. The last
+line follows with direct calculation, and considering the cyclic
+permutations.
+
+Matrix current
+--------------
+
+We can also calculate the matrix current :math:`j\mapsto{}i` between
+neighboring cells:
+
+.. math::
+
+     J_{ij}^a
+     =
+     \frac{\partial}{\partial\xi}
+     S[
+       U_{ij}
+       \mapsto
+       U_{ij}
+       +
+       T^a U_{ij}\xi
+       \,,
+       \;
+       U_{ji}
+       \mapsto
+       U_{ji}
+       -
+       U_{ji} T^a \xi
+     ]\rvert_{\xi=0}
+
+where :math:`T^a` is an appropriate matrix generator, and the
+transformation is inserted to only one of the links.
+
+The local gauge invariance now implies that there is an exact discrete
+continuity equation
+
+.. math::
+
+     \sum_{j\in{}N_i} J_{ij}^a
+     =
+     \frac{\partial}{\partial\xi}
+     S[Q(i)\mapsto{}e^{\xi T^a}Q(i)e^{-\xi T^a}]\rvert_{\xi=0}
+     =
+     \frac{\partial}{\partial\xi}
+     S[\Omega(i)\mapsto{}e^{-\xi T^a}\Omega(i)e^{\xi T^a}]\rvert_{\xi=0}
+     =
+     R_{i}^a
+     \,.
+
+This is one main advantage of the procedure, in addition to its
+invariance vs. gauge fixing.
+
+However, the definition of the matrix current is *asymmetric*, so that
+:math:`J_{ij}\ne{}-J_{ji}` as they differ by a quantity of order
+:math:`\mathcal{O}(h^2)`. The asymmetry arises when :math:`T^a` does not
+commute with :math:`U_{ij}`, so it is an issue only for the SU(2)
+component. The problem is only in the interpretation.
+
+The currents measured at sites *i* and *j* are different: since the
+parallel transport of spin between these two locations can imply
+rotation due to the gauge field, there’s no reason why we should have
+:math:`J_{ij}=-J_{ji}`, except in the continuum limit where the two
+points become close to each other.
+
+Implementation
+--------------
+
+The above is sufficient for numerical implementation: it is not
+necessary to find out the saddle point equation symbolically. It can be
+deduced by parametrizing first so that :math:`Q^2=1` and then using
+automatic differentiation (AD). Only a routine evaluating the value of
+the action is then needed, in terms of the AD dual variables.
+
+When looking for the saddle point in equilibrium, the problem becomes
+even simpler; *Q* has then additional restrictions, and we are looking
+for the minimum of the free energy *S*.
diff --git a/doc/index.rst b/doc/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..05258158e0095e09c65add082a00fc9bddcd96d9
--- /dev/null
+++ b/doc/index.rst
@@ -0,0 +1,20 @@
+.. usadelndsoc documentation master file, created by
+   sphinx-quickstart on Tue Apr 23 13:16:43 2024.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+usadelndsoc
+===========
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   discretization
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/doc/pandoc/convert.sh b/doc/pandoc/convert.sh
deleted file mode 100755
index 1c4608d37004105df39508e0af2a3324668ff34d..0000000000000000000000000000000000000000
--- a/doc/pandoc/convert.sh
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/bin/bash
-PDIR=$(dirname "$0")
-pandoc --from markdown --to html \
-       --embed-resources --standalone \
-       --lua-filter $PDIR/gitlab-math.lua \
-       --lua-filter $PDIR/fix-links.lua \
-       --katex \
-       -o "${1/.md/.html}" "$1"
diff --git a/doc/pandoc/fix-links.lua b/doc/pandoc/fix-links.lua
deleted file mode 100644
index 36ac7dc1274d995f103ab5649a6f48a8ab91def0..0000000000000000000000000000000000000000
--- a/doc/pandoc/fix-links.lua
+++ /dev/null
@@ -1,44 +0,0 @@
--- see if the file exists
-function file_exists(file)
-  local f = io.open(file, "rb")
-  if f then f:close() end
-  return f ~= nil
-end
-
--- 1. Replaces .md with .html
--- 2. Replaces absolute paths with relative ones
--- 3. Replaces folder links with links to their index / readme
-function fix_link (url) 
-  -- Replace md with html
-  fixed_url = url:gsub("%.md", ".html")
-  -- Replace project-root-relative (i.e. /path/to/thing.md) paths with relative ones
-  fixed_url = fixed_url:gsub("^/", (
-    function(s)
-      basedir = io.popen("realpath --relative-to=\"$(dirname \"" .. s:sub(2) .. "\")\" .", 'r'):read('*a')
-      return basedir:gsub("%s", "") .. "/" .. s:sub(2)
-    end
-  ))
-  -- Allowed readme / index names
-  dir_name = {} -- values are functions that return the path to the new link
-  dir_name["index.md"] = (function(s) return s:gsub("/?$", "/index.html") end)
-  dir_name["INDEX.md"] = (function(s) return s:gsub("/?$", "/INDEX.html") end)
-  dir_name["readme.md"] = (function(s) return s:gsub("/?$", "/readme.html") end)
-  dir_name["README.md"] = (function(s) return s:gsub("/?$", "/README.html") end)
-  -- Is the url pointing to a folder?
-  file_or_folder = io.popen("[ -d '" .. fixed_url .. "' ] && echo 'folder'", 'r'):read('*a')
-  -- If it's a folder, find the readme / index and replace the link with that. 
-  if file_or_folder:gsub("%s", "") == "folder" then
-    for index,mk_path in pairs(dir_name) do
-      if file_exists(fixed_url .. "/" .. index) then
-        fixed_url = mk_path(fixed_url)
-        break
-      end
-    end
-  end
-  return fixed_url
-end
-
-
-function Link(link) link.target = fix_link(link.target); return link end
-function Image(img) img.src = fix_link(img.src); return src end
--- return {{Meta = Meta}, {Link = Link, Image = Image}}
\ No newline at end of file
diff --git a/doc/pandoc/gitlab-math.lua b/doc/pandoc/gitlab-math.lua
deleted file mode 100644
index 616103117e2a236ec9933297b4d244acc07e40c4..0000000000000000000000000000000000000000
--- a/doc/pandoc/gitlab-math.lua
+++ /dev/null
@@ -1,18 +0,0 @@
-function Math(el)
-    if el.mathtype == "InlineMath" then
-        if el.text:sub(1,1) == '`' and el.text:sub(#el.text) == '`' then
-            local text = el.text:sub(2,#el.text-1)
-            return pandoc.Math(el.mathtype, text)
-        else
-            local cont = pandoc.read(el.text)
-            return el
-            --return { pandoc.Str("$") } .. cont.blocks[1].content .. { pandoc.Str("$") }
-        end
-    end
-end
-
-function CodeBlock(el)
-    if el.classes[1] == "math" then
-        return pandoc.Para({ pandoc.Math("DisplayMath", el.text) })
-    end
-end
diff --git a/meson.build b/meson.build
index ddd2aefc199c69fd7243e69441d001bc52d76af7..bafe4e1221a4a22df6a38d05af85ac3194294dd1 100644
--- a/meson.build
+++ b/meson.build
@@ -41,7 +41,6 @@ incdir_numpy = run_command(py3,
   ],
   check: true
 ).stdout().strip()
-cc = meson.get_compiler('cpp')
 numpy_dep = declare_dependency(
   compile_args : ['-DNPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION'],
   include_directories : [include_directories(incdir_numpy)],
@@ -66,3 +65,4 @@ deps = [cppad_dep, eigen3_dep]
 #executable('main2', ['src/xtest.f95', 'src/adjac/adjac.f95'])
 
 subdir('usadelndsoc')
+subdir('doc')
diff --git a/usadelndsoc/meson.build b/usadelndsoc/meson.build
index 3da9aeddd24f4b82e8996f38dd07c8933a26937e..21fc6d798e7b2811242c8fadabfb3718f92c3aee 100644
--- a/usadelndsoc/meson.build
+++ b/usadelndsoc/meson.build
@@ -1,4 +1,4 @@
-py3.extension_module(
+usadelndsoc_pymod = py3.extension_module(
   '_core',
   ['../src/core.cpp'],
   dependencies : deps + [numpy_dep, py3_dep, pybind11_dep],