diff --git a/README.rst b/README.rst
index 0b5d5f4eab859467fa0b032c2be29b35432f1e5f..a3f754413443fab303b93ba357e66384d956db6a 100644
--- a/README.rst
+++ b/README.rst
@@ -1,3 +1,4 @@
+================================
 Big Fluid and Particle Simulator
 ================================
 
@@ -20,6 +21,9 @@ I plan on adding documentation on the procedure as when other people
 show interest in using the code, or when time permits.
 
 
+.. _sec-installation:
+
+------------
 Installation
 ------------
 
@@ -54,6 +58,7 @@ Also, in order to run the C++ code you need to have an MPI compiler
 installed, the HDF5 C library as well as FFTW 3 (at least 3.3 I think).
 
 
+--------
 Comments
 --------
 
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index bc5a5e6ff3752b5de7e4443e6189f204eb0621e9..aa11892d21ecb1e2af658a3d7885c81e0d84d44b 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -699,6 +699,12 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             self,
             field_name = 'vorticity',
             iteration = 0):
+        """read the Fourier representation of a vector field.
+
+        Read the binary file containing iteration ``iteration`` of the
+        field ``field_name``, and return it as a properly shaped
+        ``numpy.memmap`` object.
+        """
         return np.memmap(
                 os.path.join(self.work_dir,
                              self.simname + '_{0}_i{1:0>5x}'.format('c' + field_name, iteration)),
diff --git a/bfps/tools.py b/bfps/tools.py
index 964c84329a080c6c7741f3ea56707f2fc13fee02..0ecfc768f6a972cfb5520b2998bc6eced110e46e 100644
--- a/bfps/tools.py
+++ b/bfps/tools.py
@@ -113,6 +113,7 @@ def padd_with_zeros(
     :param odtype: data type to use --- in principle conversion between
                   single and double precision can be performed with this
                   function as well.
+                  If ``None``, then use ``a.dtype``.
     :type n0: int
     :type n1: int
     :type n2: int
diff --git a/documentation/_static/api.rst b/documentation/_static/api.rst
index ac3f757b4d6e157c8e835fb2648030e73032059e..088df2ffa83ca3d6413e18c23c68542bdabd4a90 100644
--- a/documentation/_static/api.rst
+++ b/documentation/_static/api.rst
@@ -1,19 +1,23 @@
+===
 API
 ===
 
-bfps.NavierStokes
------------------
+----------
+bfps.tools
+----------
 
-.. autoclass:: NavierStokes.NavierStokes
+.. automodule:: tools
     :members:
     :undoc-members:
     :inherited-members:
     :show-inheritance:
 
-bfps.tools
-----------
 
-.. automodule:: tools
+-----------------
+bfps.NavierStokes
+-----------------
+
+.. autoclass:: NavierStokes.NavierStokes
     :members:
     :undoc-members:
     :inherited-members:
diff --git a/documentation/_static/development.rst b/documentation/_static/development.rst
index 20553e8da417b04b1919ff43db7095bfe55ca9f7..d2aa9da37c0822aa82b635c7684b8f990146e079 100644
--- a/documentation/_static/development.rst
+++ b/documentation/_static/development.rst
@@ -1,7 +1,9 @@
+===========
 Development
 ===========
 
 
+---------------------
 Versioning guidelines
 ---------------------
 
@@ -47,6 +49,7 @@ projects, these may change.
    forked from `vX.Y`, and then rule 3 is adapted accordingly.
 
 
+------------
 Code testing
 ------------
 
diff --git a/documentation/_static/overview.rst b/documentation/_static/overview.rst
index b694f002ff03d5072a92cff0f9c61fbaca8433cc..1d5d11cbffbd5e3933ff2815b63c73442b964808 100644
--- a/documentation/_static/overview.rst
+++ b/documentation/_static/overview.rst
@@ -1,6 +1,8 @@
+=====================
 Overview and Tutorial
 =====================
 
+---------
 Equations
 ---------
 
@@ -22,6 +24,126 @@ In fact, the code solves the vorticity formulation of these equations:
     \mathbf{\omega} \cdot \nabla \mathbf{u} +
     \nu \Delta \mathbf{\omega} + \nabla \times \mathbf{f}
 
+-----------
+Conventions
+-----------
+
+The C++ backend is based on ``FFTW``, therefore the Fourier
+representations are *transposed*.
+In brief, this is the way the fields are represented on disk and in
+memory (both in the C++ backend and in Python postprocessing):
+
+    * real space representations of 3D vector fields consist of
+      contiguous arrays, with the shape ``(nz, ny, nx, 3)``:
+      :math:`n_z \times n_y \times n_x` triplets, where :math:`z` is the
+      slowest coordinate, :math:`x` the fastest; each triplet is then
+      the sequence of :math:`x` component, :math:`y` component and
+      :math:`z` component.
+
+    * Fourier space representations of 3D vector fields consist of
+      contiguous arrays, with the shape ``(ny, nz, nx/2+1, 3)``:
+      :math:`k_y` is the slowest coordinate, :math:`k_x` the fastest;
+      each triplet of 3 complex numbers is then the :math:`(x, y, z)`
+      components, as ``FFTW`` requires for the correspondence with the
+      real space representations.
+
+:func:`read_cfield <NavierStokes.NavierStokes.read_cfield>` will return
+a properly shaped ``numpy.array`` containing a snapshot of the Fourier
+representation of a 3D field.
+
+If you'd like to construct the corresponding wave numbers, you can
+follow this procedure:
+
+.. code:: python
+
+    import numpy as np
+    from bfps import NavierStokes
+
+    c = NavierStokes(
+            work_dir = '/location/of/simulation/data',
+            simname = 'simulation_name_goes_here')
+    df = c.get_data_file()
+    kx = df['kspace/kx'].value
+    ky = df['kspace/ky'].value
+    kz = df['kspace/kz'].value
+    df.close()
+    kval = np.zeros(kz.shape + ky.shape + kx.shape + (3,),
+                    dtype = kx.dtype)
+    kval[..., 0] = kx[None, None, :]
+    kval[..., 1] = ky[:, None, None]
+    kval[..., 2] = kz[None, :, None]
+
+``kval`` will have the same shape as the result of
+:func:`read_cfield <NavierStokes.NavierStokes.read_cfield>`.
+Obviously, the machine being used should have enough RAM to hold the
+field...
+
+--------
 Tutorial
 --------
 
+First DNS
+---------
+
+Installing ``bfps`` is not trivial, and the instructions are in
+:ref:`sec-installation`.
+After installing, you should have a new executable script
+available, called ``bfps``, that you can execute.
+Just executing it will run a small test DNS on a real space grid of size
+:math:`32 \times 32 \times 32`, in the ``N0032`` folder in your current
+folder, with the simulation name ``test``.
+So, open a console, and type ``bfps``:
+
+.. code:: bash
+
+    # depending on how curious you are, you may have a look at the
+    # options first:
+    bfps --help
+    # or you may just run it:
+    bfps
+
+The simulation itself should not take more than a few seconds, since
+this is just a :math:`32^3` simulation run for 8 iterations.
+First thing you can do afterwards is open up a python console, and type
+the following:
+
+.. code:: python
+
+    import numpy as np
+    from bfps import NavierStokes
+
+    c = NavierStokes(
+            work_dir = '/location/of/simulation/data',
+            simname = 'simulation_name_goes_here')
+    c.compute_statistics()
+    print ('Rlambda = {0:.0f}, kMeta = {1:.4f}, CFL = {2:.4f}'.format(
+            c.statistics['Rlambda'],
+            c.statistics['kMeta'],
+            (c.parameters['dt']*c.statistics['vel_max'] /
+             (2*np.pi/c.parameters['nx']))))
+    print ('Tint = {0:.4e}, tauK = {1:.4e}'.format(c.statistics['Tint'],
+                                                   c.statistics['tauK']))
+    print ('total time simulated is = {0:.4e} Tint, {1:.4e} tauK'.format(
+            c.data_file['iteration'].value*c.parameters['dt'] / c.statistics['Tint'],
+            c.data_file['iteration'].value*c.parameters['dt'] / c.statistics['tauK']))
+
+:func:`compute_statistics <NavierStokes.NavierStokes.compute_statistics>`
+will read the data
+file generated by the DNS, compute a bunch of basic statistics, for
+example the Taylor scale Reynolds number :math:`R_\lambda` that we're
+printing in the example code.
+
+What happens is that the DNS will have generated an ``HDF5`` file
+containing a bunch of specific datasets (spectra, moments of real space
+representations, etc).
+The function
+:func:`compute_statistics <NavierStokes.NavierStokes.compute_statistics>`
+performs simple postprocessing that may however be expensive, therefore
+it also saves some data into a ``<simname>_postprocess.h5`` file, and
+then it also performs some time averages, yielding the ``statistics``
+dictionary that is used in the above code.
+
+What happened
+-------------
+
+Behind the scenes, something more complicated took place.