diff --git a/doc/RTD/source/Planetary/index.rst b/doc/RTD/source/Planetary/index.rst
index 5cc45fa0e51fbfc38e42e7cf111a74ccfe9f835e..2575db1356d8e01beb41dfcc512c2dcc7b637a2b 100644
--- a/doc/RTD/source/Planetary/index.rst
+++ b/doc/RTD/source/Planetary/index.rst
@@ -32,3 +32,4 @@ chosen from the several available equations of state.
    Hydro Scheme <hydro_scheme>
    Equations of State <equations_of_state>
    Initial Conditions <initial_conditions>
+   Removed Particles <removed_particles>
diff --git a/doc/RTD/source/Planetary/removed_particles.rst b/doc/RTD/source/Planetary/removed_particles.rst
new file mode 100755
index 0000000000000000000000000000000000000000..22faadc2810a567dbfb546853711cf1833e9605b
--- /dev/null
+++ b/doc/RTD/source/Planetary/removed_particles.rst
@@ -0,0 +1,25 @@
+.. Planetary Removed Particles
+   Jacob Kegerreis, 27th November 2020
+
+.. _planetary_removed_particles:
+   
+Removed Particles
+=================
+
+If a particle leaves a non-periodic simulation box then it is removed. 
+Currently, the particle's information is simply printed to stdout in a csv 
+format (see ``hydro_remove_part()`` in ``src/hydro/Planetary/hydro.h``).
+This could be extracted for analysis e.g. using grep:
+``grep '## Removed' -A 1 --no-group-separator output.txt > removed.txt``, 
+then read e.g. with python:
+
+.. code-block:: python
+
+    import numpy as np
+    id, x, y, z, vx, vy, vz, mass, u, P, rho, h, mat_id, time = np.loadtxt(
+        "removed.txt", delimiter=",", unpack=True
+    )
+    pos = np.transpose([x, y, z])
+    vel = np.transpose([vx, vy, vz])
+    id = id.astype(int)
+    mat_id = mat_id.astype(int)
\ No newline at end of file
diff --git a/src/cell_drift.c b/src/cell_drift.c
index 2eeb84c7ab066ac7806e6933dd873d37ceb2b520..291316b17b67364942ee0c67a9c1f6b0558d43ba 100644
--- a/src/cell_drift.c
+++ b/src/cell_drift.c
@@ -188,7 +188,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
 #endif
 
             /* One last action before death? */
-            hydro_remove_part(p, xp);
+            hydro_remove_part(p, xp, e->time);
 
             /* Remove the particle entirely */
             cell_remove_part(e, c, p, xp);
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 8e7e995780cde5d51ceb11e61a464041139db553..27c87bbf6cda1c105981dfdee007343f3fff4cee 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -476,9 +476,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 /**
  * @brief Prepares a particle for the density calculation.
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 4f622fbbdb49ce0efa6fdcb22f1ac21344fcd535..ddd20566ad05d9f78696d1152f016ec4e00dad9b 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -905,8 +905,9 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 #endif /* SWIFT_GADGET2_HYDRO_H */
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 4972cf7161c7f97506196640e79d438357b9d55a..07e732432fd37b0a83033680884a3e6b5d7ba58e 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -747,8 +747,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part* p, const struct xpart* xp) {}
+    const struct part* p, const struct xpart* xp, const double time) {}
 
 #endif /* SWIFT_GIZMO_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 8dab9acba8e1ef804b11e0ef9bb56c61eb9c1928..55440943c1a761dcd7f5106f72aa5e202f14e7b8 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -879,8 +879,9 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 #endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Phantom/hydro.h b/src/hydro/Phantom/hydro.h
index 313cac53177c7b61ff582553898d40ae8c898eb4..4d2a64575c35c4c1ede5b1180a7569d4a5e55284 100644
--- a/src/hydro/Phantom/hydro.h
+++ b/src/hydro/Phantom/hydro.h
@@ -458,9 +458,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 /**
  * @brief Prepares a particle for the density calculation.
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 85c02aafc27a0656c05acda78dfcb75757e17b1b..32eccbf8c699e896dca72369dbec94fa447ac496 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -885,13 +885,21 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {
-
-  printf("Removed particle id=%lld \n", p->id);
-  printParticle_single(p, xp);
-  fflush(stdout);
+    const struct part *p, const struct xpart *xp, const double time) {
+
+  /* Print the particle info as csv to facilitate later analysis, e.g. with
+   * grep '## Removed' -A 1 --no-group-separator output.txt > removed.txt
+   */
+  printf(
+      "## Removed particle: "
+      "id, x, y, z, vx, vy, vz, m, u, P, rho, h, mat_id, time \n"
+      "%lld, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, "
+      "%.7g, %.7g, %.7g, %.7g, %.7g, %d, %.7g \n",
+      p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->mass,
+      p->u, p->force.pressure, p->rho, p->h, p->mat_id, time);
 }
 
 #endif /* SWIFT_PLANETARY_HYDRO_H */
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 8ed236716853a8826f0d0c32af81b23632ee1886..d94905b881321110324a32d552637ad9b096cb0f 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -925,8 +925,9 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 #endif /* SWIFT_PRESSURE_ENERGY_HYDRO_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 4aacbe3b3215e293c64cc460b8a5a691bd235436..382c73247eb6afc98cfe0ad98bd7791d8ba7adde 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -924,8 +924,9 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 #endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index f1fefcf72b52c79dbf8c48b5cd51dfc35f77130d..258738b5d2b90a54d2fd4d923ec4b136e295ce96 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -843,8 +843,9 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 #endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h
index ad1533fe586d855e2cbcca61d04d9b5f3f8c2eb6..0fbdf98ed3d365b48db1e75d9464ebd55d84e94a 100644
--- a/src/hydro/SPHENIX/hydro.h
+++ b/src/hydro/SPHENIX/hydro.h
@@ -465,9 +465,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part *p, const struct xpart *xp) {}
+    const struct part *p, const struct xpart *xp, const double time) {}
 
 /**
  * @brief Prepares a particle for the density calculation.
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 907f16c6dc65bf2ef9efcd4ab3501cb2cdab562c..ae1a5cb2cb0e73a5c32ecca285455d8a31e54a27 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -962,8 +962,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_density(
  *
  * @param p The particle.
  * @param xp The extended particle data.
+ * @param time The simulation time.
  */
 __attribute__((always_inline)) INLINE static void hydro_remove_part(
-    const struct part* p, const struct xpart* xp) {}
+    const struct part* p, const struct xpart* xp, const double time) {}
 
 #endif /* SWIFT_SHADOWSWIFT_HYDRO_H */