diff --git a/src/part_type.c b/src/part_type.c
index af97bd34aaace93a9faa953c0c9345d83ca3bc34..ba7efeae6b0c789e7f130f0e0aa45544f66a3ab6 100644
--- a/src/part_type.c
+++ b/src/part_type.c
@@ -21,4 +21,4 @@
 #include "part_type.h"
 
 const char* part_type_names[swift_type_count] = {"Gas",   "DM",   "Dummy",
-                                                 "Dummy", "Star", "BH"};
+                                                 "Dummy", "Stars", "BH"};
diff --git a/src/part_type.h b/src/part_type.h
index fbe2b2aeaea37503635372b0f09f8edde4578721..901f47193fa0e72b362c8dce5199a1d0a20526c9 100644
--- a/src/part_type.h
+++ b/src/part_type.h
@@ -27,7 +27,7 @@
 enum part_type {
   swift_type_gas = 0,
   swift_type_dark_matter = 1,
-  swift_type_star = 4,
+  swift_type_stars = 4,
   swift_type_black_hole = 5,
   swift_type_count
 } __attribute__((packed));
diff --git a/src/runner.c b/src/runner.c
index b1c5c955b7bb9c8f11ae5f69fb06dc4d055a00a7..a9c47e2dd8790ba90121a2785654209fa5e4986a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -268,16 +268,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
 	/* We now have a particle whose smoothing length has converged */
 
-        /* As of here, particle force variables will be set. */
-
-        /* Compute variables required for the force loop */
-        stars_prepare_force(sp, cosmo);
-
-        /* The particle force values are now set.  Do _NOT_
-           try to read any particle density variables! */
-
-        /* Prepare the particle for the force loop over neighbours */
-        stars_reset_acceleration(sp);
+        /* Compute the stellar evolution  */
+        stars_evolve_spart(sp, cosmo);
 
       }
 
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index 6707e0652455600d237cd9e0e7f4ca69627c803e..f8baf1441c5aacae79847b7fdeaf29e82e06c44a 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -27,14 +27,14 @@
  * @param c cell
  * @param timer 1 if the time is to be recorded.
  */
-void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
+void runner_doself_stars_density(struct runner *r, struct cell *c, int timer) {
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
 
   TIMER_TIC;
 
   /* Anything to do here? */
-  if (!cell_is_active_star(c, e)) return;
+  if (!cell_is_active_stars(c, e)) return;
 
   /* Cosmological terms */
   const float a = cosmo->a;
@@ -72,19 +72,17 @@ void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
 
 #ifdef SWIFT_DEBUG_CHECKS
       /* Check that particles have been drifted to the current time */
-      /* if (si->ti_drift != e->ti_current) */
-      /*   error("Particle si not drifted to current time"); */
       if (pj->ti_drift != e->ti_current)
         error("Particle pj not drifted to current time");
 #endif
 
       if (r2 > 0.f && r2 < hig2) {
-	runner_iact_nonsym_star_density(r2, dx, hi, hj, si, pj, a, H);
+	runner_iact_nonsym_stars_density(r2, dx, hi, hj, si, pj, a, H);
       }
     } /* loop over the parts in ci. */
   }   /* loop over the sparts in ci. */
 
-  TIMER_TOC(timer_doself_star_density);
+  TIMER_TOC(timer_doself_stars_density);
  
 }
 
@@ -95,14 +93,14 @@ void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
  * @param c cell
  * @param timer 1 if the time is to be recorded.
  */
-void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
+void runner_dosubpair_stars_density(struct runner *r, struct cell *restrict ci,
 				   struct cell *restrict cj) {
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
 
   /* Anything to do here? */
-  if (!cell_is_active_star(ci, e) && !cell_is_active_star(cj, e)) return;
+  if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
 
   const int scount_i = ci->scount;
   const int count_j = cj->count;
@@ -149,14 +147,12 @@ void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
 
 #ifdef SWIFT_DEBUG_CHECKS
       /* Check that particles have been drifted to the current time */
-      /* if (si->ti_drift != e->ti_current) */
-      /*   error("Particle si not drifted to current time"); */
       if (pj->ti_drift != e->ti_current)
         error("Particle pj not drifted to current time");
 #endif
 
       if (r2 < hig2)
-	runner_iact_nonsym_star_density(r2, dx, hi, hj, si, pj, a, H);
+	runner_iact_nonsym_stars_density(r2, dx, hi, hj, si, pj, a, H);
 
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -164,15 +160,15 @@ void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
 }
 
 
-void runner_dopair_star_density(struct runner *r, struct cell *restrict ci,
+void runner_dopair_stars_density(struct runner *r, struct cell *restrict ci,
 				   struct cell *restrict cj, int timer) {
 
   TIMER_TIC;
   
-  runner_dosubpair_star_density(r, ci, cj);
-  runner_dosubpair_star_density(r, cj, ci);
+  runner_dosubpair_stars_density(r, ci, cj);
+  runner_dosubpair_stars_density(r, cj, ci);
 
-  if (timer) TIMER_TOC(timer_dopair_star_density);
+  if (timer) TIMER_TOC(timer_dopair_stars_density);
 }
 
 
@@ -190,7 +186,7 @@ void runner_dopair_star_density(struct runner *r, struct cell *restrict ci,
  * @param cj The second #cell.
  * @param shift The shift vector to apply to the particles in ci.
  */
-void runner_dopair_subset_star_density(struct runner *r, struct cell *restrict ci,
+void runner_dopair_subset_stars_density(struct runner *r, struct cell *restrict ci,
 				       struct spart *restrict sparts_i, int *restrict ind,
 				       int scount, struct cell *restrict cj,
 				       const double *shift) {
@@ -238,14 +234,12 @@ void runner_dopair_subset_star_density(struct runner *r, struct cell *restrict c
 
 #ifdef SWIFT_DEBUG_CHECKS
       /* Check that particles have been drifted to the current time */
-      /* if (spi->ti_drift != e->ti_current) */
-      /*   error("Particle pi not drifted to current time"); */
       if (pj->ti_drift != e->ti_current)
         error("Particle pj not drifted to current time");
 #endif
       /* Hit or miss? */
       if (r2 < hig2) {
-        runner_iact_nonsym_star_density(r2, dx, hi, pj->h, spi, pj, a, H);
+        runner_iact_nonsym_stars_density(r2, dx, hi, pj->h, spi, pj, a, H);
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -264,7 +258,7 @@ void runner_dopair_subset_star_density(struct runner *r, struct cell *restrict c
  * @param ind The list of indices of particles in @c ci to interact with.
  * @param scount The number of particles in @c ind.
  */
-void runner_doself_subset_star_density(struct runner *r, struct cell *restrict ci,
+void runner_doself_subset_stars_density(struct runner *r, struct cell *restrict ci,
                    struct spart *restrict sparts, int *restrict ind, int scount) {
 
   const struct engine *e = r->e;
@@ -310,20 +304,18 @@ void runner_doself_subset_star_density(struct runner *r, struct cell *restrict c
 
 #ifdef SWIFT_DEBUG_CHECKS
       /* Check that particles have been drifted to the current time */
-      /* if (spi->ti_drift != e->ti_current) */
-      /*   error("Particle pi not drifted to current time"); */
       if (pj->ti_drift != e->ti_current)
         error("Particle pj not drifted to current time");
 #endif
 
       /* Hit or miss? */
       if (r2 > 0.f && r2 < hig2) {
-	runner_iact_nonsym_star_density(r2, dx, hi, hj, spi, pj, a, H);
+	runner_iact_nonsym_stars_density(r2, dx, hi, hj, spi, pj, a, H);
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
 
-  TIMER_TOC(timer_doself_subset_star_density);
+  TIMER_TOC(timer_doself_subset_stars_density);
 }
 
 
@@ -337,11 +329,11 @@ void runner_doself_subset_star_density(struct runner *r, struct cell *restrict c
  * @param ind The list of indices of particles in @c ci to interact with.
  * @param scount The number of particles in @c ind.
  */
-void runner_doself_subset_branch_star_density(struct runner *r, struct cell *restrict ci,
+void runner_doself_subset_branch_stars_density(struct runner *r, struct cell *restrict ci,
                           struct spart *restrict sparts, int *restrict ind,
                           int scount) {
 
-  runner_doself_subset_star_density(r, ci, sparts, ind, scount);
+  runner_doself_subset_stars_density(r, ci, sparts, ind, scount);
 }
 
  /**
@@ -356,7 +348,7 @@ void runner_doself_subset_branch_star_density(struct runner *r, struct cell *res
  * @param scount The number of particles in @c ind.
  * @param cj The second #cell.
  */
- void runner_dopair_subset_branch_star_density(struct runner *r, struct cell *restrict ci,
+ void runner_dopair_subset_branch_stars_density(struct runner *r, struct cell *restrict ci,
                           struct spart *restrict sparts_i, int *restrict ind,
                           int scount, struct cell *restrict cj) {
 
@@ -371,10 +363,10 @@ void runner_doself_subset_branch_star_density(struct runner *r, struct cell *res
       shift[k] = -e->s->dim[k];
   }
 
-  runner_dopair_subset_star_density(r, ci, sparts_i, ind, scount, cj, shift);
+  runner_dopair_subset_stars_density(r, ci, sparts_i, ind, scount, cj, shift);
 }
 
-void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct spart *sparts,
+void runner_dosub_subset_stars_density(struct runner *r, struct cell *ci, struct spart *sparts,
                   int *ind, int scount, struct cell *cj, int sid, int gettimer) {
 
   const struct engine *e = r->e;
@@ -383,8 +375,8 @@ void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct
   TIMER_TIC;
 
   /* Should we even bother? */
-  if (!cell_is_active_star(ci, e) &&
-      (cj == NULL || !cell_is_active_star(cj, e)))
+  if (!cell_is_active_stars(ci, e) &&
+      (cj == NULL || !cell_is_active_stars(cj, e)))
     return;
   if (ci->scount == 0 || (cj != NULL && cj->scount == 0)) return;
 
@@ -409,16 +401,16 @@ void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct
     if (cell_can_recurse_in_self_task(ci)) {
 
       /* Loop over all progeny. */
-      runner_dosub_subset_star_density(r, sub, sparts, ind, scount, NULL, -1, 0);
+      runner_dosub_subset_stars_density(r, sub, sparts, ind, scount, NULL, -1, 0);
       for (int j = 0; j < 8; j++)
         if (ci->progeny[j] != sub && ci->progeny[j] != NULL)
-          runner_dosub_subset_star_density(r, sub, sparts, ind, scount, ci->progeny[j], -1, 0);
+          runner_dosub_subset_stars_density(r, sub, sparts, ind, scount, ci->progeny[j], -1, 0);
 
     }
 
     /* Otherwise, compute self-interaction. */
     else
-      runner_doself_subset_branch_star_density(r, ci, sparts, ind, scount);
+      runner_doself_subset_branch_stars_density(r, ci, sparts, ind, scount);
   } /* self-interaction. */
 
   /* Otherwise, it's a pair interaction. */
@@ -438,496 +430,496 @@ void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct
         /* Regular sub-cell interactions of a single cell. */
         case 0: /* (  1 ,  1 ,  1 ) */
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 1: /* (  1 ,  1 ,  0 ) */
           if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 2: /* (  1 ,  1 , -1 ) */
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           break;
 
         case 3: /* (  1 ,  0 ,  1 ) */
           if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 4: /* (  1 ,  0 ,  0 ) */
           if (ci->progeny[4] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[4] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[4] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 5: /* (  1 ,  0 , -1 ) */
           if (ci->progeny[4] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           break;
 
         case 6: /* (  1 , -1 ,  1 ) */
           if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           break;
 
         case 7: /* (  1 , -1 ,  0 ) */
           if (ci->progeny[4] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           break;
 
         case 8: /* (  1 , -1 , -1 ) */
           if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+            runner_dosub_subset_stars_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
                          -1, 0);
           if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+            runner_dosub_subset_stars_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
                          -1, 0);
           break;
 
         case 9: /* (  0 ,  1 ,  1 ) */
           if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 10: /* (  0 ,  1 ,  0 ) */
           if (ci->progeny[2] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[2] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[2] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[2] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
 
         case 11: /* (  0 ,  1 , -1 ) */
           if (ci->progeny[2] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[2] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[2] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           if (ci->progeny[6] == sub && cj->progeny[5] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
+            runner_dosub_subset_stars_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
                          -1, 0);
           if (ci->progeny[6] != NULL && cj->progeny[5] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
+            runner_dosub_subset_stars_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
                          -1, 0);
           break;
 
         case 12: /* (  0 ,  0 ,  1 ) */
           if (ci->progeny[1] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[1] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[1],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[1],
                          -1, 0);
           if (ci->progeny[1] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[1] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[1],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[1],
                          -1, 0);
           if (ci->progeny[1] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[1] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[1],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[1],
                          -1, 0);
           if (ci->progeny[1] == sub && cj->progeny[6] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[6],
+            runner_dosub_subset_stars_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[6],
                          -1, 0);
           if (ci->progeny[1] != NULL && cj->progeny[6] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[1],
+            runner_dosub_subset_stars_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[1],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[3] == sub && cj->progeny[6] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[6],
+            runner_dosub_subset_stars_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[6],
                          -1, 0);
           if (ci->progeny[3] != NULL && cj->progeny[6] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[3],
+            runner_dosub_subset_stars_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[3],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[5] == sub && cj->progeny[6] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[6],
+            runner_dosub_subset_stars_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[6],
                          -1, 0);
           if (ci->progeny[5] != NULL && cj->progeny[6] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[5],
+            runner_dosub_subset_stars_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[5],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           if (ci->progeny[7] == sub && cj->progeny[6] != NULL)
-            runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[6],
+            runner_dosub_subset_stars_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[6],
                          -1, 0);
           if (ci->progeny[7] != NULL && cj->progeny[6] == sub)
-            runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[7],
+            runner_dosub_subset_stars_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[7],
                          -1, 0);
           break;
       }
@@ -935,12 +927,12 @@ void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct
     }
 
     /* Otherwise, compute the pair directly. */
-    else if (cell_is_active_star(ci, e) || cell_is_active_star(cj, e)) {
+    else if (cell_is_active_stars(ci, e) || cell_is_active_stars(cj, e)) {
 
       /* Do any of the cells need to be drifted first? */
       if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
 
-      runner_dopair_subset_branch_star_density(r, ci, sparts, ind, scount, cj);
+      runner_dopair_subset_branch_stars_density(r, ci, sparts, ind, scount, cj);
     }
 
   } /* otherwise, pair interaction. */
@@ -957,18 +949,18 @@ void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct
  * @param c #cell c
  *
  */
-void runner_doself_branch_star_density(struct runner *r, struct cell *c) {
+void runner_doself_branch_stars_density(struct runner *r, struct cell *c) {
 
   const struct engine *restrict e = r->e;
 
   /* Anything to do here? */
-  if (!cell_is_active_star(c, e)) return;
+  if (!cell_is_active_stars(c, e)) return;
 
   /* Did we mess up the recursion? */
   if (c->h_max_old * kernel_gamma > c->dmin)
     error("Cell smaller than smoothing length");
 
-  runner_doself_star_density(r, c, 1);
+  runner_doself_stars_density(r, c, 1);
 }
 
 
@@ -981,16 +973,16 @@ void runner_doself_branch_star_density(struct runner *r, struct cell *c) {
  * @param cj #cell cj
  *
  */
-void runner_dopair_branch_star_density(struct runner *r, struct cell *ci, struct cell *cj) {
+void runner_dopair_branch_stars_density(struct runner *r, struct cell *ci, struct cell *cj) {
 
   const struct engine *restrict e = r->e;
 
   /* Anything to do here? */
-  if (!cell_is_active_star(ci, e) && !cell_is_active_star(cj, e)) return;
+  if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
 
-  /* /\* Check that cells are drifted. *\/ */
-  /* if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e)) */
-  /*   error("Interacting undrifted cells."); */
+  /* Check that cells are drifted. */
+  if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
+    error("Interacting undrifted cells.");
 
   /* Get the sort ID. */
   double shift[3] = {0.0, 0.0, 0.0};
@@ -1043,7 +1035,7 @@ void runner_dopair_branch_star_density(struct runner *r, struct cell *ci, struct
   }
 #endif /* SWIFT_DEBUG_CHECKS */
 
-  runner_dopair_star_density(r, ci, cj, 1);
+  runner_dopair_stars_density(r, ci, cj, 1);
 }
 
 /**
@@ -1058,7 +1050,7 @@ void runner_dopair_branch_star_density(struct runner *r, struct cell *ci, struct
  * @todo Hard-code the sid on the recursive calls to avoid the
  * redundant computations to find the sid on-the-fly.
  */
-void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct cell *cj, int sid,
+void runner_dosub_pair_stars_density(struct runner *r, struct cell *ci, struct cell *cj, int sid,
                  int gettimer) {
 
   struct space *s = r->e->s;
@@ -1067,7 +1059,7 @@ void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct ce
   TIMER_TIC;
 
   /* Should we even bother? */
-  if (!cell_is_active_star(ci, e) && !cell_is_active_star(cj, e)) return;
+  if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
   if (ci->scount == 0 || cj->scount == 0) return;
 
   /* Get the type of pair if not specified explicitly. */
@@ -1083,200 +1075,200 @@ void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct ce
       /* Regular sub-cell interactions of a single cell. */
       case 0: /* (  1 ,  1 ,  1 ) */
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         break;
 
       case 1: /* (  1 ,  1 ,  0 ) */
         if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
         break;
 
       case 2: /* (  1 ,  1 , -1 ) */
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         break;
 
       case 3: /* (  1 ,  0 ,  1 ) */
         if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
         break;
 
       case 4: /* (  1 ,  0 ,  0 ) */
         if (ci->progeny[4] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[0], -1, 0);
         if (ci->progeny[4] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[1], -1, 0);
         if (ci->progeny[4] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[2], -1, 0);
         if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[1], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[3], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[2], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[3], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[3], -1, 0);
         break;
 
       case 5: /* (  1 ,  0 , -1 ) */
         if (ci->progeny[4] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[1], -1, 0);
         if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[3], -1, 0);
         break;
 
       case 6: /* (  1 , -1 ,  1 ) */
         if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
         break;
 
       case 7: /* (  1 , -1 ,  0 ) */
         if (ci->progeny[4] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[2], -1, 0);
         if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[3], -1, 0);
         break;
 
       case 8: /* (  1 , -1 , -1 ) */
         if (ci->progeny[4] != NULL && cj->progeny[3] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[4], cj->progeny[3], -1, 0);
         break;
 
       case 9: /* (  0 ,  1 ,  1 ) */
         if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
         break;
 
       case 10: /* (  0 ,  1 ,  0 ) */
         if (ci->progeny[2] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[0], -1, 0);
         if (ci->progeny[2] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[1], -1, 0);
         if (ci->progeny[2] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[4], -1, 0);
         if (ci->progeny[2] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[5], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[1], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[5], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[0], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[4], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[5], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[1], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[5], -1, 0);
         break;
 
       case 11: /* (  0 ,  1 , -1 ) */
         if (ci->progeny[2] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[1], -1, 0);
         if (ci->progeny[2] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[2], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[2], cj->progeny[5], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[1] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[1], -1, 0);
         if (ci->progeny[6] != NULL && cj->progeny[5] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[6], cj->progeny[5], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[6], cj->progeny[5], -1, 0);
         break;
 
       case 12: /* (  0 ,  0 ,  1 ) */
         if (ci->progeny[1] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[1], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[1], cj->progeny[0], -1, 0);
         if (ci->progeny[1] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[1], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[1], cj->progeny[2], -1, 0);
         if (ci->progeny[1] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[1], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[1], cj->progeny[4], -1, 0);
         if (ci->progeny[1] != NULL && cj->progeny[6] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[1], cj->progeny[6], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[1], cj->progeny[6], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[0], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[2], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[4], -1, 0);
         if (ci->progeny[3] != NULL && cj->progeny[6] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[3], cj->progeny[6], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[3], cj->progeny[6], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[0], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[2], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[4], -1, 0);
         if (ci->progeny[5] != NULL && cj->progeny[6] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[5], cj->progeny[6], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[5], cj->progeny[6], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[0] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[0], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[2] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[2], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[4] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[4], -1, 0);
         if (ci->progeny[7] != NULL && cj->progeny[6] != NULL)
-          runner_dosub_pair_star_density(r, ci->progeny[7], cj->progeny[6], -1, 0);
+          runner_dosub_pair_stars_density(r, ci->progeny[7], cj->progeny[6], -1, 0);
         break;
     }
 
   }
 
   /* Otherwise, compute the pair directly. */
-  else if (cell_is_active_star(ci, e) || cell_is_active_star(cj, e)) {
+  else if (cell_is_active_stars(ci, e) || cell_is_active_stars(cj, e)) {
 
     /* Make sure both cells are drifted to the current timestep. */
     if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
@@ -1291,7 +1283,7 @@ void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct ce
       error("Interacting unsorted cell.");
 
     /* Compute the interactions. */
-    runner_dopair_branch_star_density(r, ci, cj);
+    runner_dopair_branch_stars_density(r, ci, cj);
   }
 
   if (gettimer) TIMER_TOC(TIMER_DOSUB_PAIR);
@@ -1305,12 +1297,12 @@ void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct ce
  * @param ci The first #cell.
  * @param gettimer Do we have a timer ?
  */
-void runner_dosub_self_star_density(struct runner *r, struct cell *ci, int gettimer) {
+void runner_dosub_self_stars_density(struct runner *r, struct cell *ci, int gettimer) {
 
   TIMER_TIC;
 
   /* Should we even bother? */
-  if (ci->scount == 0 || !cell_is_active_star(ci, r->e)) return;
+  if (ci->scount == 0 || !cell_is_active_stars(ci, r->e)) return;
 
   /* Recurse? */
   if (cell_can_recurse_in_self_task(ci)) {
@@ -1318,18 +1310,18 @@ void runner_dosub_self_star_density(struct runner *r, struct cell *ci, int getti
     /* Loop over all progeny. */
     for (int k = 0; k < 8; k++)
       if (ci->progeny[k] != NULL) {
-        runner_dosub_self_star_density(r, ci->progeny[k], 0);
+        runner_dosub_self_stars_density(r, ci->progeny[k], 0);
         for (int j = k + 1; j < 8; j++)
           if (ci->progeny[j] != NULL)
-            runner_dosub_pair_star_density(r, ci->progeny[k], ci->progeny[j], -1, 0);
+            runner_dosub_pair_stars_density(r, ci->progeny[k], ci->progeny[j], -1, 0);
       }
   }
 
   /* Otherwise, compute self-interaction. */
  else {
 
-    runner_doself_branch_star_density(r, ci);
+    runner_doself_branch_stars_density(r, ci);
   }
 
-  if (gettimer) TIMER_TOC(timer_dosub_self_star_density);
+  if (gettimer) TIMER_TOC(timer_dosub_self_stars_density);
 }
diff --git a/src/scheduler.c b/src/scheduler.c
index 349a746003a582a05f87c4cdd641e70c548429d6..391ebe47876490d852bbf70a76bae95930318b03 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -241,7 +241,8 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
   int density_cluster[4] = {0};
   int gradient_cluster[4] = {0};
   int force_cluster[4] = {0};
-  int star_density_cluster[4] = {0};
+  int gravity_cluster[5] = {0};
+  int stars_density_cluster[4] = {0};
 
   /* Check whether we need to construct a group of tasks */
   for (int type = 0; type < task_type_count; ++type) {
@@ -262,8 +263,8 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
             force_cluster[k] = 1;
           if (type == task_type_self + k && subtype == task_subtype_grav)
             gravity_cluster[k] = 1;
-	  if (type == task_type_self + k && subtype == task_subtype_star_density)
-	    star_density_cluster[k] = 1;
+	  if (type == task_type_self + k && subtype == task_subtype_stars_density)
+	    stars_density_cluster[k] = 1;
         }
         if (type == task_type_grav_mesh) gravity_cluster[2] = 1;
         if (type == task_type_grav_long_range) gravity_cluster[3] = 1;
@@ -318,9 +319,9 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
   fprintf(f, "\t subgraph cluster4{\n");
   fprintf(f, "\t\t label=\"\";\n");
   for (int k = 0; k < 4; ++k)
-    if (star_density_cluster[k])
+    if (stars_density_cluster[k])
       fprintf(f, "\t\t \"%s %s\";\n", taskID_names[task_type_self + k],
-              subtaskID_names[task_subtype_star_density]);
+              subtaskID_names[task_subtype_stars_density]);
   fprintf(f, "\t};\n");
 
   /* Be clean */
@@ -988,7 +989,7 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
       scheduler_splittask_gravity(t, s);
     } else if (t->type == task_type_grav_mesh) {
       /* For future use */
-    } else if (t->subtype == task_subtype_star_density) {
+    } else if (t->subtype == task_subtype_stars_density) {
       /* For future use */
     } else {
 #ifdef SWIFT_DEBUG_CHECKS
@@ -1338,7 +1339,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
       case task_type_extra_ghost:
         if (t->ci == t->ci->super_hydro) cost = wscale * count_i;
         break;
-      case task_type_star_ghost:
+      case task_type_stars_ghost:
         if (t->ci == t->ci->super_hydro) cost = wscale * scount_i;
 	break;
       case task_type_drift_part:
@@ -1542,7 +1543,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
         break;
       case task_type_sort:
       case task_type_ghost:
-      case task_type_star_ghost:
       case task_type_drift_part:
         qid = t->ci->super_hydro->owner;
         break;
@@ -1551,6 +1551,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
         break;
       case task_type_kick1:
       case task_type_kick2:
+      case task_type_stars_ghost:
       case task_type_timestep:
         qid = t->ci->super->owner;
         break;
diff --git a/src/space.c b/src/space.c
index b4d9d3ab35ce64c10e6d0af642caeff38c476b55..e3701dc9249d8b7170b3f0f99bd4e57b5138c1b0 100644
--- a/src/space.c
+++ b/src/space.c
@@ -182,9 +182,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->ghost_in = NULL;
     c->ghost_out = NULL;
     c->ghost = NULL;
-    c->star_ghost_in = NULL;
-    c->star_ghost_out = NULL;
-    c->star_ghost = NULL;
+    c->stars_ghost_in = NULL;
+    c->stars_ghost_out = NULL;
+    c->stars_ghost = NULL;
     c->kick1 = NULL;
     c->kick2 = NULL;
     c->timestep = NULL;
@@ -658,13 +658,13 @@ void space_rebuild(struct space *s, int verbose) {
       /* Swap the link with part/spart */
       if (s->gparts[k].type == swift_type_gas) {
         s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
-      } else if (s->gparts[k].type == swift_type_star) {
+      } else if (s->gparts[k].type == swift_type_stars) {
         s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
       }
       if (s->gparts[nr_gparts].type == swift_type_gas) {
         s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
             &s->gparts[nr_gparts];
-      } else if (s->gparts[nr_gparts].type == swift_type_star) {
+      } else if (s->gparts[nr_gparts].type == swift_type_stars) {
         s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
             &s->gparts[nr_gparts];
       }
@@ -1565,7 +1565,7 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
         memswap(&ind[j], &target_cid, sizeof(int));
         if (gparts[j].type == swift_type_gas) {
           parts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
-        } else if (gparts[j].type == swift_type_star) {
+        } else if (gparts[j].type == swift_type_stars) {
           sparts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
         }
       }
@@ -1573,7 +1573,7 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
       ind[k] = target_cid;
       if (gparts[k].type == swift_type_gas) {
         parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
-      } else if (gparts[k].type == swift_type_star) {
+      } else if (gparts[k].type == swift_type_stars) {
         sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
       }
     }
@@ -2389,7 +2389,7 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
       xp->v_full[2] = gp->v_full[2];
     }
 
-    else if (gp->type == swift_type_star) {
+    else if (gp->type == swift_type_stars) {
 
       /* Get it's stellar friend */
       struct spart *sp = &s->sparts[-gp->id_or_neg_offset];
@@ -2588,7 +2588,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
   /* Initialise the rest */
   for (int k = 0; k < count; k++) {
 
-    star_first_init_spart(&sp[k]);
+    stars_first_init_spart(&sp[k]);
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (sp[k].gpart && sp[k].gpart->id_or_neg_offset != -(k + delta))
@@ -2604,7 +2604,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
 /**
  * @brief Initialises all the s-particles by setting them into a valid state
  *
- * Calls star_first_init_spart() on all the particles
+ * Calls stars_first_init_spart() on all the particles
  */
 void space_first_init_sparts(struct space *s, int verbose) {
   const ticks tic = getticks();
@@ -2673,7 +2673,7 @@ void space_init_sparts_mapper(void *restrict map_data, int scount,
                              void *restrict extra_data) {
 
   struct spart *restrict sparts = (struct spart *)map_data;
-  for (int k = 0; k < scount; k++) star_init_spart(&sparts[k]);
+  for (int k = 0; k < scount; k++) stars_init_spart(&sparts[k]);
 }
 
 /**
@@ -2736,10 +2736,10 @@ void space_convert_quantities(struct space *s, int verbose) {
  * @param dim Spatial dimensions of the domain.
  * @param parts Array of Gas particles.
  * @param gparts Array of Gravity particles.
- * @param sparts Array of star particles.
+ * @param sparts Array of stars particles.
  * @param Npart The number of Gas particles in the space.
  * @param Ngpart The number of Gravity particles in the space.
- * @param Nspart The number of star particles in the space.
+ * @param Nspart The number of stars particles in the space.
  * @param periodic flag whether the domain is periodic or not.
  * @param replicate How many replications along each direction do we want?
  * @param generate_gas_in_ics Are we generating gas particles from the gparts?
@@ -3415,7 +3415,7 @@ void space_struct_restore(struct space *s, FILE *stream) {
                         NULL, "sparts");
   }
 
-  /* Need to reconnect the gravity parts to their hydro and star particles. */
+  /* Need to reconnect the gravity parts to their hydro and stars particles. */
   /* Re-link the parts. */
   if (s->nr_parts > 0 && s->nr_gparts > 0)
     part_relink_parts_to_gparts(s->gparts, s->nr_gparts, s->parts);
diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index 07df0ca5458ad5f243a88f2e2625d9fd9a58b7ef..65baba7b41c6e112ad6ec26c23c5505f685e10fe 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -132,25 +132,18 @@ __attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
 }
 
 /**
- * @brief Prepare a particle for the force calculation.
+ * @brief Evolve the stellar properties of a #spart.
  *
+ * This function allows for example to compute the SN rate before sending
+ * this information to a different MPI rank.
  *
  * @param sp The particle to act upon
  * @param cosmo The current cosmological model.
+ * @param stars_properties The #stars_props
  */
-__attribute__((always_inline)) INLINE static void stars_prepare_force(
-    struct spart *restrict sp, const struct cosmology *cosmo) {}
+__attribute__((always_inline)) INLINE static void stars_evolve_spart(
+    struct spart *restrict sp, const struct stars_props *stars_properties,
+    const struct cosmology *cosmo) {}
 
 
-/**
- * @brief Reset acceleration fields of a particle
- *
- * Resets all star acceleration and time derivative fields in preparation
- * for the sums taking place in the variaous force tasks
- *
- * @param sp The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_reset_acceleration(
-    struct spart *restrict p) {}
-
 #endif /* SWIFT_DEFAULT_STARS_H */
diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
index fae3ce514803872722b38eb4b852dfdaef48560b..8ae42389e1570b62152361d9d8bdc4afa36e451c 100644
--- a/src/stars/Default/stars_part.h
+++ b/src/stars/Default/stars_part.h
@@ -50,11 +50,14 @@ struct spart {
   /*! Particle time bin */
   timebin_t time_bin;
 
-  /* Number of neighbours. */
-  float wcount;
+  struct {
+    /* Number of neighbours. */
+    float wcount;
 
-  /* Number of neighbours spatial derivative. */
-  float wcount_dh;
+    /* Number of neighbours spatial derivative. */
+    float wcount_dh;
+
+  } density;
 
 #ifdef SWIFT_DEBUG_CHECKS