diff --git a/configure.ac b/configure.ac
index 9cadad241f37e86e2020a4e460761ea7c6253414..0162b89d3076dd581ad301878a802af73ab03d41 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1510,6 +1510,9 @@ case "$with_feedback" in
thermal)
AC_DEFINE([FEEDBACK_THERMAL], [1], [Thermal Blastwave])
;;
+ const)
+ AC_DEFINE([FEEDBACK_CONST], [1], [Constant feedback model])
+ ;;
none)
;;
diff --git a/src/Makefile.am b/src/Makefile.am
index 889b0a3affa32303e408c2bc025663fb54b0b4ff..cd9accb47be9b7725002dfcdcd62d980abaf1f86 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -136,6 +136,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
stars.h stars_io.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 \
+ stars/const/stars.h stars/const/stars_iact.h stars/const/stars_io.h \
+ stars/const/stars_debug.h stars/const/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/stars.h b/src/stars.h
index 3e921239a29d862aba998c138623eb1cb81a37b9..a34f49e850e2f498186c1123dac5e3e0ea90b8fb 100644
--- a/src/stars.h
+++ b/src/stars.h
@@ -24,7 +24,12 @@
/* So far only one model here */
/* Straight-forward import */
+#if defined(FEEDBACK_CONST)
+#include "./stars/const/stars.h"
+#include "./stars/const/stars_iact.h"
+#else
#include "./stars/Default/stars.h"
#include "./stars/Default/stars_iact.h"
+#endif
#endif
diff --git a/src/stars/const/stars.h b/src/stars/const/stars.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e54cf7cb48342a274f8b8709c51864f2c2852f8
--- /dev/null
+++ b/src/stars/const/stars.h
@@ -0,0 +1,150 @@
+/*******************************************************************************
+ * 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_CONST_STARS_H
+#define SWIFT_CONST_STARS_H
+
+#include <float.h>
+#include "minmax.h"
+
+/**
+ * @brief Computes the gravity time-step of a given star particle.
+ *
+ * @param sp Pointer to the s-particle data.
+ */
+__attribute__((always_inline)) INLINE static float stars_compute_timestep(
+ const struct spart* const sp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Initialises the s-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_first_init_spart(
+ struct spart* sp) {
+
+ sp->time_bin = 0;
+}
+
+/**
+ * @brief Prepares a s-particle for its interactions
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_init_spart(
+ struct spart* sp) {
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ sp->ids_ngbs_density[i] = -1;
+ sp->num_ngb_density = 0;
+#endif
+
+ sp->density.wcount = 0.f;
+ sp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param sp The particle.
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
+ struct spart* restrict sp) {}
+
+/**
+ * @brief Finishes the calculation of (non-gravity) forces acting on stars
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_end_force(
+ struct spart* sp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param sp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void stars_kick_extra(
+ struct spart* sp, float dt) {}
+
+/**
+ * @brief Finishes the calculation of density on stars
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void stars_end_density(
+ struct spart* sp, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = sp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ 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->density.wcount *= h_inv_dim;
+ sp->density.wcount_dh *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #spart has 0
+ * ngbs.
+ *
+ * @param sp The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
+ struct spart* restrict sp, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = sp->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ sp->density.wcount = kernel_root * h_inv_dim;
+ sp->density.wcount_dh = 0.f;
+}
+
+/**
+ * @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_evolve_spart(
+ struct spart* restrict sp, const struct stars_props* stars_properties,
+ const struct cosmology* cosmo) {}
+
+#endif /* SWIFT_CONST_STARS_H */
diff --git a/src/stars/const/stars_debug.h b/src/stars/const/stars_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..53a54b028ad92fce26a41ba46569db745a51544d
--- /dev/null
+++ b/src/stars/const/stars_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_CONST_STARS_DEBUG_H
+#define SWIFT_CONST_STARS_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void stars_debug_particle(
+ const struct spart* p) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v_full=[%.3e,%.3e,%.3e] p->mass=%.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->mass, p->ti_begin, p->ti_end);
+}
+
+#endif /* SWIFT_CONST_STARS_DEBUG_H */
diff --git a/src/stars/const/stars_iact.h b/src/stars/const/stars_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..59751d33c3c3c8a77cb038da48a54323919c948c
--- /dev/null
+++ b/src/stars/const/stars_iact.h
@@ -0,0 +1,136 @@
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param si First sparticle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
+ struct spart *restrict si,
+ const struct part *restrict pj, float a,
+ float H) {
+
+ float wi, wi_dx;
+
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Compute the kernel function */
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ si->density.wcount += wi;
+ si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ /* Update ngb counters */
+ if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
+ si->ids_ngbs_density[si->num_ngb_density] = pj->id;
+ ++si->num_ngb_density;
+#endif
+}
+
+/**
+ * @brief Feedback interaction between two particles (non-symmetric).
+ * Used for updating properties of gas particles neighbouring a star particle
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (si - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param si First (star) particle.
+ * @param pj Second (gas) particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
+ const struct spart *restrict si,
+ struct part *restrict pj, float a, float H, float omega_normalisation) {
+ float wi, wi_dx;
+
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Compute the kernel function */
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ // ALEXEI: is just plain kernel_eval fine here?
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute weighting for distributing various properties (TODO: better comment?) */
+ // ALEXEI: come up with better name for omega_frac?
+ float omega_frac = wi/hydro_get_physical_density(pj,cosmo)*omega_normalisation;
+
+ /* Update particle mass */
+ float current_mass = hydro_get_mass(pj);
+ float new_mass = current_mass + si->to_distribute.mass*omega_frac;
+ hydro_set_mass(pj,new_mass);
+
+ // ALEXEI: do we want to use the smoothed mass fraction?
+ /* Update total metallicity */
+ float current_metal_mass_total = pj->chemistry_data.metal_mass_fraction_total * current_mass;
+ float new_metal_mass_total = current_metal_mass_total + si->to_distribute.chemistry_data.metal_mass_fraction_total *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.metal_mass_fraction_total = new_metal_mass_total/new_mass;
+
+ /* Update mass fraction of each tracked element */
+ for (int elem = 0; elem < chemistry_element_count; elem++) {
+ float current_metal_mass = pj->chemistry_data.metal_mass_fraction[elem] * current_mass;
+ float new_metal_mass = current_metal_mass + si->to_distribute.chemistry_data.metal_mass_fraction[elem] *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass/new_mass;
+ }
+
+ /* Update iron mass fraction from SNIa */
+ float current_iron_from_SNIa_mass = pj->chemistry_data.iron_mass_fraction_from_SNIa * current_mass;
+ float new_iron_from_SNIa_mass = current_iron_from_SNIa_mass + si->to_distribute.chemistry_data.iron_mass_fraction_from_SNIa *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.iron_mass_fraction_from_SNIa = new_iron_from_SNIa_mass/new_mass;
+
+ /* Update mass fraction from SNIa (TODO: MAKE SURE IT REALLY IS A FRACTION) */
+ float current_mass_fraction_from_SNIa = pj->chemistry_data.mass_from_SNIa * current_mass;
+ float new_mass_fraction_from_SNIa = current_mass_fraction_from_SNIa + si->to_distribute.chemistry_data.mass_from_SNIa *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.mass_from_SNIa = new_mass_fraction_from_SNIa/new_mass;
+
+ /* Update metal mass fraction from SNIa */
+ float current_metal_mass_fraction_from_SNIa = pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass;
+ float new_metal_mass_fraction_from_SNIa = current_metal_mass_fraction_from_SNIa + si->to_distribute.chemistry_data.metal_mass_fraction_from_SNIa *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.metal_mass_fraction_from_SNIa = new_metal_mass_fraction_from_SNIa/new_mass;
+
+ /* Update mass fraction from SNII (TODO: MAKE SURE IT REALLY IS A FRACTION) */
+ float current_mass_fraction_from_SNII = pj->chemistry_data.mass_from_SNII * current_mass;
+ float new_mass_fraction_from_SNII = current_mass_fraction_from_SNII + si->to_distribute.chemistry_data.mass_from_SNII *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.mass_from_SNII = new_mass_fraction_from_SNII/new_mass;
+
+ /* Update metal mass fraction from SNII */
+ float current_metal_mass_fraction_from_SNII = pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass;
+ float new_metal_mass_fraction_from_SNII = current_metal_mass_fraction_from_SNII + si->to_distribute.chemistry_data.metal_mass_fraction_from_SNII *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.metal_mass_fraction_from_SNII = new_metal_mass_fraction_from_SNII/new_mass;
+
+ /* Update mass fraction from AGB (TODO: MAKE SURE IT REALLY IS A FRACTION) */
+ float current_mass_fraction_from_AGB = pj->chemistry_data.mass_from_AGB * current_mass;
+ float new_mass_fraction_from_AGB = current_mass_fraction_from_AGB + si->to_distribute.chemistry_data.mass_from_AGB *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.mass_from_AGB = new_mass_fraction_from_AGB/new_mass;
+
+ /* Update metal mass fraction from AGB */
+ float current_metal_mass_fraction_from_AGB = pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass;
+ float new_metal_mass_fraction_from_AGB = current_metal_mass_fraction_from_AGB + si->to_distribute.chemistry_data.metal_mass_fraction_from_AGB *
+ si->to_distribute.mass * omega_frac;
+ pj->chemistry_data.metal_mass_fraction_from_AGB = new_metal_mass_fraction_from_AGB/new_mass;
+}
diff --git a/src/stars/const/stars_io.h b/src/stars/const/stars_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..42fae067dbfbed7a46fbfc3c78a344b8150209a6
--- /dev/null
+++ b/src/stars/const/stars_io.h
@@ -0,0 +1,199 @@
+/*******************************************************************************
+ * 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_CONST_STARS_IO_H
+#define SWIFT_CONST_STARS_IO_H
+
+#include "io_properties.h"
+#include "stars_part.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.
+ */
+INLINE static void stars_read_particles(struct spart *sparts,
+ struct io_props *list,
+ int *num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 5;
+
+ /* 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);
+ list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_LENGTH, sparts, h);
+}
+
+/**
+ * @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.
+ */
+INLINE static void stars_write_particles(const struct spart *sparts,
+ struct io_props *list,
+ int *num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 5;
+
+ /* 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);
+ list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ sparts, h);
+}
+
+/**
+ * @brief Initialize the global properties of the stars scheme.
+ *
+ * By default, takes the values provided by the hydro.
+ *
+ * @param sp The #stars_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param p The already read-in properties of the hydro scheme.
+ */
+INLINE static void stars_props_init(struct stars_props *sp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *p) {
+
+ /* Kernel properties */
+ 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", p->h_tolerance);
+
+ /* Get derived properties */
+ sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm;
+ const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance);
+ sp->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
+ kernel_norm;
+
+ /* Maximal smoothing length */
+ sp->h_max = parser_get_opt_param_float(params, "Stars:h_max", p->h_max);
+
+ /* Number of iterations to converge h */
+ sp->max_smoothing_iterations = parser_get_opt_param_int(
+ params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
+
+ /* Time integration properties */
+ const float max_volume_change =
+ parser_get_opt_param_float(params, "Stars:max_volume_change", -1);
+ if (max_volume_change == -1)
+ sp->log_max_h_change = p->log_max_h_change;
+ else
+ sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+}
+
+/**
+ * @brief Print the global properties of the stars scheme.
+ *
+ * @param sp The #stars_props.
+ */
+INLINE static void stars_props_print(const struct stars_props *sp) {
+
+ /* Now stars */
+ message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
+ sp->eta_neighbours, sp->target_neighbours);
+
+ message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).",
+ sp->h_tolerance, sp->delta_neighbours);
+
+ message(
+ "Stars integration: Max change of volume: %.2f "
+ "(max|dlog(h)/dt|=%f).",
+ pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
+
+ if (sp->h_max != FLT_MAX)
+ message("Maximal smoothing length allowed: %.4f", sp->h_max);
+
+ message("Maximal iterations in ghost task set to %d",
+ sp->max_smoothing_iterations);
+}
+
+#if defined(HAVE_HDF5)
+INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
+ const struct stars_props *sp) {
+
+ io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
+ io_write_attribute_f(h_grpstars, "Kernel target N_ngb",
+ sp->target_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
+ io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
+ sp->h_tolerance);
+ io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
+ io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
+ sp->log_max_h_change);
+ io_write_attribute_f(h_grpstars, "Volume max change time-step",
+ pow_dimension(expf(sp->log_max_h_change)));
+ io_write_attribute_i(h_grpstars, "Max ghost iterations",
+ sp->max_smoothing_iterations);
+}
+#endif
+
+/**
+ * @brief Write a #stars_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+INLINE static void stars_props_struct_dump(const struct stars_props *p,
+ FILE *stream) {
+ restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
+ "starsprops", "stars props");
+}
+
+/**
+ * @brief Restore a stars_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+INLINE static void stars_props_struct_restore(const struct stars_props *p,
+ FILE *stream) {
+ restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
+ "stars props");
+}
+
+#endif /* SWIFT_DEFAULT_STAR_IO_H */
diff --git a/src/stars/const/stars_part.h b/src/stars/const/stars_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..32006a15b67b9e97be7aac1a86ec45d369bcc1a0
--- /dev/null
+++ b/src/stars/const/stars_part.h
@@ -0,0 +1,115 @@
+/*******************************************************************************
+ * 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>
+
+/**
+ * @brief Particle fields for the star particles.
+ *
+ * All quantities related to gravity are stored in the associate #gpart.
+ */
+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];
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff_sort[3];
+
+ /*! Particle velocity. */
+ float v[3];
+
+ /*! Star mass */
+ float mass;
+
+ /* Particle cutoff radius. */
+ float h;
+
+ /*! Particle time bin */
+ timebin_t time_bin;
+
+ struct {
+ /* Number of neighbours. */
+ float wcount;
+
+ /* Number of neighbours spatial derivative. */
+ float wcount_dh;
+
+ } density;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ /*! List of interacting particles in the density SELF and PAIR */
+ long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
+
+ /*! Number of interactions in the density SELF and PAIR */
+ int num_ngb_density;
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Contains all the constants and parameters of the stars scheme
+ */
+struct stars_props {
+
+ /*! Resolution parameter */
+ float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
+ float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
+ float delta_neighbours;
+
+ /*! Maximal smoothing length */
+ float h_max;
+
+ /*! Maximal number of iterations to converge h */
+ int max_smoothing_iterations;
+
+ /*! Maximal change of h over one time-step */
+ float log_max_h_change;
+};
+
+#endif /* SWIFT_DEFAULT_STAR_PART_H */