diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 2a5aeba7d1db0b1e1e56a9a6eed3059aba6a09ff..0df1f91194b6d1e7e98cb1b75be7d3eaaca7fc32 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -762,6 +762,7 @@ WARN_LOGFILE           =
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
 INPUT		       += @top_srcdir@/src/hydro/Minimal
 INPUT		       += @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
 INPUT		       += @top_srcdir@/src/cooling/const_du
diff --git a/examples/main.c b/examples/main.c
index 11163b42523fa5b1de1438ad8e67dde0fe9c88ef..d8049c3a505002b5c840c61b6b1898accb335e83 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -303,6 +303,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) {
     message("sizeof(struct part)  is %4zi bytes.", sizeof(struct part));
     message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(struct spart) is %4zi bytes.", sizeof(struct spart));
     message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
     message("sizeof(struct task)  is %4zi bytes.", sizeof(struct task));
     message("sizeof(struct cell)  is %4zi bytes.", sizeof(struct cell));
@@ -368,7 +369,8 @@ int main(int argc, char *argv[]) {
 
   struct part *parts = NULL;
   struct gpart *gparts = NULL;
-  size_t Ngas = 0, Ngpart = 0;
+  struct spart *sparts = NULL;
+  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
   double dim[3] = {0., 0., 0.};
   int periodic = 0;
   int flag_entropy_ICs = 0;
@@ -384,8 +386,8 @@ int main(int argc, char *argv[]) {
                  MPI_INFO_NULL, dry_run);
 #endif
 #else
-  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
-                 &periodic, &flag_entropy_ICs, dry_run);
+  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
+                 &Nspart, &periodic, &flag_entropy_ICs, dry_run);
 #endif
   if (myrank == 0) {
     clocks_gettime(&toc);
@@ -400,6 +402,7 @@ int main(int argc, char *argv[]) {
     free(gparts);
     gparts = NULL;
     for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
+    for (size_t k = 0; k < Nspart; ++k) sparts[k].gpart = NULL;
     Ngpart = 0;
   }
   if (!with_hydro) {
@@ -411,23 +414,26 @@ int main(int argc, char *argv[]) {
   }
 
   /* Get the total number of particles across all nodes. */
-  long long N_total[2] = {0, 0};
+  long long N_total[3] = {0, 0, 0};
 #if defined(WITH_MPI)
   long long N_long[2] = {Ngas, Ngpart};
   MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
 #else
   N_total[0] = Ngas;
   N_total[1] = Ngpart;
+  N_total[2] = Nspart;
 #endif
   if (myrank == 0)
-    message("Read %lld gas particles and %lld gparts from the ICs.", N_total[0],
-            N_total[1]);
+    message(
+        "Read %lld gas particles, %lld star particles and %lld gparts from the "
+        "ICs.",
+        N_total[0], N_total[2], N_total[1]);
 
   /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
   struct space s;
-  space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic,
-             with_self_gravity, talking, dry_run);
+  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
+             periodic, with_self_gravity, talking, dry_run);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -529,6 +535,8 @@ int main(int argc, char *argv[]) {
     return 0;
   }
 
+  engine_dump_snapshot(&e);
+
 #ifdef WITH_MPI
   /* Split the space. */
   engine_split(&e, &initial_partition);
diff --git a/src/common_io.c b/src/common_io.c
index 1f1ec401547c81e137b4e7d836ab58cb87280d8b..a9981e3024d76c445a8d69b03512a67072296c05 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -597,7 +597,7 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
  *
  * @param parts The array of #part freshly read in.
  * @param gparts The array of #gpart freshly read in with all the DM particles
- *at the start
+ * at the start
  * @param Ngas The number of gas particles read in.
  * @param Ndm The number of DM particles read in.
  */
@@ -624,6 +624,41 @@ void duplicate_hydro_gparts(struct part* const parts,
   }
 }
 
+/**
+ * @brief Copy every #spart into the corresponding #gpart and link them.
+ *
+ * This function assumes that the DM particles and gas particles are all at
+ * the start of the gparts array and adds the star particles afterwards
+ *
+ * @param sparts The array of #spart freshly read in.
+ * @param gparts The array of #gpart freshly read in with all the DM and gas
+ * particles at the start.
+ * @param Nstars The number of stars particles read in.
+ * @param Ndm The number of DM and gas particles read in.
+ */
+void duplicate_star_gparts(struct spart* const sparts,
+                           struct gpart* const gparts, size_t Nstars,
+                           size_t Ndm) {
+
+  for (size_t i = 0; i < Nstars; ++i) {
+
+    /* Duplicate the crucial information */
+    gparts[i + Ndm].x[0] = sparts[i].x[0];
+    gparts[i + Ndm].x[1] = sparts[i].x[1];
+    gparts[i + Ndm].x[2] = sparts[i].x[2];
+
+    gparts[i + Ndm].v_full[0] = sparts[i].v[0];
+    gparts[i + Ndm].v_full[1] = sparts[i].v[1];
+    gparts[i + Ndm].v_full[2] = sparts[i].v[2];
+
+    gparts[i + Ndm].mass = sparts[i].mass;
+
+    /* Link the particles */
+    gparts[i + Ndm].id_or_neg_offset = -i;
+    sparts[i].gpart = &gparts[i + Ndm];
+  }
+}
+
 /**
  * @brief Copy every DM #gpart into the dmparts array.
  *
diff --git a/src/common_io.h b/src/common_io.h
index 7aedee0f2624dcff916a8398e244009a87109915..bf1840d497c46f58568d1bed7cb3409f60e047ee 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -75,6 +75,9 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm);
 void duplicate_hydro_gparts(struct part* const parts,
                             struct gpart* const gparts, size_t Ngas,
                             size_t Ndm);
+void duplicate_star_gparts(struct spart* const sparts,
+                           struct gpart* const gparts, size_t Nstars,
+                           size_t Ndm);
 
 void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
 
diff --git a/src/part.h b/src/part.h
index 0bf4359f891619b0900f8aa9f17b2a2a71127579..269e12a9c6b4adbd1046f1a008f1b21ab3214128 100644
--- a/src/part.h
+++ b/src/part.h
@@ -36,6 +36,7 @@
 /* Some constants. */
 #define part_align 128
 #define xpart_align 128
+#define spart_align 128
 #define gpart_align 128
 
 /* Import the right hydro particle definition */
@@ -62,6 +63,9 @@
 /* Import the right gravity particle definition */
 #include "./gravity/Default/gravity_part.h"
 
+/* Import the right star particle definition */
+#include "./stars/Default/star_part.h"
+
 void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset);
 void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts);
 #ifdef WITH_MPI
diff --git a/src/single_io.c b/src/single_io.c
index ceeba4eb80a47c3feed7e898deb5e1fe7e427c0c..4b78eeb5accc49f8f8de2be800eab9268e77a410 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -45,6 +45,7 @@
 #include "io_properties.h"
 #include "kernel_hydro.h"
 #include "part.h"
+#include "stars_io.h"
 #include "units.h"
 
 /*-----------------------------------------------------------------------------
@@ -312,10 +313,12 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
  * @param fileName The file to read.
  * @param internal_units The system units used internally
  * @param dim (output) The dimension of the volume.
- * @param parts (output) Array of Gas particles.
+ * @param parts (output) Array of #part particles.
  * @param gparts (output) Array of #gpart particles.
+ * @param sparts (output) Array of #spart particles.
  * @param Ngas (output) number of Gas particles read.
  * @param Ngparts (output) The number of #gpart read.
+ * @param Nstars (output) The number of #spart read.
  * @param periodic (output) 1 if the volume is periodic, 0 if not.
  * @param flag_entropy (output) 1 if the ICs contained Entropy in the
  * InternalEnergy
@@ -332,8 +335,9 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
  */
 void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
-                    size_t* Ngas, size_t* Ngparts, int* periodic,
-                    int* flag_entropy, int dry_run) {
+                    struct spart** sparts, size_t* Ngas, size_t* Ngparts,
+                    size_t* Nstars, int* periodic, int* flag_entropy,
+                    int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -439,15 +443,22 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
         units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
 
   /* Allocate memory to store SPH particles */
-  *Ngas = N[0];
+  *Ngas = N[GAS];
   if (posix_memalign((void*)parts, part_align, *Ngas * sizeof(struct part)) !=
       0)
     error("Error while allocating memory for SPH particles");
   bzero(*parts, *Ngas * sizeof(struct part));
 
+  /* Allocate memory to store star particles */
+  *Nstars = N[STAR];
+  if (posix_memalign((void*)sparts, spart_align,
+                     *Nstars * sizeof(struct spart)) != 0)
+    error("Error while allocating memory for star particles");
+  bzero(*sparts, *Nstars * sizeof(struct spart));
+
   /* Allocate memory to store all particles */
-  Ndm = N[1];
-  *Ngparts = N[1] + N[0];
+  Ndm = N[DM];
+  *Ngparts = N[GAS] + N[DM] + N[STAR];
   if (posix_memalign((void*)gparts, gpart_align,
                      *Ngparts * sizeof(struct gpart)) != 0)
     error("Error while allocating memory for gravity particles");
@@ -491,6 +502,11 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
         darkmatter_read_particles(*gparts, list, &num_fields);
         break;
 
+      case STAR:
+        Nparticles = *Nstars;
+        star_read_particles(*sparts, list, &num_fields);
+        break;
+
       default:
         message("Particle Type %d not yet supported. Particles ignored", ptype);
     }
@@ -507,9 +523,12 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
   /* Prepare the DM particles */
   if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
 
-  /* Now duplicate the hydro particle into gparts */
+  /* Duplicate the hydro particles into gparts */
   if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
+  /* Duplicate the star particles into gparts */
+  if (!dry_run) duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
+
   /* message("Done Reading particles..."); */
 
   /* Clean up */
@@ -541,18 +560,20 @@ void write_output_single(struct engine* e, const char* baseName,
 
   hid_t h_file = 0, h_grp = 0;
   const size_t Ngas = e->s->nr_parts;
+  const size_t Nstars = e->s->nr_sparts;
   const size_t Ntot = e->s->nr_gparts;
   int periodic = e->s->periodic;
   int numFiles = 1;
   struct part* parts = e->s->parts;
   struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
+  struct spart* sparts = e->s->sparts;
   static int outputCount = 0;
 
   /* Number of unassociated gparts */
-  const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
+  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
 
-  long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
+  long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0, 0, Nstars, 0};
 
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
@@ -729,6 +750,11 @@ void write_output_single(struct engine* e, const char* baseName,
         darkmatter_write_particles(dmparts, list, &num_fields);
         break;
 
+      case STAR:
+        N = Nstars;
+        star_write_particles(sparts, list, &num_fields);
+        break;
+
       default:
         error("Particle Type %d not yet supported. Aborting", ptype);
     }
diff --git a/src/single_io.h b/src/single_io.h
index 51a30a7bc6af7f3aaf5708a3d2df14982e026e3e..dc60b86b162188c4a6d15f0aeb9758ea8a098984 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -31,7 +31,8 @@
 
 void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
-                    size_t* Ngas, size_t* Ndm, int* periodic, int* flag_entropy,
+                    struct spart** sparts, size_t* Ngas, size_t* Ndm,
+                    size_t* Nstars, int* periodic, int* flag_entropy,
                     int dry_run);
 
 void write_output_single(struct engine* e, const char* baseName,
diff --git a/src/space.c b/src/space.c
index 935677a9ebed97acfde8341ec1545ef4f33a56c0..eebef1a0f99fd4e5d22f7875530c4c64ea23ae1f 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1770,8 +1770,9 @@ void space_init_gparts(struct space *s) {
  */
 void space_init(struct space *s, const struct swift_params *params,
                 double dim[3], struct part *parts, struct gpart *gparts,
-                size_t Npart, size_t Ngpart, int periodic, int gravity,
-                int verbose, int dry_run) {
+                struct spart *sparts, size_t Npart, size_t Ngpart,
+                size_t Nspart, int periodic, int gravity, int verbose,
+                int dry_run) {
 
   /* Clean-up everything */
   bzero(s, sizeof(struct space));
@@ -1789,6 +1790,9 @@ void space_init(struct space *s, const struct swift_params *params,
   s->nr_gparts = Ngpart;
   s->size_gparts = Ngpart;
   s->gparts = gparts;
+  s->nr_sparts = Nspart;
+  s->size_sparts = Nspart;
+  s->sparts = sparts;
   s->nr_queues = 1; /* Temporary value until engine construction */
 
   /* Decide on the minimal top-level cell size */
@@ -1848,6 +1852,11 @@ void space_init(struct space *s, const struct swift_params *params,
       gparts[k].x[1] += shift[1];
       gparts[k].x[2] += shift[2];
     }
+    for (size_t k = 0; k < Nspart; k++) {
+      sparts[k].x[0] += shift[0];
+      sparts[k].x[1] += shift[1];
+      sparts[k].x[2] += shift[2];
+    }
   }
 
   if (!dry_run) {
@@ -1879,9 +1888,23 @@ void space_init(struct space *s, const struct swift_params *params,
           if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
             error("Not all g-particles are within the specified domain.");
     }
+
+    /* Same for the sparts */
+    if (periodic) {
+      for (size_t k = 0; k < Nspart; k++)
+        for (int j = 0; j < 3; j++) {
+          while (sparts[k].x[j] < 0) sparts[k].x[j] += dim[j];
+          while (sparts[k].x[j] >= dim[j]) sparts[k].x[j] -= dim[j];
+        }
+    } else {
+      for (size_t k = 0; k < Nspart; k++)
+        for (int j = 0; j < 3; j++)
+          if (sparts[k].x[j] < 0 || sparts[k].x[j] >= dim[j])
+            error("Not all s-particles are within the specified domain.");
+    }
   }
 
-  /* Allocate the extra parts array. */
+  /* Allocate the extra parts array for the gas particles. */
   if (Npart > 0) {
     if (posix_memalign((void *)&s->xparts, xpart_align,
                        Npart * sizeof(struct xpart)) != 0)
@@ -1893,6 +1916,7 @@ void space_init(struct space *s, const struct swift_params *params,
   space_init_parts(s);
   space_init_xparts(s);
   space_init_gparts(s);
+  // space_init_sparts(s);  // MATTHIEU
 
   /* Init the space lock. */
   if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
diff --git a/src/space.h b/src/space.h
index 53cf2d0c8fa548ae19aa7452abb38c3e3e028165..17f16ecff3a313464752bb5c5f050cbe2e27b0f4 100644
--- a/src/space.h
+++ b/src/space.h
@@ -108,6 +108,9 @@ struct space {
   /*! The total number of g-parts in the space. */
   size_t nr_gparts, size_gparts;
 
+  /*! The total number of g-parts in the space. */
+  size_t nr_sparts, size_sparts;
+
   /*! The particle data (cells have pointers to this). */
   struct part *parts;
 
@@ -117,6 +120,9 @@ struct space {
   /*! The g-particle data (cells have pointers to this). */
   struct gpart *gparts;
 
+  /*! The s-particle data (cells have pointers to this). */
+  struct spart *sparts;
+
   /*! General-purpose lock for this space. */
   swift_lock_type lock;
 
@@ -152,8 +158,9 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
                  double *shift);
 void space_init(struct space *s, const struct swift_params *params,
                 double dim[3], struct part *parts, struct gpart *gparts,
-                size_t Npart, size_t Ngpart, int periodic, int gravity,
-                int verbose, int dry_run);
+                struct spart *sparts, size_t Npart, size_t Ngpart,
+                size_t Nspart, int periodic, int gravity, int verbose,
+                int dry_run);
 void space_sanitize(struct space *s);
 void space_map_cells_pre(struct space *s, int full,
                          void (*fun)(struct cell *c, void *data), void *data);
diff --git a/src/stars/Default/star.h b/src/stars/Default/star.h
new file mode 100644
index 0000000000000000000000000000000000000000..9e0ca81edff06b8a32afb185f24a88b41dc87da7
--- /dev/null
+++ b/src/stars/Default/star.h
@@ -0,0 +1,105 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *               2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_H
+#define SWIFT_DEFAULT_GRAVITY_H
+
+#include <float.h>
+#include "minmax.h"
+
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct gpart* const gp) {
+
+  const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
+                    gp->a_grav[1] * gp->a_grav[1] +
+                    gp->a_grav[2] * gp->a_grav[2];
+
+  const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
+
+  const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
+
+  return dt;
+}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+    struct gpart* gp) {
+
+  gp->ti_begin = 0;
+  gp->ti_end = 0;
+  gp->epsilon = 0.;  // MATTHIEU
+}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_gpart(
+    struct gpart* gp) {
+
+  /* Zero the acceleration */
+  gp->a_grav[0] = 0.f;
+  gp->a_grav[1] = 0.f;
+  gp->a_grav[2] = 0.f;
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ * @param const_G Newton's constant in internal units
+ */
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+    struct gpart* gp, float const_G) {
+
+  /* Let's get physical... */
+  gp->a_grav[0] *= const_G;
+  gp->a_grav[1] *= const_G;
+  gp->a_grav[2] *= const_G;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+    struct gpart* gp, float dt, float half_dt) {}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_H */
diff --git a/src/stars/Default/star_debug.h b/src/stars/Default/star_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..fdef8ee0e65231f1fbeb6cf508494065bef4af2e
--- /dev/null
+++ b/src/stars/Default/star_debug.h
@@ -0,0 +1,31 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_STAR_DEBUG_H
+#define SWIFT_DEFAULT_STAR_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void star_debug_particle(
+    const struct spart* p) {
+  printf(
+      "x=[%.3e,%.3e,%.3e], "
+      "v_full=[%.3e,%.3e,%.3e] \n t_begin=%d, t_end=%d\n",
+      p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
+      p->ti_begin, p->ti_end);
+}
+
+#endif /* SWIFT_DEFAULT_STAR_DEBUG_H */
diff --git a/src/stars/Default/star_io.h b/src/stars/Default/star_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..96bbdce6d83dc241d05e7dd1754f476dc0b8e5f9
--- /dev/null
+++ b/src/stars/Default/star_io.h
@@ -0,0 +1,72 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_STAR_IO_H
+#define SWIFT_DEFAULT_STAR_IO_H
+
+#include "io_properties.h"
+
+/**
+ * @brief Specifies which s-particle fields to read from a dataset
+ *
+ * @param sparts The s-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void star_read_particles(struct spart* sparts, struct io_props* list,
+                         int* num_fields) {
+
+  /* Say how much we want to read */
+  *num_fields = 4;
+
+  /* List what we want to read */
+  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+                                UNIT_CONV_LENGTH, sparts, x);
+  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+                                UNIT_CONV_SPEED, sparts, v);
+  list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+                                sparts, mass);
+  list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, sparts, id);
+}
+
+/**
+ * @brief Specifies which s-particle fields to write to a dataset
+ *
+ * @param sparts The s-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+void star_write_particles(struct spart* sparts, struct io_props* list,
+                          int* num_fields) {
+
+  /* Say how much we want to read */
+  *num_fields = 4;
+
+  /* List what we want to read */
+  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+                                 sparts, x);
+  list[1] =
+      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
+  list[2] =
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
+  list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
+                                 sparts, id);
+}
+
+#endif /* SWIFT_DEFAULT_STAR_IO_H */
diff --git a/src/stars/Default/star_part.h b/src/stars/Default/star_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ee5a9452549f5f02aaf8131204b5018ccb1d029
--- /dev/null
+++ b/src/stars/Default/star_part.h
@@ -0,0 +1,53 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_STAR_PART_H
+#define SWIFT_DEFAULT_STAR_PART_H
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+/* Star particle. */
+struct spart {
+
+  /* Particle ID. */
+  long long id;
+
+  /* Pointer to corresponding gravity part. */
+  struct gpart* gpart;
+
+  /* Particle position. */
+  double x[3];
+
+  /* Offset between current position and position at last tree rebuild. */
+  float x_diff[3];
+
+  /* Particle velocity. */
+  float v[3];
+
+  float mass;
+
+  /* Particle time of beginning of time-step. */
+  int ti_begin;
+
+  /* Particle time of end of time-step. */
+  int ti_end;
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_DEFAULT_STAR_PART_H */
diff --git a/src/stars_io.h b/src/stars_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..18a13ec19163008f1c8e9f64cf544ddf812db655
--- /dev/null
+++ b/src/stars_io.h
@@ -0,0 +1,26 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_STAR_IO_H
+#define SWIFT_STAR_IO_H
+
+#include "./const.h"
+
+#include "./stars/Default/star_io.h"
+
+#endif /* SWIFT_STAR_IO_H */