diff --git a/doc/Makefile b/doc/Makefile
index 08a0d326ad065ed1f248688a15974bc23939e097..e9190d0a0e8616f2f77552a7b47ec07fe39c85a6 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -13,11 +13,9 @@ help:
 
 .PHONY: help Makefile
 
-doxygen:
-	doxygen Doxyfile
-	breathe-apidoc -f -g struct,file doxygen_build/xml -o "$(SOURCEDIR)/Doxygen" -p CSDS
-
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 %: Makefile
+	doxygen Doxyfile
+	breathe-apidoc -f -g struct,file doxygen_build/xml -o "$(SOURCEDIR)/Doxygen" -p CSDS
 	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/doc/source/Examples/Energy.png b/doc/source/Examples/Energy.png
new file mode 100644
index 0000000000000000000000000000000000000000..f068bba105a3bff9fce4e3a75926f5b0a7aadedd
Binary files /dev/null and b/doc/source/Examples/Energy.png differ
diff --git a/doc/source/Examples/index.rst b/doc/source/Examples/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..aa773af475bc4014f2e2eef6dfacc332ba6cc87b
--- /dev/null
+++ b/doc/source/Examples/index.rst
@@ -0,0 +1,45 @@
+.. Examples
+   Loic Hausammann 2021
+
+Examples
+========
+
+Three different examples are available in the repository.
+The first one is described in :ref:`running`.
+Therefore I will only described the two others.
+
+Generating a Snapshot
+---------------------
+
+When running a deep analysis of a single given time, it might be faster
+(or more convenient) to use a snapshot.
+The file ``examples/create_snapshot.py`` can generate such files.
+It takes three parameters: the time (``--time``), the output file (``--output``) and
+the logfiles.
+Then it simply loops over all the existing fields thanks to ``get_list_fields`` and write
+them into an HDF5 file.
+
+
+Simple Orbital Simulation
+-------------------------
+
+The example in ``examples/SimpleOrbits`` consists in a
+single particle orbiting around a point mass (see ``README`` for more details).
+The parameters for the snapshots and the CSDS are set in order to have a large
+time resolution for the snapshots while a poor resolution for the CSDS.
+Even in such unfavorable case, the CSDS is still able to manage an accurate evolution
+of the kinetic energy (see figure below).
+The places where the CSDS line touches the snapshot line correspond to
+where the records are written.
+
+.. figure:: Energy.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Kinetic energy as function of time
+
+    This figure is produced by ``examples/SimpleOrbits`` and shows the evolution of
+    kinetic energy. As the energy is conserved, an oscillation is present due to the
+    conversion between kinetic and potential energy.
+    As it can be seen, in this example, the CSDS at low time resolution is as accurate
+    than the snapshots at high resolution.
diff --git a/doc/source/NewScheme/index.rst b/doc/source/NewScheme/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..95971eaac4ad9153907da7967ae215b52a61b497
--- /dev/null
+++ b/doc/source/NewScheme/index.rst
@@ -0,0 +1,7 @@
+.. NewScheme
+   Loic Hausammann 2021
+
+.. _implementation:
+
+Implementing new Scheme
+=======================
diff --git a/doc/source/Parameters/example.yml b/doc/source/Parameters/example.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b2c5827d3d2ce8ea53c8ae98979f42dfb55ca9c2
--- /dev/null
+++ b/doc/source/Parameters/example.yml
@@ -0,0 +1,65 @@
+Header:
+  BoxSize: [10.0401, 10.0401, 10.0401]  # Size of the box
+  Dimension: 3                          # Number of dimensions
+  RunName: Untitled SWIFT simulation    # Name of the simulation
+  Periodic: 0                           # Are we using periodic boundaries?
+  NumberParts: 0                        # Initial number of gas particles
+  NumberSParts: 0                       # Initial Number of stars
+  NumberGParts: 1                       # Initial Number of gravity particles
+
+Cosmology:
+  Omega_cdm: 0.2305    # Cosmological parameter for cold dark matter at z=0
+  Omega_lambda: 0.724  # Same for dark energy
+  Omega_b: 0.0455      # Same for baryons
+  Omega_r: 0           # Same for radiations
+  Omega_k: 0           # Same for curvature
+  Omega_nu_0: 0        # Same for neutrino
+  w_0: -1              # First parameter for the dark energy equation of state
+  w_a: 0               # Second parameter
+  Hubble0: 70.3        # The Hubble-Lemaitre constant at z=0 in internal units
+
+InternalUnitSystem:
+  UnitMass_in_cgs: 5.97e+27             # Unit of mass in gram
+  UnitLength_in_cgs: 1.49e+13           # Unit of length in cm
+  UnitTime_in_cgs: 3.14315e+07          # Unit of time in s
+  UnitCurrent_in_cgs: 1                 # Unit of current in cgs
+  UnitTemp_in_cgs: 1                    # Unit of temperature in cgs
+
+Code:
+  Code: SWIFT                 # Name of the code that ran the simulation
+  CodeVersion: 0.9.0          # Version of the code
+  CompilerName: GCC           # Name of the compiler
+  CompilerVersion: 9.3.0      # Version of the compiler
+  GitBranch: fix_gear         # Git branch of the compiled code
+  GitRevision: v0.9.0-501-ge8e3806f-dirty   # Git revision of the compiled code
+  GitDate: 2021-06-08 10:56:18 +0200        # Date of the compilation
+  ConfigurationOptions: '--with-python=/usr --enable-csds --with-ext-potential=point-mass' # Configuration options used
+  RandomSeed: 0               # Random seed used
+
+Policy:        # The policies corresponds to the runtime argument used
+  steal: 1
+  keep: 1
+  block: 0
+  cpu tight: 0
+  mpi: 0       # Did we use MPI?
+  numa affinity: 1
+  hydro: 0     # Did we use hydrodynamics?
+  self gravity: 0  # Did we use self gravity?
+  external gravity: 1  # Did we use external gravity?
+  cosmological integration: 0  # Did we use a cosmological integration?
+  drift everything: 0
+  reconstruct multi-poles: 0
+  temperature: 0
+  cooling: 0  # Did we include radiative cooling?
+  stars: 0    # Did we use stars (feedback not included)?
+  structure finding: 0
+  star formation: 0  # Did we use star formation?
+  feedback: 0        # Did we use stellar feedback?
+  black holes: 0     # Did we use black holes?
+  fof search: 0
+  time-step limiter: 0
+  time-step sync: 0
+  csds: 1      # Output written with the CSDS
+  line of sight: 0
+  sink: 0      # Did we use sink particles?
+  rt: 0        # Did we use radiative transfer?
diff --git a/doc/source/Parameters/index.rst b/doc/source/Parameters/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1737f1928e3e01a5524680e21c38130c8c50b3dc
--- /dev/null
+++ b/doc/source/Parameters/index.rst
@@ -0,0 +1,19 @@
+.. Parameters
+   Loic Hausammann 2021
+
+Parameters File
+===============
+
+Some information about the simulation are stored within a parameter file.
+The internal units are stored in this file.
+Among all the parameters, an important block is for the cosmological parameters that helps to
+do the interpolation in cosmological simulations.
+
+The initial number of particles is used for the generation of index files.
+The CSDS uses them as initial guess on the number of particles contained in the logfile.
+Without it, the CSDS would have to allocate a small initial array and increase it multiple times until
+reaching a size large enough to store all the particles.
+This allows to quickly obtain an array large enough.
+
+
+.. literalinclude:: example.yml
diff --git a/doc/source/QuickStart/compiling_code.rst b/doc/source/QuickStart/compiling_code.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3d78b60a44278005b8e78c3f0f25772cd5f374dd
--- /dev/null
+++ b/doc/source/QuickStart/compiling_code.rst
@@ -0,0 +1,31 @@
+.. CompilingCode
+   Loic Hausammann 2021
+
+.. _compiling_code:
+
+How to Compile the Code
+=======================
+
+To compile this code, you need to download `SWIFT's repository <www.swiftsim.com>`_
+and follow the instruction on how to compile it.
+Currently, the CSDS is implemented only within a few different scheme.
+If you wish to implement a new one, please follow the implementation
+within the other schems and :ref:`implementation`.
+During the configuration step, the CSDS needs to be enabled along with python wrapper if needed:
+
+.. code-block:: bash
+
+   ./autogen.sh
+   ./configure --enable-csds --with-python=$PYTHON_PATH
+
+On a standard ubuntu, ``$PYTHON_PATH`` is given by ``/usr/``.
+It is worth to mention that the CSDS will not work with python 2.7.x
+and that the user should move to 3.0 as soon as possible.
+
+Then both SWIFT and the CSDS can be compiled:
+
+.. code-block:: bash
+
+   make -j 10
+
+where 10 is the number of threads for the compilation.
diff --git a/doc/source/QuickStart/example.png b/doc/source/QuickStart/example.png
new file mode 100644
index 0000000000000000000000000000000000000000..8484fa711b34f8487e941e372965fd87f0e2700e
Binary files /dev/null and b/doc/source/QuickStart/example.png differ
diff --git a/doc/source/QuickStart/index.rst b/doc/source/QuickStart/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..ee70ca744a9dfab064241fd9982558079dcc300e
--- /dev/null
+++ b/doc/source/QuickStart/index.rst
@@ -0,0 +1,22 @@
+.. QuickStart
+   Loic Hausammann 2021
+
+.. _quick_start:
+
+QuickStart
+==========
+
+Here you will find how to compile the CSDS and run it on your first simulations.
+The CSDS is currently fully based on SWIFT and thus cannot be used directly without this code.
+Anyway, this code can still be used without running SWIFT at the condition that the output file
+correspond.
+In the future, we are planning to make this library fully independent from SWIFT.
+
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   compiling_code
+   running_first
+   using_python
diff --git a/doc/source/QuickStart/running_first.rst b/doc/source/QuickStart/running_first.rst
new file mode 100644
index 0000000000000000000000000000000000000000..043bfdd27a2979268b3e283b66c44c2cb776bf40
--- /dev/null
+++ b/doc/source/QuickStart/running_first.rst
@@ -0,0 +1,49 @@
+.. RunningFirst
+   Loic Hausammann 2021
+
+.. _running:
+
+
+Running Your First Example
+==========================
+
+The best approach to test the CSDS is to use an example within SWIFT.
+Let's start with the compilation of SWIFT and running our first example.
+
+
+.. code-block:: bash
+
+   ./autogen.sh
+   ./configure --enable-csds --with-python=$PYTHON_PATH
+   make -j 10
+   cd examples/HydroTests/SedovBlast_3D
+   ./run.sh
+
+The file ``run.sh`` is not running with the CSDS by default and needs to add
+the runtime parameter ``--csds`` to the command calling SWIFT
+(e.g. ``../../swift --csds --hydro --limiter --threads=4 sedov.yml``).
+If everything is fine, the last message from SWIFT should be ``Done``
+and two files have appeared (``index_0000.dump`` and ``index.yml``)
+Then the reader example can be used to interpret the output:
+
+.. code-block:: bash
+
+   cd ../../../csds/examples
+   python3 reader_example.py
+
+This example should give you the internal energy as function of the radius:
+
+.. figure:: example.png
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Internal energy as function of the radius for the Sedov Blast
+
+    This figure shows the internal energy as function of the radius for the
+    Sedov Blast. At such early time, all the energy is concentrated at the
+    center.
+
+This example can be used to see the evolution of the blast thanks to the parameter
+``--time`` and will read any logfile provided (by default the Sedov Blast example).
+This example is adapted only for hydrodynamics simulations and will raise an error
+on simulations without hydrodynamics.
diff --git a/doc/source/QuickStart/using_python.rst b/doc/source/QuickStart/using_python.rst
new file mode 100644
index 0000000000000000000000000000000000000000..2ad76a7764fbd4fa46879fc5b628b2b6ef7d9c33
--- /dev/null
+++ b/doc/source/QuickStart/using_python.rst
@@ -0,0 +1,108 @@
+.. PythonAPI
+   Loic Hausammann 2021
+
+.. _python:
+
+Python API
+==========
+
+The python API is relatively simple and fully self documented.
+Here, I will explain the main features and how they work.
+As the CSDS is currently depend on SWIFT and the scheme used,
+the library should be manually imported:
+
+.. code-block:: python
+
+   import sys
+   sys.path.append("SWIFT_PATH/csds/src/.libs/")
+   import libcsds as csds
+
+This module contains a single object (``Reader``) that fully deals with
+the CSDS files.
+To open a logfile, this object is call in a with statement:
+
+.. code-block:: python
+
+   with csds.Reader(filename, verbose=0, number_threads=1, number_index=10) as reader:
+       # Use the Reader
+
+In this example, all the parameters are set with their default values.
+The Reader is able to make some operations in parallel and
+will use the number of threads given by ``number_threads``.
+The index files are used to speedup the reading and are set at
+regular interval of simulation time (e.g. every :math:`(t_e - t_b) / (N - 1)`
+where :math:`t_b` (:math:`t_e`) is the initial (final) time and :math:`N` the number of index file.
+To get the best performances, the best is to write an index file as close
+as possible and below any future requested time.
+
+Now the Reader can provide some information about the simulation:
+
+.. code-block:: python
+
+   reader.get_list_fields(part_type=None)
+   reader.gas.get_list_fields()
+
+The two calls get the fields that exists for the particles.
+The only difference is that the second call will look at only the fields
+for the gas particles.
+The parameter ``part_type`` can be used to get the same result by simply providing 0.
+In SWIFT, this type corresponds to the gas.
+Multiple types of particles can be provided through a list.
+In such cases, the fields will be restricted to the common ones.
+We do not only provide ``reader.gas`` but also for other type of particles
+(``dark_matter``, ``dark_matter_background``, ``sinks``, ``stars`` and ``black_holes``).
+
+.. code-block:: python
+
+   reader.get_time_limits()
+
+This function returns the minimal and maximal time contained within the simulation.
+Any requested time should stay within the two limits (boundary included).
+
+.. code-block:: python
+
+   reader.get_box_size()
+
+This last method returns the box size of the simulation.
+
+Let's now focus on getting the particles. This is done by a single method called ``get_data``.
+Depending on the parameters, this function will have different optimizations:
+
+.. code-block:: python
+
+   reader.get_data(fields, time=0, filter_by_ids=None, part_type=None)
+   reader.gas.get_data(fields)
+
+It is mandatory to provide the required fields. If you wish to get all
+of them, simply provide the output of ``reader.get_list_fields``.
+```get_data`` returns a dictionary with the keys given by ``fields`` and
+the values are numpy arrays.
+The arrays are sorted according to the particle types.
+For the units, all the arrays are in internal units
+(see SWIFT for more details and be especially careful with cosmological simulations).
+Obviously, the parameter ``time`` is the requested time (in internal units)
+that should be within the limits provided by ``get_time_limits``.
+``filter_by_ids`` allows to pick only a subset of the particles.
+This parameter takes a list containing an array for each type of particle (or None).
+The list should have the same size than the number of particle types in SWIFT.
+The last parameter allows to read only some type of particles (e.g. ``part_type=0`` or
+``part_type=[0, 1]``).
+The two last parameters cannot be used together.
+
+When reading only a subset of the particles (``filter_by_ids``), the particles are searched
+within the index file and then read in the logfile.
+Due to the search, this parameter is especially efficient for a number
+of particles largely below the total number.
+It is worth to explicitly mention that the function will be inefficient
+if almost all the particles are requested.
+In such cases, it is better to request all the particles and then to drop the unwanted ones.
+
+When requesting all the particles (of a single type or many types),
+the code starts with the first particle in the index file and
+then read them the one after the onter.
+If some particles are removed from the simulation, the code stops with the current one
+and goes to the next one.
+While it is possible to find the removed particles with only the index file,
+it would take too much time to do it.
+The only case where anticipating the removal of particles would provide better results is when
+most of the particles are removed.
diff --git a/doc/source/conf.py b/doc/source/conf.py
index de7111a874eb40b3cda6ff9e78bec984191aca93..8333ec31c54a32c9ec4c3ec20f89c229bfce88d8 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -184,4 +184,9 @@ breathe_projects = {
     "CSDS":"../doxygen_build/xml/",
 }
 
+breathe_domain_by_extension = {
+    "h": "c",
+    "c": "c",
+}
+
 breathe_default_project = "CSDS"
diff --git a/doc/source/index.rst b/doc/source/index.rst
index c4424c9eb76b67f7cd481055bdf8f0fa2bf79a62..010650134b3af314f5d8ff0110d83dea23ccae71 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -6,9 +6,54 @@
 Welcome to Continuous Simulation Data Stream's documentation!
 =============================================================
 
+The Continuous Simulation Data Stream (CSDS) is a new output system
+that is especially designed for gravity based simulations.
+In such simulations, gravity is producing large differences of timescale between different area.
+An individual time step strategy has been designed a long time ago to speed up the simulations.
+Currently, the snapshots consists in writing the complete state of the simulation
+at regular interval.
+This approach does not take into account the timescales and thus produce large outputs.
+The CSDS has been designed to fully take into account the timescales to reduce the storage space.
+Our approach consists in using a single file describing the whole evolution of the simulation.
+In it, the particles/cells are written independently and according to their own time step.
+To reconstruct the state of the simulation at a given time, the particles are then
+simply interpolated from the nearest writing.
+Thanks to our approach, we do not only reduce the storage space, but also ensure to
+be able to recover the simulation at any time and not only a few discrete times.
+As the particles/cells are written at their own time scale, this approach ensures to have
+the same relative accuracy for all the particles.
+For more details on the technique, please look at `Hausammann, Gonnet and Schaller (2021)`.
+For simplicity, in the rest of the documentation, I will only mention the particles, but
+everything can be applied to cells.
+
+Our strategy for simulations requiring distributed memory (e.g. MPI) is
+to simply consider one output file per rank and thus each one can be considered
+independently from the others. It means that the arrays provided by the code
+needs to be concatenated together in order to get the full state.
+
+A few terms will be use within this documentation and need to be defined.
+The outputs from the CSDS are composed of three different (set of) files.
+The most important one is the logfile (``.dump``) and contains the whole simulation.
+Then comes the parameter file (``.yml``) that contains some information about the simulation
+(e.g. cosmological parameters, box size, ...).
+The last files are the index files (``.index``).
+They are generated by the reader and not during the simulation.
+They contain part of the information inside the logfile and
+are just a way to quickly access this large file.
+Within the logfile, everything is written in the form of records.
+Those small chunk of data can represent either a time step or a particle.
+
+
+The CSDS is licensed under the LGPL version 3.0.
+
 .. toctree::
    :maxdepth: 2
 
+   
+   QuickStart/index.rst
+   Examples/index.rst
+   Parameters/index.rst
+   NewScheme/index.rst
    Doxygen/index.rst
 
 Indices and tables
diff --git a/examples/SimpleOrbits/plotSolution.py b/examples/SimpleOrbits/plotSolution.py
index 9360b9fb6ffcd6e528cf7a5415855bdf6e5bf5ef..db4a8cceb25fa3df36378a17afa07eab4a3347e1 100644
--- a/examples/SimpleOrbits/plotSolution.py
+++ b/examples/SimpleOrbits/plotSolution.py
@@ -104,7 +104,7 @@ def doSnapshots():
 
         r = np.sum(pos**2, axis=1)**0.5
         v2 = np.sum(vel**2, axis=1)
-        E[i, :] = 0.5 * v2 - G * M / r
+        E[i, :] = 0.5 * v2
 
         # Get the pos / vel of the required particle
         ind = ids == id_focus
@@ -131,10 +131,11 @@ def doStatistics():
     """
     Do the plots with the energy output file.
     """
-    data = np.genfromtxt("energy.txt", names=True)
+    data = np.genfromtxt("statistics.txt")
 
-    times = data["Time"]
-    E = data["E_tot"]
+    print("You might need to modify the following variables")
+    times = data[:, 1]
+    E = data[:, 13]
     plt.figure(fig_1.number)
     plotRelative(times, E, "-", label="Statistics")
 
@@ -158,8 +159,10 @@ def doCSDS():
 
         for i, t in enumerate(times):
             # Get the next particles
-            pos, vel, ids = reader.get_particle_data(
+            out = reader.get_data(
                 ["Coordinates", "Velocities", "ParticleIDs"], t)
+            pos, vel = out["Coordinates"], out["Velocities"]
+            ids = out["ParticleIDs"]
             sort = np.argsort(ids)
             ids = ids[sort]
             rel_pos = pos[sort, :] - center
@@ -168,7 +171,7 @@ def doCSDS():
             # Compute the derived values
             r = np.sum(rel_pos**2, axis=1)**0.5
             v2 = np.sum(vel**2, axis=1)
-            E[i, :] = 0.5 * v2 - G * M / r
+            E[i, :] = 0.5 * v2
             ind = ids == id_focus
             p[i, :] = rel_pos[ind, :]
             v[i, :] = vel[ind, :]
@@ -196,7 +199,7 @@ doCSDS()
 # add text
 plt.figure(fig_1.number)
 plt.xlabel("Time [yr]")
-plt.ylabel(r"$\frac{E - E(t=0)}{E(t=0)}$")
+plt.ylabel(r"$\frac{E_{kin} - E_{kin}(t=0)}{E_{kin}(t=0)}$")
 plt.legend(ncol=2)
 plt.savefig("Energy.png")
 
diff --git a/examples/SimpleOrbits/simple_orbits.yml b/examples/SimpleOrbits/simple_orbits.yml
index 22fbffc2d37d539dbc2081dfb4d1ef148df35878..cf39708117a66057d93e039814d1dfa5bc9d2bb9 100644
--- a/examples/SimpleOrbits/simple_orbits.yml
+++ b/examples/SimpleOrbits/simple_orbits.yml
@@ -40,7 +40,7 @@ PointMassPotential:
 
 # Parameters governing the CSDS snapshot system
 CSDS:
-  delta_step:           4000     # Update the particle log every this many updates
+  delta_step:           1000     # Update the particle log every this many updates
   basename:             index  # Common part of the filenames
   initial_buffer_size:  0.01      # (Optional) Buffer size in GB
   buffer_scale:	        10     # (Optional) When buffer size is too small, update it with required memory times buffer_scale
diff --git a/src/csds_python_wrapper.c b/src/csds_python_wrapper.c
index c514fda7456922a2d422df6c9e8a35b17575c47e..97dbb4be2b4aa8d6bbe6bc72af8e35c12d95c4c4 100644
--- a/src/csds_python_wrapper.c
+++ b/src/csds_python_wrapper.c
@@ -270,7 +270,7 @@ csds_loader_create_output(void **output, const int *field_indices,
 
   struct csds_python_field python_fields[100];
 
-  /* Create the python list */
+  /* Create the python dictionary */
   PyObject *dict = PyDict_New();
   struct csds_python_field *current_field;
 
@@ -1006,7 +1006,7 @@ static PyTypeObject PyObjectReader_Type = {
         "-------\n\n"
         "get_time_limits\n"
         "    Provides the initial and final time of the simulation\n"
-        "get_particle_data\n"
+        "get_data\n"
         "    Reads some particle's data from the logfile\n\n"
         "Examples\n"
         "--------\n\n"
@@ -1017,7 +1017,7 @@ static PyTypeObject PyObjectReader_Type = {
         ">>>    if \"Coordinates\" not in fields:\n"
         ">>>        raise Exception(\"Field Coordinates not present in the "
         "logfile.\")\n"
-        ">>>    pos, ent = reader.get_particle_data([\"Coordinates\", "
+        ">>>    pos, ent = reader.get_data([\"Coordinates\", "
         "\"Entropies\"]"
         ", 0.5 * (t0 + t1))\n",
     .tp_init = (initproc)Reader_init,