diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c index 87f8e5cc48614b44660b87529a31d64aaa354ca3..4cd6a45a1056a67168ac91b0bb961dbbc2528b2e 100644 --- a/src/feedback/EAGLE/feedback.c +++ b/src/feedback/EAGLE/feedback.c @@ -114,7 +114,7 @@ double eagle_feedback_energy_fraction(const struct spart* sp, /* Calculate f_E */ const double Z_term = pow(max(Z_smooth, 1e-6) / Z_0, n_Z); - const double n_term = pow(max(n_birth, 1e-6) / n_0, -n_n); + const double n_term = pow(n_birth / n_0, -n_n); const double denonimator = 1. + Z_term * n_term; return f_E_min + (f_E_max - f_E_min) / denonimator; @@ -181,6 +181,11 @@ INLINE static void compute_SNII_feedback( delta_u = f_E * E_SNe * N_SNe / ngb_gas_mass; } +#ifdef SWIFT_DEBUG_CHECKS + if (f_E < feedback_props->f_E_min || f_E > feedback_props->f_E_max) + error("f_E is not in the valid range! f_E=%f sp->id=%lld", f_E, sp->id); +#endif + /* Store all of this in the star for delivery onto the gas */ sp->f_E = f_E; sp->feedback_data.to_distribute.SNII_heating_probability = prob; @@ -309,12 +314,17 @@ INLINE static void determine_bin_yield_SNII( INLINE static void evolve_SNIa(const float log10_min_mass, const float log10_max_mass, const struct feedback_props* props, - struct spart* sp, float star_age_Gyr, - float dt_Gyr) { + struct spart* sp, double star_age_Gyr, + double dt_Gyr) { /* Check if we're outside the mass range for SNIa */ if (log10_min_mass >= props->log10_SNIa_max_mass_msun) return; +#ifdef SWIFT_DEBUG_CHECKS + if (dt_Gyr < 0.) error("Negative time-step length!"); + if (star_age_Gyr < 0.) error("Negative age!"); +#endif + /* If the max mass is outside the mass range update it to be the maximum * and use updated values for the star's age and timestep in this function */ if (log10_max_mass > props->log10_SNIa_max_mass_msun) { @@ -323,15 +333,10 @@ INLINE static void evolve_SNIa(const float log10_min_mass, const float max_mass = exp10f(props->log10_SNIa_max_mass_msun); const float lifetime_Gyr = lifetime_in_Gyr(max_mass, Z, props); - dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr; + dt_Gyr = max(star_age_Gyr + dt_Gyr - lifetime_Gyr, 0.); star_age_Gyr = lifetime_Gyr; } -#ifdef SWIFT_DEBUG_CHECKS - if (dt_Gyr < 0.) error("Negative time-step length!"); - if (star_age_Gyr < 0.) error("Negative age!"); -#endif - /* Compute the number of SNIa */ const float num_SNIa = eagle_feedback_number_of_SNIa( sp, star_age_Gyr, star_age_Gyr + dt_Gyr, props); @@ -704,9 +709,9 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, /* Convert dt and stellar age from internal units to Gyr. */ const double Gyr_in_cgs = 1e9 * 365. * 24. * 3600.; const double time_to_cgs = units_cgs_conversion_factor(us, UNIT_CONV_TIME); - const float conversion_factor = time_to_cgs / Gyr_in_cgs; - const float dt_Gyr = dt * conversion_factor; - const float star_age_Gyr = age * conversion_factor; + const double conversion_factor = time_to_cgs / Gyr_in_cgs; + const double dt_Gyr = dt * conversion_factor; + const double star_age_Gyr = age * conversion_factor; /* Get the metallicity */ const float Z = sp->chemistry_data.metal_mass_fraction_total;