diff --git a/src/Makefile.am b/src/Makefile.am
index 032a6dfaf23329bc56fdb55baffab968c4501be5..5a356ba6707eee0664b9cb7817c8f2195d5512b2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -50,7 +50,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h \
- star_formation_struct.h star_formation.h velociraptor_struct.h velociraptor_io.h
+ star_formation_struct.h star_formation.h velociraptor_struct.h velociraptor_io.h \
+ random.h
# source files for EAGLE cooling
EAGLE_COOLING_SOURCES =
diff --git a/src/random.h b/src/random.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a8f45e69bc5c1b1d771a3fd6b5bab494f034008
--- /dev/null
+++ b/src/random.h
@@ -0,0 +1,64 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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_RANDOM_H
+#define SWIFT_RANDOM_H
+
+/* COde configuration */
+#include "../config.h"
+
+/* Standard header */
+#include <stdlib.h>
+
+/**
+ * @brief The categories of random number generated.
+ */
+enum random_number_type {
+ random_number_star_formation = 0,
+ random_number_stellar_feedback,
+ random_number_stellar_enrichment,
+ random_number_BH_feedback
+};
+
+/**
+ * @brief Returns a pseudo-random number in the range [0, 1[.
+ *
+ * We generate numbers that are always reproducible for a given particle ID and
+ * simulation time. If more than one number per time-step per particle is
+ * needed, additional randomness can be obtained by using the type argument.
+ *
+ * @param id The ID of the particle for which to generate a number.
+ * @param ti_current The time (on the time-line) for which to generate a number.
+ * @prarm type The #random_number_type to generate.
+ */
+INLINE static double random_unit_interval(const long long int id,
+ const integertime_t ti_current,
+ const enum random_number_type type) {
+
+ /* Range used for the seeds. Best if prime */
+ static const long long seed_range = 131071LL;
+ static const double RAND_MAX_inv = 1. / ((double)RAND_MAX);
+
+ /* Calculate the seed */
+ unsigned int seed = (id + ti_current + (long long)type) % seed_range;
+
+ /* Generate a random number between 0 and 1. */
+ return rand_r(&seed) / RAND_MAX_inv;
+}
+
+#endif /* SWIFT_RANDOM_H */
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 88c09186d26da5e4dbb12f2a8e907928948bfff0..ceb440250ba22f49933527eb3da04c7414285be9 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -29,6 +29,7 @@
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
+#include "random.h"
#include "stars.h"
#include "units.h"
@@ -296,23 +297,25 @@ INLINE static int star_formation_convert_to_star(
const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
const double n_H = physical_density * X_H / phys_const->const_proton_mass;
+ /* Are we above the threshold for automatic star formation? */
+ if (physical_density >
+ starform->max_gas_density * phys_const->const_proton_mass) {
+ return 1;
+ }
+
/* Pressure on the EOS for this particle */
const double pressure = EOS_pressure(n_H, starform);
+ /* Calculate the specific star formation rate */
double SFRpergasmass;
if (hydro_get_physical_density(p, cosmo) <
starform->KS_high_den_thresh * phys_const->const_proton_mass) {
- /* Calculate the star formation rate */
+
SFRpergasmass =
starform->SF_normalization * pow(pressure, starform->SF_power_law);
- } else if (hydro_get_physical_density(p, cosmo) >
- starform->max_gas_density * phys_const->const_proton_mass) {
- /* We give the star formation tracers values of the mass and -1 for sSFR
- * to be able to trace them, we don't want random variables there*/
- xp->sf_data.SFR = p->mass;
- xp->sf_data.sSFR = -1.f;
- return 1.f;
+
} else {
+
SFRpergasmass = starform->SF_high_den_normalization *
pow(pressure, starform->SF_high_den_power_law);
}
@@ -320,18 +323,16 @@ INLINE static int star_formation_convert_to_star(
/* Store the SFR */
xp->sf_data.SFR = SFRpergasmass * p->mass;
xp->sf_data.sSFR = SFRpergasmass;
+
/* Calculate the propability of forming a star */
const double prop = SFRpergasmass * dt_star;
- /* Calculate the seed */
- unsigned int seed = (p->id + e->ti_current) % 8191;
-
- /* Generate a random number between 0 and 1. */
- const double randomnumber =
- rand_r(&seed); // MATTHIEU: * starform->inv_RAND_MAX;
+ /* Get a random number between 0 and 1 for star formation */
+ const double random_number = random_unit_interval(
+ p->id, e->ti_current, random_number_star_formation);
- /* Calculate if we form a star */
- return (prop > randomnumber);
+ /* Have we been lucky and need to form a star? */
+ return (prop > random_number);
}
/* Check if it is the first time steps after star formation */