diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py index f83607e03a162df5ecf83b675b2e0586bb71268a..73e4878e8e00a35fe19c359652be0d57153dea62 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py @@ -123,6 +123,7 @@ stars_pos = f["/PartType4/Coordinates"][:, :] stars_BirthDensity = f["/PartType4/BirthDensity"][:] stars_BirthTime = f["/PartType4/BirthTime"][:] stars_XH = f["/PartType4/ElementAbundance"][:,0] +stars_BirthTemperature = f["/PartType4/BirthTemperature"][:] # Centre the box gas_pos[:, 0] -= centre[0] @@ -207,9 +208,17 @@ figure() subplot(111, xscale="linear", yscale="linear") hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5]) xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0) -ylabel("${\\rm Probability}$", labelpad=-7) +ylabel("${\\rm Probability}$", labelpad=3) savefig("BirthDensity.png", dpi=200) +# Histogram of the birth temperature +figure() +subplot(111, xscale="linear", yscale="linear") +hist(np.log10(stars_BirthTemperature),density=True,bins=20,range=[3.5,5.0]) +xlabel("${\\rm Stellar~birth~temperature}~[{\\rm K}]$", labelpad=0) +ylabel("${\\rm Probability}$", labelpad=3) +savefig("BirthTemperature.png", dpi=200) + # Plot of the specific star formation rate in the galaxy rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100) rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100) diff --git a/src/runner.c b/src/runner.c index ff2564d8b4b2fb4d6699b295c805cbfd2e6cbd41..3ff1c1274c09936c7b77e939a5c0b6b88750858d 100644 --- a/src/runner.c +++ b/src/runner.c @@ -735,7 +735,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { /* Copy the properties of the gas particle to the star particle */ star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo, - with_cosmology); + with_cosmology, phys_const, + hydro_props, us, cooling); /* Update the Star formation history */ star_formation_logger_log_new_spart(sp, &c->stars.sfh); diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h index 40161aded43872dbd66c8007892f62e8eafadb7a..4d1013c0f6dd13552cb533aa942f47dd588a469c 100644 --- a/src/star_formation/EAGLE/star_formation.h +++ b/src/star_formation/EAGLE/star_formation.h @@ -210,7 +210,7 @@ INLINE static int star_formation_is_star_forming( const struct hydro_props* restrict hydro_props, const struct unit_system* restrict us, const struct cooling_function_data* restrict cooling, - const struct entropy_floor_properties* restrict entropy_floor) { + const struct entropy_floor_properties* restrict entropy_floor_props) { /* Minimal density (converted from critical density) for star formation */ const double rho_crit_times_min_over_den = @@ -242,17 +242,15 @@ INLINE static int star_formation_is_star_forming( /* Check if it exceeded the minimum density */ if (n_H < density_threshold) return 0; - /* Calculate the temperature */ - const double temperature = cooling_get_temperature(phys_const, hydro_props, - us, cosmo, cooling, p, xp); + /* Calculate the entropy of the particle */ + const double entropy = hydro_get_physical_entropy(p, xp, cosmo); - /* Temperature on the equation of state */ - const double temperature_eos = - entropy_floor_temperature(p, cosmo, entropy_floor); + /* Calculate the entropy EOS of the particle */ + const double entropy_eos = entropy_floor(p, cosmo, entropy_floor_props); /* Check the Scahye & Dalla Vecchia 2012 EOS-based temperature critrion */ - return (temperature < - temperature_eos * starform->ten_to_temperature_margin_threshold_dex); + return (entropy < + entropy_eos * starform->ten_to_temperature_margin_threshold_dex); } /** @@ -382,7 +380,11 @@ INLINE static void star_formation_update_part_not_SFR( INLINE static void star_formation_copy_properties( const struct part* p, const struct xpart* xp, struct spart* sp, const struct engine* e, const struct star_formation* starform, - const struct cosmology* cosmo, const int with_cosmology) { + const struct cosmology* cosmo, const int with_cosmology, + const struct phys_const* phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cooling_function_data* restrict cooling) { /* Store the current mass */ sp->mass = hydro_get_mass(p); @@ -406,6 +408,10 @@ INLINE static void star_formation_copy_properties( /* Store the birth density in the star particle */ sp->birth_density = hydro_get_physical_density(p, cosmo); + /* Store the birth temperature in the star particle */ + sp->birth_temperature = cooling_get_temperature(phys_const, hydro_props, us, + cosmo, cooling, p, xp); + /* Flag that this particle has not done feedback yet */ sp->f_E = -1.f; } diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index 05bb584a5330002e610df9c1866c672f5ae1977f..c479feb5c66328f9fab8bf62593ca66b6658b79e 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -129,7 +129,11 @@ INLINE static void star_formation_update_part_not_SFR( INLINE static void star_formation_copy_properties( const struct part* p, const struct xpart* xp, struct spart* sp, const struct engine* e, const struct star_formation* starform, - const struct cosmology* cosmo, const int with_cosmology) {} + const struct cosmology* cosmo, const int with_cosmology, + const struct phys_const* phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cooling_function_data* restrict cooling) {} /** * @brief initialization of the star formation law diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h index 25c0bcee6a8d2aa2f4be5479556eebc50027aa72..dfe645718d689841f89cf592194d435af299a642 100644 --- a/src/star_formation/none/star_formation.h +++ b/src/star_formation/none/star_formation.h @@ -132,7 +132,11 @@ INLINE static void star_formation_update_part_not_SFR( INLINE static void star_formation_copy_properties( const struct part* p, const struct xpart* xp, struct spart* sp, const struct engine* e, const struct star_formation* starform, - const struct cosmology* cosmo, const int with_cosmology) {} + const struct cosmology* cosmo, const int with_cosmology, + const struct phys_const* phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cooling_function_data* restrict cooling) {} /** * @brief initialization of the star formation law diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h index ac663cc4b3ba6c96ba3230511f14d3235b347f3d..0c03bee3007066c7c51c7ce0fb3d88d37a1b2ae3 100644 --- a/src/stars/EAGLE/stars_io.h +++ b/src/stars/EAGLE/stars_io.h @@ -64,7 +64,7 @@ INLINE static void stars_write_particles(const struct spart *sparts, int *num_fields) { /* Say how much we want to write */ - *num_fields = 9; + *num_fields = 10; /* List what we want to write */ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, @@ -85,6 +85,9 @@ INLINE static void stars_write_particles(const struct spart *sparts, birth_time); list[8] = io_make_output_field("FeedbackEnergyFraction", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts, f_E); + list[9] = + io_make_output_field("BirthTemperature", FLOAT, 1, UNIT_CONV_TEMPERATURE, + sparts, birth_temperature); } /** diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h index f8ba330b2f09ac0966a35dfd6e2468b7193c38ed..bed54bb175026ae930668fb102dcfef9c9af76d7 100644 --- a/src/stars/EAGLE/stars_part.h +++ b/src/stars/EAGLE/stars_part.h @@ -85,6 +85,9 @@ struct spart { /*! Birth density */ float birth_density; + /*! Birth temperature */ + float birth_temperature; + /*! Feedback energy fraction */ float f_E;