Prevent rounding errors when computing stellar ages in the EAGLE feedback model.

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;