diff --git a/src/entropy_floor/EAGLE/entropy_floor.h b/src/entropy_floor/EAGLE/entropy_floor.h index 41d35fa0484cc1ee491a3c6293893ad5d2b5583f..5bfa3fc28c3474ee1aa87e294b417b5321ddce3c 100644 --- a/src/entropy_floor/EAGLE/entropy_floor.h +++ b/src/entropy_floor/EAGLE/entropy_floor.h @@ -136,6 +136,48 @@ static INLINE float entropy_floor( return gas_entropy_from_pressure(rho, pressure); } +static INLINE float entropy_floor_temperature( + const struct part *p, const struct cosmology *cosmo, + const struct entropy_floor_properties *props){ + + /* Physical density in internal units */ + const float rho = hydro_get_physical_density(p, cosmo); + + /* Critical density at this redshift. + * Recall that this is 0 in a non-cosmological run */ + const float rho_crit = cosmo->critical_density; + const float rho_crit_baryon = cosmo->Omega_b * rho_crit; + + /* Physical temperature */ + float temperature = 0.f; + + /* Are we in the regime of the Jeans equation of state? */ + if ((rho >= rho_crit_baryon * props->Jeans_over_density_threshold) && + (rho >= props->Jeans_density_threshold)) { + + const float jeans_slope = props->Jeans_gamma_effective - 1.f; + + const float temperature_Jeans = props->Jeans_temperature_norm * + pow(rho * props->Jeans_density_threshold_inv, jeans_slope); + + temperature = max(temperature, temperature_Jeans); + } + + /* Are we in the regime of the Cool equation of state? */ + if ((rho >= rho_crit_baryon * props->Cool_over_density_threshold) && + (rho >= props->Cool_density_threshold)) { + + const float cool_slope = props->Cool_gamma_effective - 1.f; + + const float temperature_Cool = props->Cool_temperature_norm * + pow(rho * props->Cool_density_threshold_inv, cool_slope); + + temperature = max(temperature, temperature_Cool); + } + + return temperature; +} + /** * @brief Initialise the entropy floor by reading the parameters and converting * to internal units. diff --git a/src/runner.c b/src/runner.c index dfbb9e2149246d880a9d0e303c3bc885ce40afa3..551a706734c790cfbd67f041b9805bed2da4b5f4 100644 --- a/src/runner.c +++ b/src/runner.c @@ -577,6 +577,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { const struct hydro_props *restrict hydro_props = e->hydro_properties; const struct unit_system *restrict us = e->internal_units; struct cooling_function_data *restrict cooling = e->cooling_func; + const struct entropy_floor_properties *entropy = e->entropy_floor; const double time_base = e->time_base; const integertime_t ti_current = e->ti_current; const int current_stars_count = c->stars.count; @@ -604,7 +605,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { /* Is this particle star forming? */ if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo, - hydro_props, us, cooling)) { + hydro_props, us, cooling, entropy)) { /* Time-step size for this particle */ double dt_star; diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h index b72bb38babaca51b3875147d04c46f2de95de1a7..c8cd837396e4ec635b9035d46dbb331ed94182da 100644 --- a/src/star_formation/EAGLE/star_formation.h +++ b/src/star_formation/EAGLE/star_formation.h @@ -32,6 +32,7 @@ #include "random.h" #include "stars.h" #include "units.h" +#include "entropy_floor.h" /** * @file src/star_formation/EAGLE/star_formation.h @@ -224,7 +225,8 @@ INLINE static int star_formation_is_star_forming( const struct cosmology* cosmo, const struct hydro_props* restrict hydro_props, const struct unit_system* restrict us, - const struct cooling_function_data* restrict cooling) { + const struct cooling_function_data* restrict cooling, + const struct entropy_floor_properties* restrict entropy) { /* Minimal density (converted from critical density) for star formation */ const double rho_crit_times_min_over_den = @@ -261,7 +263,7 @@ INLINE static int star_formation_is_star_forming( us, cosmo, cooling, p, xp); /* Temperature on the equation of state */ - const double temperature_eos = EOS_temperature(n_H, starform); + const double temperature_eos = entropy_floor_temperature(p,cosmo,entropy); /* Check the Scahye & Dalla Vecchia 2012 EOS-based temperature critrion */ return (temperature <