Update GEAR feedback

src/feedback/GEAR/feedback.c

@@ -225,7 +225,13 @@ void feedback_evolve_spart(struct spart* restrict sp,
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
+
+  if (star_age_beg_step < -1e-6) {
+    error("Negative age for a star");
+  }
 #endif
+  const double star_age_beg_step_safe =
+      star_age_beg_step < 0 ? 0 : star_age_beg_step;
 
   /* Reset the feedback */
   feedback_reset_feedback(sp, feedback_props);
@@ -238,7 +244,8 @@ void feedback_evolve_spart(struct spart* restrict sp,
 
   /* Compute the stellar evolution */
   stellar_evolution_evolve_spart(sp, &feedback_props->stellar_model, cosmo, us,
-                                 phys_const, ti_begin, star_age_beg_step, dt);
+                                 phys_const, ti_begin, star_age_beg_step_safe,
+                                 dt);
 
   /* Transform the number of SN to the energy */
   sp->feedback_data.energy_ejected =

src/feedback/GEAR/stellar_evolution.c

@@ -187,10 +187,12 @@ void stellar_evolution_compute_discrete_feedback_properties(
     const float log_m_end_step, const float m_beg_step, const float m_end_step,
     const float m_init, const int number_snia, const int number_snii) {
 
+  /* Limit the mass within the imf limits */
+  const float m_beg_lim = min(m_beg_step, sm->imf.mass_max);
+  const float m_end_lim = max(m_end_step, sm->imf.mass_min);
+
   /* Compute the average mass */
-  const float m_avg =
-      initial_mass_function_get_integral_imf(&sm->imf, m_end_step, m_beg_step) /
-      initial_mass_function_get_integral_xi(&sm->imf, m_end_step, m_beg_step);
+  const float m_avg = 0.5 * (m_beg_lim + m_end_lim);
   const float log_m_avg = log10(m_avg);
 
   /* Compute the mass ejected */