diff --git a/src/engine.c b/src/engine.c
index 84e5655dd53e42b75104c9b68744ede74ac378c3..29d35e087bc82aaae3ee392428cb05d420be288b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2724,7 +2724,7 @@ void engine_step(struct engine *e) {
     mask |= 1 << task_type_sub_pair;
     mask |= 1 << task_type_ghost;
     mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
-                                                      does not harm */
+                                                       does not harm */
 
     submask |= 1 << task_subtype_density;
     submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 258115f48a8c5939599ba3271e19c0d88c8f8b89..b70ed690d368cfa797c513ff724ffddb2d0210df 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -72,6 +72,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  * @brief Prepares a particle for the volume calculation.
  *
  * Simply makes sure all necessary variables are initialized to zero.
+ * Initializes the Voronoi cell.
  *
  * @param p The particle to act upon
  */
@@ -88,19 +89,13 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
 /**
  * @brief Finishes the volume calculation.
  *
- * Multiplies the density and number of neighbours by the appropiate constants
- * and adds the self-contribution term. Calculates the volume and uses it to
- * update the primitive variables (based on the conserved variables). The latter
- * should only be done for active particles. This is okay, since this method is
- * only called for active particles.
+ * Calls the finalize method on the Voronoi cell, which calculates the volume
+ * and centroid of the cell. We use the return value of this function to set
+ * a new value for the smoothing length and possibly force another iteration
+ * of the volume calculation for this particle. We then use the volume to
+ * convert conserved variables into primitive variables.
  *
- * Multiplies the components of the matrix E with the appropriate constants and
- * inverts it. Initializes the variables used during the gradient loop. This
- * cannot be done in hydro_prepare_force, since that method is called for all
- * particles, and not just the active ones. If we would initialize the
- * variables there, gradients for passive particles would be zero, while we
- * actually use the old gradients in the flux calculation between active and
- * passive particles.
+ * This method also initializes the gradient variables (if gradients are used).
  *
  * @param p The particle to act upon.
  * @param The current physical time.
@@ -200,8 +195,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
 /**
  * @brief Reset acceleration fields of a particle
  *
- * This is actually not necessary for GIZMO, since we just set the accelerations
- * after the flux calculation.
+ * This is actually not necessary for Shadowswift, since we just set the
+ * accelerations after the flux calculation.
  *
  * @param p The particle to act upon.
  */
@@ -256,7 +251,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
 /**
  * @brief Extra operations to be done during the drift
  *
- * Not used for GIZMO.
+ * Not used for Shadowswift.
  *
  * @param p Particle to act upon.
  * @param xp The extended particle data to act upon.
@@ -272,7 +267,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
  *
  * We use the new conserved variables to calculate the new velocity of the
  * particle, and use that to derive the change of the velocity over the particle
- * time step.
+ * time step. We also add a correction to the velocity which steers the
+ * generator towards the centroid of the cell.
  *
  * If the particle time step is zero, we set the accelerations to zero. This
  * should only happen at the start of the simulation.
@@ -356,7 +352,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
 /**
  * @brief Extra operations done during the kick
  *
- * Not used for GIZMO.
+ * Not used for Shadowswift.
  *
  * @param p Particle to act upon.
  * @param xp Extended particle data to act upon.
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
index fa6c25dd301fd1ed400964a07a9b61a2f2b89131..d867152a99bcc877bbed0e72e8df98561f7430d0 100644
--- a/src/hydro/Shadowswift/hydro_debug.h
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
diff --git a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
index 4a7df6560f8409a279e349c3226dcbc60b554b05..9977c36d94d9647da071b867de5a8d101457de6d 100644
--- a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
+++ b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
@@ -50,6 +50,20 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
   hydro_slope_limit_cell_init(p);
 }
 
+/**
+ * @brief Add the gradient estimate for a single quantity due to a particle pair
+ * to the total gradient for that quantity
+ *
+ * This corresponds to one term of equation (21) in Springel (2010).
+ *
+ * @param qL Value of the quantity on the left.
+ * @param qR Value of the quantity on the right.
+ * @param cLR Vector pointing from the midpoint of the particle pair to the
+ * geometrical centroid of the face in between the particles.
+ * @param xLR Vector pointing from the right particle to the left particle.
+ * @param A Surface area of the face in between the particles.
+ * @param grad Current value of the gradient for the quantity (is updated).
+ */
 __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
     float qL, float qR, float *cLR, float *xLR, float rLR, float A,
     float *grad) {
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 40d5da95958fd2109571e083b5c5ecb53de110b9..1c1ec9e1ad1f833d3bf73986026d7d0c96d8ff5f 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -1,8 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *               2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -25,14 +23,9 @@
 #include "voronoi_algorithm.h"
 
 /**
- * @brief Calculate the volume interaction between particle i and particle j
+ * @brief Calculate the Voronoi cell by interacting particle pi and pj
  *
- * The volume is in essence the same as the weighted number of neighbours in a
- * classical SPH density calculation.
- *
- * We also calculate the components of the matrix E, which is used for second
- * order accurate gradient calculations and for the calculation of the interface
- * surface areas.
+ * This method wraps around voronoi_cell_interact().
  *
  * @param r2 Squared distance between particle i and particle j.
  * @param dx Distance vector between the particles (dx = pi->x - pj->x).
@@ -54,15 +47,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 }
 
 /**
- * @brief Calculate the volume interaction between particle i and particle j:
- * non-symmetric version
- *
- * The volume is in essence the same as the weighted number of neighbours in a
- * classical SPH density calculation.
+ * @brief Calculate the Voronoi cell by interacting particle pi with pj
  *
- * We also calculate the components of the matrix E, which is used for second
- * order accurate gradient calculations and for the calculation of the interface
- * surface areas.
+ * This method wraps around voronoi_cell_interact().
  *
  * @param r2 Squared distance between particle i and particle j.
  * @param dx Distance vector between the particles (dx = pi->x - pj->x).
@@ -125,11 +112,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
  * active), both runner_iact_force and runner_iact_nonsym_force call this
  * method, with an appropriate mode.
  *
- * This method calculates the surface area of the interface between particle i
- * and particle j, as well as the interface position and velocity. These are
- * then used to reconstruct and predict the primitive variables, which are then
- * fed to a Riemann solver that calculates a flux. This flux is used to update
- * the conserved variables of particle i or both particles.
+ * This method retrieves the oriented surface area and face midpoint for the
+ * Voronoi face between pi and pj (if it exists). It uses the midpoint position
+ * to reconstruct the primitive quantities (if gradients are used) at the face
+ * and then uses the face quantities to estimate a flux through the face using
+ * a Riemann solver.
  *
  * This method also calculates the maximal velocity used to calculate the time
  * step.
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index b3e01bf1206cbf8ccdbed9183a01b754e560dd19..ffde3095624f80f429fac6de0571152fa5b1f168 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -1,6 +1,6 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
@@ -55,14 +55,35 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_DENSITY, parts, primitives.rho);
 }
 
+/**
+ * @brief Get the internal energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Internal energy of the particle
+ */
 float convert_u(struct engine* e, struct part* p) {
   return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
 }
 
+/**
+ * @brief Get the entropic function of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Entropic function of the particle
+ */
 float convert_A(struct engine* e, struct part* p) {
   return p->primitives.P / pow_gamma(p->primitives.rho);
 }
 
+/**
+ * @brief Get the total energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Total energy of the particle
+ */
 float convert_Etot(struct engine* e, struct part* p) {
   float momentum2;
 
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index f0a360c3e46c55f88f49285fa810d8fc562fa538..46f9ec4509e553ddca179db2e11f28fd7e109b4f 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -1,8 +1,8 @@
 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *                    Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *               2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
index 2e62bce4fc49c97593d73b71ab1c9c5df2353783..a7b3f7511fa7cd247fd9d2399cd200d8d943630e 100644
--- a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
@@ -80,6 +80,17 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
   pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
 }
 
+/**
+ * @brief Apply the cell wide slope limiter to the gradient of a single quantity
+ *
+ * This corresponds to equation (B2) in Hopkins (2015).
+ *
+ * @param grad Gradient to slope limit
+ * @param qval Value of the quantity at the cell generator
+ * @param qmin Minimal value of the quantity among all cell neighbours
+ * @param qmax Maximal value of the quantity among all cell neighbours
+ * @param maxr Maximal distance between the generator and all of its neighbours
+ */
 __attribute__((always_inline)) INLINE static void
 hydro_slope_limit_cell_quantity(float* grad, float qval, float qmin, float qmax,
                                 float maxr) {
diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h
index c58b2474facb3b228a497fb9285abd5f7ebb02dd..222c97b617512b4fcaa1bfd935674bd5b7d8eca5 100644
--- a/src/hydro/Shadowswift/voronoi1d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h
@@ -27,6 +27,16 @@
 #include "inline.h"
 #include "voronoi1d_cell.h"
 
+/**
+ * @brief Initialize a 1D Voronoi cell
+ *
+ * Sets the positions of left and right neighbours to very large values, the
+ * generator position to the given particle position, and all other quantities
+ * to zero.
+ *
+ * @param cell 1D Voronoi cell to initialize.
+ * @param x Position of the generator of the cell.
+ */
 __attribute__((always_inline)) INLINE void voronoi_cell_init(
     struct voronoi_cell *cell, double *x) {
   cell->x = x[0];
@@ -38,6 +48,19 @@ __attribute__((always_inline)) INLINE void voronoi_cell_init(
   cell->centroid = 0.0f;
 }
 
+/**
+ * @brief Interact a 1D Voronoi cell with a particle with given relative
+ * position and ID
+ *
+ * This method checks if the given relative position is closer to the cell
+ * generator than the current left or right neighbour and updates neighbours
+ * accordingly.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param dx Relative position of the interacting generator w.r.t. the cell
+ * generator (in fact: dx = generator - neighbour).
+ * @param id ID of the interacting neighbour.
+ */
 __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     struct voronoi_cell *cell, float *dx, unsigned long long id) {
 
@@ -61,6 +84,22 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
   }
 }
 
+/**
+ * @brief Finalize a 1D Voronoi cell
+ *
+ * Calculates the relative positions of the midpoints of the faces (which in
+ * this case are just the midpoints of the segments connecting the generator
+ * with the two neighbours) w.r.t. the generator, and the cell volume (length)
+ * and centroid (midpoint of the segment connecting the midpoints of the faces).
+ * This function returns the maximal radius at which a particle could still
+ * change the structure of the cell, i.e. twice the largest distance between
+ * the cell generator and one of its faces. If the cell has been interacted with
+ * all neighbours within this radius, we know for sure that the cell is
+ * complete.
+ *
+ * @param cell 1D Voronoi cell.
+ * @return Maximal radius that could still change the structure of the cell.
+ */
 __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
     struct voronoi_cell *cell) {
 
@@ -77,6 +116,28 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
   return max_radius;
 }
 
+/**
+ * @brief Get the oriented surface area and midpoint of the face between a
+ * 1D Voronoi cell and the given neighbour
+ *
+ * This function also checks if the given neighbour is in fact a neighbour of
+ * this cell. Since we perform gradient and flux calculations for all neighbour
+ * pairs within the smoothing length, which assumes the cell to be spherical,
+ * it can happen that this is not the case. It is the responsibility of the
+ * routine that calls this function to check for a zero return value and
+ * deal with it appropriately.
+ *
+ * For this specific case, we simply check if the neighbour is the left or
+ * right neighbour and set the surface area to the unit vector pointing either
+ * to the right or to the left. The midpoint is set to the relative position
+ * vector of the appropriate face.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param ngb ID of a particle that is possibly a neighbour of this cell.
+ * @param A Array to store the oriented surface area in.
+ * @param midpoint Array to store the relative position of the face in.
+ * @return 0 if the given neighbour is not a neighbour, 1 otherwise.
+ */
 __attribute__((always_inline)) INLINE int voronoi_get_face(
     struct voronoi_cell *cell, unsigned long long ngb, float *A,
     float *midpoint) {
@@ -106,6 +167,15 @@ __attribute__((always_inline)) INLINE int voronoi_get_face(
   return 1;
 }
 
+/**
+ * @brief Get the centroid of a 1D Voronoi cell
+ *
+ * We store only the relevant coordinate of the centroid, but need to return
+ * a 3D vector.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param centroid Array to store the centroid in.
+ */
 __attribute__((always_inline)) INLINE void voronoi_get_centroid(
     struct voronoi_cell *cell, float *centroid) {
 
@@ -114,10 +184,4 @@ __attribute__((always_inline)) INLINE void voronoi_get_centroid(
   centroid[2] = 0.0f;
 }
 
-__attribute__((always_inline)) INLINE int voronoi_is_neighbour(
-    struct voronoi_cell *cell, unsigned long long ngb) {
-
-  return (ngb == cell->idL || ngb == cell->idR);
-}
-
 #endif  // SWIFT_VORONOI1D_ALGORITHM_H
diff --git a/src/hydro/Shadowswift/voronoi1d_cell.h b/src/hydro/Shadowswift/voronoi1d_cell.h
index 2dc9bcf4e9a30153d68c1615f10cb5accbaa033e..b2866095939993d4fa09c1ec59ecaa8e47b31798 100644
--- a/src/hydro/Shadowswift/voronoi1d_cell.h
+++ b/src/hydro/Shadowswift/voronoi1d_cell.h
@@ -20,6 +20,7 @@
 #ifndef SWIFT_VORONOI1D_CELL_H
 #define SWIFT_VORONOI1D_CELL_H
 
+/* 1D Voronoi cell */
 struct voronoi_cell {
 
   /* The position of the generator of the cell. */