diff --git a/src/Makefile.am b/src/Makefile.am
index 0f61fb108d8d8acbd420c994266f4803a3f69d3e..59b77dfb4b3bb5ece47e6b970d5d35c7708cd793 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -126,8 +126,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
                  riemann/riemann_checks.h \
 	 	 stars.h stars_io.h \
-		 stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
-		 stars/Default/star_debug.h stars/Default/star_part.h  \
+		 stars/Default/stars.h stars/Default/stars_iact.h stars/Default/stars_io.h \
+		 stars/Default/stars_debug.h stars/Default/stars_part.h  \
 	         potential/none/potential.h potential/point_mass/potential.h \
                  potential/isothermal/potential.h potential/disc_patch/potential.h \
                  potential/sine_wave/potential.h \
diff --git a/src/cell.c b/src/cell.c
index ef4ad7f5cabb5c5885835e17101dd5c7ffabe7f9..c4af3916443a42c82d8dc5c8e8671d1e49a5c0a9 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1115,7 +1115,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
           if (gparts[j].type == swift_type_gas) {
             parts[-gparts[j].id_or_neg_offset - parts_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 - sparts_offset].gpart =
                 &gparts[j];
           }
@@ -1125,7 +1125,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
         gbuff[k] = temp_buff;
         if (gparts[k].type == swift_type_gas) {
           parts[-gparts[k].id_or_neg_offset - parts_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 - sparts_offset].gpart =
               &gparts[k];
         }
diff --git a/src/cell.h b/src/cell.h
index 49fd3d7008e297c4038bf3d2ec4b4b5659e5692e..f29a9d3d34cfde3c1493e4778cc3b2faa908e8a5 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -591,6 +591,30 @@ cell_can_recurse_in_self_hydro_task(const struct cell *c) {
   return c->split && (kernel_gamma * c->h_max_old < 0.5f * c->dmin);
 }
 
+/**
+ * @brief Can a sub-pair star task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_pair_stars_task(const struct cell *c) {
+
+  return 0;
+}
+
+/**
+ * @brief Can a sub-self stars task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_self_stars_task(const struct cell *c) {
+
+  return 0;
+}
+
 /**
  * @brief Can a pair hydro task associated with a cell be split into smaller
  * sub-tasks.
diff --git a/src/common_io.c b/src/common_io.c
index 68311107575a89ce8a2990a8e0f7a8eeb5d2d644..521f08a359e803d48f56028db1308428976ac36d 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -747,7 +747,7 @@ void io_duplicate_hydro_sparts_mapper(void* restrict data, int Nstars,
     gparts[i + Ndm].mass = sparts[i].mass;
 
     /* Set gpart type */
-    gparts[i + Ndm].type = swift_type_star;
+    gparts[i + Ndm].type = swift_type_stars;
 
     /* Link the particles */
     gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
@@ -768,7 +768,7 @@ void io_duplicate_hydro_sparts_mapper(void* restrict data, int Nstars,
  * @param Nstars The number of stars particles read in.
  * @param Ndm The number of DM and gas particles read in.
  */
-void io_duplicate_star_gparts(struct threadpool* tp, struct spart* const sparts,
+void io_duplicate_stars_gparts(struct threadpool* tp, struct spart* const sparts,
                               struct gpart* const gparts, size_t Nstars,
                               size_t Ndm) {
 
@@ -853,8 +853,8 @@ void io_check_output_fields(const struct swift_params* params,
         darkmatter_write_particles(&gp, list, &num_fields);
         break;
 
-      case swift_type_star:
-        star_write_particles(&sp, list, &num_fields);
+      case swift_type_stars:
+        stars_write_particles(&sp, list, &num_fields);
         break;
 
       default:
@@ -939,8 +939,8 @@ void io_write_output_field_parameter(const char* filename) {
         darkmatter_write_particles(NULL, list, &num_fields);
         break;
 
-      case swift_type_star:
-        star_write_particles(NULL, list, &num_fields);
+      case swift_type_stars:
+        stars_write_particles(NULL, list, &num_fields);
         break;
 
       default:
diff --git a/src/common_io.h b/src/common_io.h
index 152b40a8d7c931b3398f4f04d3a61e9cf7f1836c..ce7cba61ffaf4653cbb69f1cb1e6dc5be2cfd61c 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -104,7 +104,7 @@ void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
 void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
                                struct gpart* const gparts, size_t Ngas,
                                size_t Ndm);
-void io_duplicate_star_gparts(struct threadpool* tp, struct spart* const sparts,
+void io_duplicate_stars_gparts(struct threadpool* tp, struct spart* const sparts,
                               struct gpart* const gparts, size_t Nstars,
                               size_t Ndm);
 
diff --git a/src/engine.c b/src/engine.c
index 292facd4d1aef1432f37cf23236a3231343d3dcd..990ba8e4decdea58c066d2872d45b32c72cb93ab 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -818,7 +818,7 @@ static void engine_redistribute_relink_mapper(void *map_data, int num_elements,
       }
 
       /* Does this gpart have a star partner ? */
-      else if (s->gparts[k].type == swift_type_star) {
+      else if (s->gparts[k].type == swift_type_stars) {
 
         const ptrdiff_t partner_index =
             offset_sparts - s->gparts[k].id_or_neg_offset;
@@ -2015,7 +2015,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
     for (size_t k = 0; k < offset_gparts; k++) {
       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];
       }
     }
@@ -2111,7 +2111,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
               &s->parts[offset_parts + count_parts - gp->id_or_neg_offset];
           gp->id_or_neg_offset = s->parts - p;
           p->gpart = gp;
-        } else if (gp->type == swift_type_star) {
+        } else if (gp->type == swift_type_stars) {
           struct spart *sp =
               &s->sparts[offset_sparts + count_sparts - gp->id_or_neg_offset];
           gp->id_or_neg_offset = s->sparts - sp;
@@ -3734,7 +3734,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Activate the star density */
       else if (t->type == task_type_self &&
                t->subtype == task_subtype_stars_density) {
-        if (cell_is_active_star(ci, e)) {
+        if (cell_is_active_stars(ci, e)) {
 	  scheduler_activate(s, t);
           cell_activate_drift_part(ci, s);
 	}
@@ -3743,7 +3743,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Store current values of dx_max and h_max. */
       else if (t->type == task_type_sub_self &&
                t->subtype == task_subtype_stars_density) {
-        if (cell_is_active_star(ci, e)) {
+        if (cell_is_active_stars(ci, e)) {
           scheduler_activate(s, t);
           cell_activate_subcell_hydro_tasks(ci, NULL, s);
         }
@@ -4120,10 +4120,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         scheduler_activate(s, t);
     }
 
+    }
+
     /* Star ghost tasks ? */
     else if (t->type == task_type_stars_ghost ||
              t->type == task_type_stars_ghost_in || t->type == task_type_stars_ghost_out) {
-      if (cell_is_active_star(t->ci, e)) scheduler_activate(s, t);
+      if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
     }
 
     /* Time-step? */
diff --git a/src/kick.h b/src/kick.h
index 50ecaea498bdd401cc0ac27525ed27986a344c59..8c1c55519c311888ea1271f2285841731d7407ba 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -150,7 +150,7 @@ __attribute__((always_inline)) INLINE static void kick_spart(
   sp->gpart->v_full[2] = sp->v[2];
 
   /* Kick extra variables */
-  star_kick_extra(sp, dt_kick_grav);
+  stars_kick_extra(sp, dt_kick_grav);
 }
 
 #endif /* SWIFT_KICK_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 5974070fac645472bc3f68da1a9f505d60813e11..7578be12678b6c83be4c69115697f67c90b11149 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -779,12 +779,12 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
     bzero(*parts, *Ngas * sizeof(struct part));
   }
 
-  /* Allocate memory to store star particles */
+  /* Allocate memory to store stars particles */
   if (with_stars) {
-    *Nstars = N[swift_type_star];
+    *Nstars = N[swift_type_stars];
     if (posix_memalign((void**)sparts, spart_align,
                        *Nstars * sizeof(struct spart)) != 0)
-      error("Error while allocating memory for star particles");
+      error("Error while allocating memory for stars particles");
     bzero(*sparts, *Nstars * sizeof(struct spart));
   }
 
@@ -793,7 +793,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
     Ndm = N[1];
     *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
                N[swift_type_dark_matter] +
-               (with_stars ? N[swift_type_star] : 0);
+               (with_stars ? N[swift_type_stars] : 0);
     if (posix_memalign((void**)gparts, gpart_align,
                        *Ngparts * sizeof(struct gpart)) != 0)
       error("Error while allocating memory for gravity particles");
@@ -842,10 +842,10 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
         }
         break;
 
-      case swift_type_star:
+      case swift_type_stars:
         if (with_stars) {
           Nparticles = *Nstars;
-          star_read_particles(*sparts, list, &num_fields);
+          stars_read_particles(*sparts, list, &num_fields);
         }
         break;
 
@@ -878,9 +878,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
     /* Duplicate the hydro particles into gparts */
     if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
 
-    /* Duplicate the star particles into gparts */
+    /* Duplicate the stars particles into gparts */
     if (with_stars)
-      io_duplicate_star_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+      io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
 
     threadpool_clean(&tp);
   }
@@ -1105,8 +1105,8 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
         darkmatter_write_particles(gparts, list, &num_fields);
         break;
 
-      case swift_type_star:
-        star_write_particles(sparts, list, &num_fields);
+      case swift_type_stars:
+        stars_write_particles(sparts, list, &num_fields);
         break;
 
       default:
@@ -1363,9 +1363,9 @@ void write_output_parallel(struct engine* e, const char* baseName,
         darkmatter_write_particles(dmparts, list, &num_fields);
         break;
 
-      case swift_type_star:
+      case swift_type_stars:
         Nparticles = Nstars;
-        star_write_particles(sparts, list, &num_fields);
+        stars_write_particles(sparts, list, &num_fields);
         break;
 
       default:
diff --git a/src/part.c b/src/part.c
index 050e10e9cdd0ab56adcd34ba3e6f2d35c274f14a..db221dbd4bf9ff2b829b02fbae673aa1c65f339e 100644
--- a/src/part.c
+++ b/src/part.c
@@ -88,7 +88,7 @@ void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
 void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
                                   struct spart *sparts) {
   for (size_t k = 0; k < N; k++) {
-    if (gparts[k].type == swift_type_star) {
+    if (gparts[k].type == swift_type_stars) {
       sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
     }
   }
@@ -108,7 +108,7 @@ void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
   for (size_t k = 0; k < N; k++) {
     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];
     }
   }
@@ -171,11 +171,11 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
         error("Linked particles are not at the same time !");
     }
 
-    else if (gparts[k].type == swift_type_star) {
+    else if (gparts[k].type == swift_type_stars) {
 
       /* Check that it is linked */
       if (gparts[k].id_or_neg_offset > 0)
-        error("Star gpart not linked to anything !");
+        error("Stars gpart not linked to anything !");
 
       /* Find its link */
       const struct spart *spart = &sparts[-gparts[k].id_or_neg_offset];
diff --git a/src/runner.c b/src/runner.c
index a9c47e2dd8790ba90121a2785654209fa5e4986a..cf70cef08ee6ab597f20efab91cc21d0d343f8db 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -202,7 +202,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
         float h_new;
         int has_no_neighbours = 0;
 
-        if (sp->wcount == 0.f) { /* No neighbours case */
+        if (sp->density.wcount == 0.f) { /* No neighbours case */
 
           /* Flag that there were no neighbours */
           has_no_neighbours = 1;
@@ -215,12 +215,12 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
           stars_end_density(sp, cosmo);
 
           /* Compute one step of the Newton-Raphson scheme */
-          const float n_sum = sp->wcount * h_old_dim;
+          const float n_sum = sp->density.wcount * h_old_dim;
           const float n_target = stars_eta_dim;
           const float f = n_sum - n_target;
           const float f_prime =
-              sp->wcount_dh * h_old_dim +
-              hydro_dimension * sp->wcount * h_old_dim_minus_one;
+              sp->density.wcount_dh * h_old_dim +
+              hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
 
           /* Avoid floating point exception from f_prime = 0 */
           h_new = h_old - f / (f_prime + FLT_MIN);
@@ -269,7 +269,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 	/* We now have a particle whose smoothing length has converged */
 
         /* Compute the stellar evolution  */
-        stars_evolve_spart(sp, cosmo);
+        stars_evolve_spart(sp, stars_properties, cosmo);
 
       }
 
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index f8baf1441c5aacae79847b7fdeaf29e82e06c44a..0ad42e1137d3507eef511464c053e89acb4083ff 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -398,7 +398,7 @@ void runner_dosub_subset_stars_density(struct runner *r, struct cell *ci, struct
   if (cj == NULL) {
 
     /* Recurse? */
-    if (cell_can_recurse_in_self_task(ci)) {
+    if (cell_can_recurse_in_self_stars_task(ci)) {
 
       /* Loop over all progeny. */
       runner_dosub_subset_stars_density(r, sub, sparts, ind, scount, NULL, -1, 0);
@@ -417,8 +417,8 @@ void runner_dosub_subset_stars_density(struct runner *r, struct cell *ci, struct
   else {
 
     /* Recurse? */
-    if (cell_can_recurse_in_pair_task(ci) &&
-        cell_can_recurse_in_pair_task(cj)) {
+    if (cell_can_recurse_in_pair_stars_task(ci) &&
+        cell_can_recurse_in_pair_stars_task(cj)) {
 
       /* Get the type of pair if not specified explicitly. */
       double shift[3] = {0.0, 0.0, 0.0};
@@ -1067,7 +1067,7 @@ void runner_dosub_pair_stars_density(struct runner *r, struct cell *ci, struct c
   sid = space_getsid(s, &ci, &cj, shift);
 
   /* Recurse? */
-  if (cell_can_recurse_in_pair_task(ci) && cell_can_recurse_in_pair_task(cj)) {
+  if (cell_can_recurse_in_pair_stars_task(ci) && cell_can_recurse_in_pair_stars_task(cj)) {
 
     /* Different types of flags. */
     switch (sid) {
@@ -1305,7 +1305,7 @@ void runner_dosub_self_stars_density(struct runner *r, struct cell *ci, int gett
   if (ci->scount == 0 || !cell_is_active_stars(ci, r->e)) return;
 
   /* Recurse? */
-  if (cell_can_recurse_in_self_task(ci)) {
+  if (cell_can_recurse_in_self_stars_task(ci)) {
 
     /* Loop over all progeny. */
     for (int k = 0; k < 8; k++)
diff --git a/src/serial_io.c b/src/serial_io.c
index 8c9298951311536adb52b96638b2472e089047c7..87d02165f7ef5880b9ec819cabd43df0cda3455f 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -580,12 +580,12 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     bzero(*parts, *Ngas * sizeof(struct part));
   }
 
-  /* Allocate memory to store star particles */
+  /* Allocate memory to store stars particles */
   if (with_stars) {
-    *Nstars = N[swift_type_star];
+    *Nstars = N[swift_type_stars];
     if (posix_memalign((void*)sparts, spart_align,
                        *Nstars * sizeof(struct spart)) != 0)
-      error("Error while allocating memory for star particles");
+      error("Error while allocating memory for stars particles");
     bzero(*sparts, *Nstars * sizeof(struct spart));
   }
 
@@ -594,7 +594,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     Ndm = N[1];
     *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
                N[swift_type_dark_matter] +
-               (with_stars ? N[swift_type_star] : 0);
+               (with_stars ? N[swift_type_stars] : 0);
     if (posix_memalign((void*)gparts, gpart_align,
                        *Ngparts * sizeof(struct gpart)) != 0)
       error("Error while allocating memory for gravity particles");
@@ -655,10 +655,10 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
             }
             break;
 
-          case swift_type_star:
+          case swift_type_stars:
             if (with_stars) {
               Nparticles = *Nstars;
-              star_read_particles(*sparts, list, &num_fields);
+              stars_read_particles(*sparts, list, &num_fields);
             }
             break;
 
@@ -700,9 +700,9 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     /* Duplicate the hydro particles into gparts */
     if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
 
-    /* Duplicate the star particles into gparts */
+    /* Duplicate the stars particles into gparts */
     if (with_stars)
-      io_duplicate_star_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+      io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
 
     threadpool_clean(&tp);
   }
@@ -1039,9 +1039,9 @@ void write_output_serial(struct engine* e, const char* baseName,
             darkmatter_write_particles(dmparts, list, &num_fields);
             break;
 
-          case swift_type_star:
+          case swift_type_stars:
             Nparticles = Nstars;
-            star_write_particles(sparts, list, &num_fields);
+            stars_write_particles(sparts, list, &num_fields);
             break;
 
           default:
diff --git a/src/single_io.c b/src/single_io.c
index 2e6ce4111f5258e53dbf4c9e366b534556cb1b05..491dd53d59c58d6157e590641b62dbfad29be981 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -471,10 +471,10 @@ void read_ic_single(const char* fileName,
 
   /* Allocate memory to store star particles */
   if (with_stars) {
-    *Nstars = N[swift_type_star];
+    *Nstars = N[swift_type_stars];
     if (posix_memalign((void**)sparts, spart_align,
                        *Nstars * sizeof(struct spart)) != 0)
-      error("Error while allocating memory for star particles");
+      error("Error while allocating memory for stars particles");
     bzero(*sparts, *Nstars * sizeof(struct spart));
   }
 
@@ -483,7 +483,7 @@ void read_ic_single(const char* fileName,
     Ndm = N[swift_type_dark_matter];
     *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
                N[swift_type_dark_matter] +
-               (with_stars ? N[swift_type_star] : 0);
+               (with_stars ? N[swift_type_stars] : 0);
     if (posix_memalign((void**)gparts, gpart_align,
                        *Ngparts * sizeof(struct gpart)) != 0)
       error("Error while allocating memory for gravity particles");
@@ -532,10 +532,10 @@ void read_ic_single(const char* fileName,
         }
         break;
 
-      case swift_type_star:
+      case swift_type_stars:
         if (with_stars) {
           Nparticles = *Nstars;
-          star_read_particles(*sparts, list, &num_fields);
+          stars_read_particles(*sparts, list, &num_fields);
         }
         break;
 
@@ -568,7 +568,7 @@ void read_ic_single(const char* fileName,
 
     /* Duplicate the star particles into gparts */
     if (with_stars)
-      io_duplicate_star_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
+      io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
 
     threadpool_clean(&tp);
   }
@@ -855,9 +855,9 @@ void write_output_single(struct engine* e, const char* baseName,
         darkmatter_write_particles(dmparts, list, &num_fields);
         break;
 
-      case swift_type_star:
+      case swift_type_stars:
         N = Nstars;
-        star_write_particles(sparts, list, &num_fields);
+        stars_write_particles(sparts, list, &num_fields);
         break;
 
       default:
diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index 65baba7b41c6e112ad6ec26c23c5505f685e10fe..fe4412548305e1b4cee9460d3a71587adf6db8e8 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -60,8 +60,8 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
   sp->num_ngb_density = 0;
 #endif
 
-  sp->wcount = 0.f;
-  sp->wcount_dh = 0.f;
+  sp->density.wcount = 0.f;
+  sp->density.wcount_dh = 0.f;
 }
 
 /**
@@ -107,8 +107,8 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
   const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
 
   /* Finish the calculation by inserting the missing h-factors */
-  sp->wcount *= h_inv_dim;
-  sp->wcount_dh *= h_inv_dim_plus_one;
+  sp->density.wcount *= h_inv_dim;
+  sp->density.wcount_dh *= h_inv_dim_plus_one;
 }
 
 
@@ -127,8 +127,8 @@ __attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   /* Re-set problematic values */
-  sp->wcount = kernel_root * h_inv_dim;
-  sp->wcount_dh = 0.f;
+  sp->density.wcount = kernel_root * h_inv_dim;
+  sp->density.wcount_dh = 0.f;
 }
 
 /**
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index 6bd1c49053d582f816dc97b02daffe9ad23bf640..2aa3d5f3c866044d08da01071a1eed5c65b54931 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -26,8 +26,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_stars_densi
   kernel_deval(ui, &wi, &wi_dx);
 
   /* Compute contribution to the number of neighbours */
-  si->wcount += wi;
-  si->wcount_dh -= (hydro_dimension * wi + ui * wi_dx);    
+  si->density.wcount += wi;
+  si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);    
 
 #ifdef DEBUG_INTERACTIONS_STARS
   /* Update ngb counters */
diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
index 5b113378586bc0f3b3292fbae522e53c825998b0..4fd28a813d735e7af570171fd0be098a84b85517 100644
--- a/src/stars/Default/stars_io.h
+++ b/src/stars/Default/stars_io.h
@@ -90,7 +90,8 @@ INLINE static void stars_props_init(struct stars_props *sp,
 		      const struct hydro_props *p) {
 
   /* Kernel properties */
-  sp->eta_neighbours = parser_get_param_float(params, "Stars:resolution_eta");
+  sp->eta_neighbours = parser_get_opt_param_float(params, "Stars:resolution_eta",
+						  p->eta_neighbours);
 
   /* Tolerance for the smoothing length Newton-Raphson scheme */
   sp->h_tolerance = parser_get_opt_param_float(params, "Stars:h_tolerance",
diff --git a/src/task.c b/src/task.c
index f3ab4ec7c0fa4d9b2621985eaccfac02ee8dca88..2b1b27863a57e5ed1d39189117fa064126d189ff 100644
--- a/src/task.c
+++ b/src/task.c
@@ -57,13 +57,13 @@ const char *taskID_names[task_type_count] = {
     "send",          "recv",          "grav_long_range",
     "grav_mm",       "grav_down_in",  "grav_down",
     "grav_mesh",     "cooling",       "sourceterms",
-    "star_ghost_in", "star_ghost",    "star_ghost_out"};
+    "stars_ghost_in","stars_ghost",   "stars_ghost_out"};
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {
     "none", "density", "gradient", "force", "grav",      "external_grav",
     "tend", "xv",      "rho",      "gpart", "multipole", "spart",
-    "star_density"};
+    "stars_density"};
 
 #ifdef WITH_MPI
 /* MPI communicators for the subtypes. */
@@ -122,7 +122,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
       return task_action_part;
       break;
 
-    case task_type_star_ghost:
+    case task_type_stars_ghost:
       return task_action_spart;
       break;
 
@@ -138,7 +138,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
           return task_action_part;
           break;
 
-        case task_subtype_star_density:
+        case task_subtype_stars_density:
 	  return task_action_all;
 	  break;
 
diff --git a/src/task.h b/src/task.h
index fd7bb985b0894793973b7b4814f8c96c62663eb1..c0793a453e5206701eb0512203d020762fb2e3f6 100644
--- a/src/task.h
+++ b/src/task.h
@@ -66,9 +66,9 @@ enum task_types {
   task_type_grav_mesh,
   task_type_cooling,
   task_type_sourceterms,
-  task_type_star_ghost_in,
-  task_type_star_ghost,
-  task_type_star_ghost_out,
+  task_type_stars_ghost_in,
+  task_type_stars_ghost,
+  task_type_stars_ghost_out,
   task_type_count
 } __attribute__((packed));
 
@@ -88,7 +88,7 @@ enum task_subtypes {
   task_subtype_gpart,
   task_subtype_multipole,
   task_subtype_spart,
-  task_subtype_star_density,
+  task_subtype_stars_density,
   task_subtype_count
 } __attribute__((packed));
 
diff --git a/src/timers.c b/src/timers.c
index b28eccf4d0685195f39d528212a7e0b2eda246c9..51d0e5f6dc4fd9e4e0567592750b8de45ecda06b 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -86,13 +86,13 @@ const char* timers_names[timer_count] = {
     "locktree",
     "runners",
     "step",
-    "doself_star_density",
-    "dopair_star_density",
-    "do_star_ghost",
-    "doself_subset_star_density",
-    "dopair_subset_star_density",
-    "dosubpair_star_density",
-    "dosub_self_star_density",
+    "doself_stars_density",
+    "dopair_stars_density",
+    "do_stars_ghost",
+    "doself_subset_stars_density",
+    "dopair_subset_stars_density",
+    "dosubpair_stars_density",
+    "dosub_self_stars_density",
 };
 
 /* File to store the timers */
diff --git a/src/timers.h b/src/timers.h
index 04e2d441c4f8864402af849ab2bc2d917031f312..aba6ae33e3b8268c088694863967b96851715153 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -87,13 +87,13 @@ enum {
   timer_locktree,
   timer_runners,
   timer_step,
-  timer_doself_star_density,
-  timer_dopair_star_density,
-  timer_dostar_ghost,
-  timer_doself_subset_star_density,
-  timer_dopair_subset_star_density,
-  timer_dosubpair_star_density,
-  timer_dosub_self_star_density,
+  timer_doself_stars_density,
+  timer_dopair_stars_density,
+  timer_dostars_ghost,
+  timer_doself_subset_stars_density,
+  timer_dopair_subset_stars_density,
+  timer_dosubpair_stars_density,
+  timer_dosub_self_stars_density,
   timer_count,
 };
 
diff --git a/src/timestep.h b/src/timestep.h
index d065df4c444cb880a74688be97245c285a391817..02e32cd8d1eb37215101398f1fd884ebd7701018 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -181,7 +181,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
     const struct spart *restrict sp, const struct engine *restrict e) {
 
   /* Stellar time-step */
-  float new_dt_star = star_compute_timestep(sp);
+  float new_dt_stars = stars_compute_timestep(sp);
 
   /* Gravity time-step */
   float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
@@ -196,7 +196,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
         sp->gpart, a_hydro, e->gravity_properties, e->cosmology);
 
   /* Take the minimum of all */
-  float new_dt = min3(new_dt_star, new_dt_self, new_dt_ext);
+  float new_dt = min3(new_dt_stars, new_dt_self, new_dt_ext);
 
   /* Apply the maximal displacement constraint (FLT_MAX  if non-cosmological)*/
   new_dt = min(new_dt, e->dt_max_RMS_displacement);
diff --git a/src/tools.c b/src/tools.c
index f3126300769a60edcebc0b456f6d1cb62ce0a1c1..45b9e90799ad28c15b0ed91172907d9797eb7f5b 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -336,7 +336,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 }
 
-void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj) {
+void pairs_all_stars_density(struct runner *r, struct cell *ci, struct cell *cj) {
 
   float r2, dx[3];
   const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
@@ -370,7 +370,7 @@ void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj)
       /* Hit or miss? */
       if (r2 < hig2) {
         /* Interact */
-        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);
       }
     }
   }
@@ -400,7 +400,7 @@ void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj)
       /* Hit or miss? */
       if (r2 < hjg2) {
         /* Interact */
-        runner_iact_nonsym_star_density(r2, dx, hj, pi->h, spj, pi, a, H);
+        runner_iact_nonsym_stars_density(r2, dx, hj, pi->h, spj, pi, a, H);
       }
     }
   }
@@ -499,7 +499,7 @@ void self_all_force(struct runner *r, struct cell *ci) {
   }
 }
 
-void self_all_star_density(struct runner *r, struct cell *ci) {
+void self_all_stars_density(struct runner *r, struct cell *ci) {
   float r2, hi, hj, hig2, dxi[3];
   struct spart *spi;
   struct part *pj;
@@ -533,7 +533,7 @@ void self_all_star_density(struct runner *r, struct cell *ci) {
       /* Hit or miss? */
       if (r2 > 0.f && r2 < hig2) {
         /* Interact */
-        runner_iact_nonsym_star_density(r2, dxi, hi, hj, spi, pj, a, H);
+        runner_iact_nonsym_stars_density(r2, dxi, hi, hj, spi, pj, a, H);
       }
 
     }
diff --git a/src/tools.h b/src/tools.h
index cb49b5cb8a1811e9f926a06b35e3e4caa3e859db..1781919b8b1966ffa451a02e3efe385b8c2a66a5 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -40,8 +40,8 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_density(struct runner *r, struct cell *ci);
 void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_force(struct runner *r, struct cell *ci);
-void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj);
-void self_all_star_density(struct runner *r, struct cell *ci);
+void pairs_all_stars_density(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_stars_density(struct runner *r, struct cell *ci);
 
 void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
 
diff --git a/tests/test27cellsStars.c b/tests/test27cellsStars.c
index cb5fd0397d3bb904e01f90d34a9930da99247763..fb84515048e017ec90fa136fb883da1293fb0ff4 100644
--- a/tests/test27cellsStars.c
+++ b/tests/test27cellsStars.c
@@ -30,20 +30,20 @@
 /* Local headers. */
 #include "swift.h"
 
-#define DOSELF1 runner_doself_branch_star_density
-#define DOSELF1_SUBSET runner_doself_subset_branch_star_density
+#define DOSELF1 runner_doself_branch_stars_density
+#define DOSELF1_SUBSET runner_doself_subset_branch_stars_density
 #ifdef TEST_DOSELF_SUBSET
-#define DOSELF1_NAME "runner_doself_subset_branch_star_density"
+#define DOSELF1_NAME "runner_doself_subset_branch_stars_density"
 #else
-#define DOSELF1_NAME "runner_doself1_branch_star_density"
+#define DOSELF1_NAME "runner_doself1_branch_stars_density"
 #endif
 
-#define DOPAIR1 runner_dopair_branch_star_density
-#define DOPAIR1_SUBSET runner_dopair_subset_branch_star_density
+#define DOPAIR1 runner_dopair_branch_stars_density
+#define DOPAIR1_SUBSET runner_dopair_subset_branch_stars_density
 #ifdef TEST_DOPAIR_SUBSET
-#define DOPAIR1_NAME "runner_dopair_subset_branch_star_density"
+#define DOPAIR1_NAME "runner_dopair_subset_branch_stars_density"
 #else
-#define DOPAIR1_NAME "runner_dopair_branch_star_density"
+#define DOPAIR1_NAME "runner_dopair_branch_stars_density"
 #endif
 
 #define NODE_ID 0
@@ -63,10 +63,10 @@
  * of the inter-particle separation.
  * @param h_pert The perturbation to apply to the smoothing length.
  */
-struct cell *make_cell(size_t n, size_t n_star, double *offset, double size, double h,
+struct cell *make_cell(size_t n, size_t n_stars, double *offset, double size, double h,
                        long long *partId, long long *spartId, double pert, double h_pert) {
   const size_t count = n * n * n;
-  const size_t scount = n_star * n_star * n_star;
+  const size_t scount = n_stars * n_stars * n_stars;
   float h_max = 0.f;
   struct cell *cell = (struct cell *)malloc(sizeof(struct cell));
   bzero(cell, sizeof(struct cell));
@@ -121,26 +121,26 @@ struct cell *make_cell(size_t n, size_t n_star, double *offset, double size, dou
   bzero(cell->sparts, scount * sizeof(struct spart));
 
   struct spart *spart = cell->sparts;
-  for (size_t x = 0; x < n_star; ++x) {
-    for (size_t y = 0; y < n_star; ++y) {
-      for (size_t z = 0; z < n_star; ++z) {
+  for (size_t x = 0; x < n_stars; ++x) {
+    for (size_t y = 0; y < n_stars; ++y) {
+      for (size_t z = 0; z < n_stars; ++z) {
         spart->x[0] =
             offset[0] +
-            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_stars;
         spart->x[1] =
             offset[1] +
-            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_stars;
         spart->x[2] =
             offset[2] +
-            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_stars;
 
 	spart->v[0] = 0;
 	spart->v[1] = 0;
 	spart->v[2] = 0;
 	if (h_pert)
-	  spart->h = size * h * random_uniform(1.f, h_pert) / (float)n_star;
+	  spart->h = size * h * random_uniform(1.f, h_pert) / (float)n_stars;
 	else
-	  spart->h = size * h / (float)n_star;
+	  spart->h = size * h / (float)n_stars;
 	h_max = fmaxf(h_max, spart->h);
 	spart->id = ++(*spartId);
 
@@ -196,7 +196,7 @@ void clean_up(struct cell *ci) {
  */
 void zero_particle_fields(struct cell *c) {
   for (int pid = 0; pid < c->scount; pid++) {
-    star_init_spart(&c->sparts[pid]);
+    stars_init_spart(&c->sparts[pid]);
   }
 }
 
@@ -205,11 +205,11 @@ void zero_particle_fields(struct cell *c) {
  */
 void end_calculation(struct cell *c, const struct cosmology *cosmo) {
   for (int pid = 0; pid < c->scount; pid++) {
-    star_end_density(&c->sparts[pid], cosmo);
+    stars_end_density(&c->sparts[pid], cosmo);
 
     /* Recover the common "Neighbour number" definition */
-    c->sparts[pid].wcount *= pow_dimension(c->sparts[pid].h);
-    c->sparts[pid].wcount *= kernel_norm;
+    c->sparts[pid].density.wcount *= pow_dimension(c->sparts[pid].h);
+    c->sparts[pid].density.wcount *= kernel_norm;
   }
 }
 
@@ -233,8 +233,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             "%6llu %10f %10f %10f %13e %13e\n",
             main_cell->sparts[pid].id, main_cell->sparts[pid].x[0],
             main_cell->sparts[pid].x[1], main_cell->sparts[pid].x[2],
-            main_cell->sparts[pid].wcount,
-            main_cell->sparts[pid].wcount_dh);
+            main_cell->sparts[pid].density.wcount,
+            main_cell->sparts[pid].density.wcount_dh);
   }
 
   /* Write all other cells */
@@ -254,7 +254,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               "%6llu %10f %10f %10f %13e %13e\n",
               cj->sparts[pjd].id, cj->sparts[pjd].x[0], cj->sparts[pjd].x[1],
               cj->sparts[pjd].x[2],
-              cj->sparts[pjd].wcount, cj->sparts[pjd].wcount_dh);
+              cj->sparts[pjd].density.wcount, cj->sparts[pjd].density.wcount_dh);
         }
       }
     }
@@ -263,15 +263,15 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
 }
 
 /* Just a forward declaration... */
-void runner_dopair_branch_star_density(struct runner *r, struct cell *ci,
+void runner_dopair_branch_stars_density(struct runner *r, struct cell *ci,
                                    struct cell *cj);
-void runner_doself_branch_star_density(struct runner *r, struct cell *c);
-void runner_dopair_subset_branch_star_density(struct runner *r,
+void runner_doself_branch_stars_density(struct runner *r, struct cell *c);
+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);
-void runner_doself_subset_branch_star_density(struct runner *r,
+void runner_doself_subset_branch_stars_density(struct runner *r,
                                          struct cell *restrict ci,
                                          struct spart *restrict sparts,
                                          int *restrict ind, int scount);
@@ -339,8 +339,8 @@ int main(int argc, char *argv[]) {
     printf(
         "\nUsage: %s -n PARTICLES_PER_AXIS -N SPARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
         "\nGenerates 27 cells, filled with particles on a Cartesian grid."
-        "\nThese are then interacted using runner_dopair_star_density() and "
-        "runner_doself_star_density()."
+        "\nThese are then interacted using runner_dopair_stars_density() and "
+        "runner_doself_stars_density()."
         "\n\nOptions:"
         "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
         "\n-p                 - Random fractional change in h, h=h*random(1,p)"
@@ -374,11 +374,11 @@ int main(int argc, char *argv[]) {
   hp.max_smoothing_iterations = 1;
   hp.CFL_condition = 0.1;
 
-  struct stars_props star_p;
-  star_p.eta_neighbours = h;
-  star_p.h_tolerance = 1e0;
-  star_p.h_max = FLT_MAX;
-  star_p.max_smoothing_iterations = 1;
+  struct stars_props stars_p;
+  stars_p.eta_neighbours = h;
+  stars_p.h_tolerance = 1e0;
+  stars_p.h_max = FLT_MAX;
+  stars_p.max_smoothing_iterations = 1;
 
   struct engine engine;
   engine.s = &space;
@@ -386,7 +386,7 @@ int main(int argc, char *argv[]) {
   engine.ti_current = 8;
   engine.max_active_bin = num_time_bins;
   engine.hydro_properties = &hp;
-  engine.stars_properties = &star_p;
+  engine.stars_properties = &stars_p;
   engine.nodeID = NODE_ID;
 
   struct cosmology cosmo;
@@ -508,10 +508,10 @@ int main(int argc, char *argv[]) {
 
   /* Run all the brute-force pairs */
   for (int j = 0; j < 27; ++j)
-    if (cells[j] != main_cell) pairs_all_star_density(&runner, main_cell, cells[j]);
+    if (cells[j] != main_cell) pairs_all_stars_density(&runner, main_cell, cells[j]);
 
   /* And now the self-interaction */
-  self_all_star_density(&runner, main_cell);
+  self_all_stars_density(&runner, main_cell);
 
   const ticks toc = getticks();