From ed6c0ba5d173e938eae240ebd9a01a13d501473d Mon Sep 17 00:00:00 2001
From: Loic Hausammann <loic_hausammann@hotmail.com>
Date: Fri, 20 Mar 2020 08:58:59 +0100
Subject: [PATCH] GEAR: discrete yields work correctly

---
 src/feedback/GEAR/initial_mass_function.c    |  23 +-
 src/feedback/GEAR/stellar_evolution.c        |  68 +++--
 src/feedback/GEAR/stellar_evolution_struct.h |  29 +-
 src/feedback/GEAR/supernovae_ii.c            | 277 +++++++++++++++----
 src/feedback/GEAR/supernovae_ii.h            |  19 +-
 5 files changed, 309 insertions(+), 107 deletions(-)

diff --git a/src/feedback/GEAR/initial_mass_function.c b/src/feedback/GEAR/initial_mass_function.c
index c0c5e2c81d..86b6ad369e 100644
--- a/src/feedback/GEAR/initial_mass_function.c
+++ b/src/feedback/GEAR/initial_mass_function.c
@@ -193,6 +193,13 @@ float initial_mass_function_get_coefficient(
 float initial_mass_function_get_integral_xi(
     const struct initial_mass_function *imf, float m1, float m2) {
 
+  /* Ensure the masses to be withing the limits */
+  m1 = min(m1, imf->mass_max);
+  m1 = max(m1, imf->mass_min);
+
+  m2 = min(m2, imf->mass_max);
+  m2 = max(m2, imf->mass_min);
+
   int k = -1;
   /* Find the correct part */
   for (int i = 0; i < imf->n_parts; i++) {
@@ -259,16 +266,14 @@ float initial_mass_function_get_imf(const struct initial_mass_function *imf,
  * @return The integral of the mass fraction.
  */
 float initial_mass_function_get_integral_imf(
-    const struct initial_mass_function *imf, const float m1, const float m2) {
+    const struct initial_mass_function *imf, float m1, float m2) {
 
-#ifdef SWIFT_DEBUG_CHECKS
-  if (m1 > imf->mass_max || m1 < imf->mass_min)
-    error("Mass 1 below or above limits expecting %g < %g < %g.", imf->mass_min,
-          m1, imf->mass_max);
-  if (m2 > imf->mass_max || m2 < imf->mass_min)
-    error("Mass 2 below or above limits expecting %g < %g < %g.", imf->mass_min,
-          m2, imf->mass_max);
-#endif
+  /* Ensure the masses to be withing the limits */
+  m1 = min(m1, imf->mass_max);
+  m1 = max(m1, imf->mass_min);
+
+  m2 = min(m2, imf->mass_max);
+  m2 = max(m2, imf->mass_min);
 
   for (int i = 0; i < imf->n_parts; i++) {
     if (m1 <= imf->mass_limits[i + 1]) {
diff --git a/src/feedback/GEAR/stellar_evolution.c b/src/feedback/GEAR/stellar_evolution.c
index c21d8f7e87..c42baae7ed 100644
--- a/src/feedback/GEAR/stellar_evolution.c
+++ b/src/feedback/GEAR/stellar_evolution.c
@@ -77,7 +77,7 @@ int stellar_evolution_compute_integer_number_supernovae(
   const float frac_sn = number_supernovae_f - number_supernovae_i;
 
   /* Get the integer number of SN */
-  return number_supernovae_i + ((rand_sn < frac_sn) ? 1 : 0);
+  return number_supernovae_i + (rand_sn < frac_sn);
 }
 
 /**
@@ -110,11 +110,11 @@ void stellar_evolution_compute_continuous_feedback_properties(
     supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia_f;
 
   /* SNII */
-  const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed(
+  const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
           &sm->snii, log_m_end_step, log_m_beg_step); 
 
   /* Sum the contributions from SNIa and SNII */
-  sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass
+  sp->feedback_data.mass_ejected = mass_frac_snii * sp->sf_data.birth_mass
     + mass_snia * phys_const->const_solar_mass;
 
   if (sp->mass <= sp->feedback_data.mass_ejected) {
@@ -132,11 +132,11 @@ void stellar_evolution_compute_continuous_feedback_properties(
 
   /* Compute the SNII yields */
   float snii_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-  supernovae_ii_get_yields(&sm->snii, log_m_end_step, log_m_beg_step,
+  supernovae_ii_get_yields_from_integral(&sm->snii, log_m_end_step, log_m_beg_step,
                            snii_yields);
 
   /* Compute the mass fraction of non processed elements */
-  const float non_processed = supernovae_ii_get_ejected_mass_fraction(
+  const float non_processed = supernovae_ii_get_ejected_mass_fraction_from_integral(
       &sm->snii, log_m_end_step, log_m_beg_step);
 
   /* Set the yields */
@@ -174,6 +174,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
  * (solMass)
  * @param m_end_step Mass of a star ending its life at the end of the step
  * (solMass)
+ * @param m_init Birth mass in solMass.
  * @param number_snia Number of SNIa produced by the stellar particle.
  * @param number_snii Number of SNII produced by the stellar particle.
  *
@@ -184,22 +185,20 @@ 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) {
 
-  /* Get the normalization to the average */
-  const float normalization =
-    number_snii == 0 ? 0. :
-    number_snii / (supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) *
-		   m_init);
+  /* 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 log_m_avg = log10(m_avg);
 
   /* Compute the mass ejected */
   /* SNIa */
   const float mass_snia =
-    (number_snia == 0) ? 0 :
-    (supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia);
+    supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia;
 
   /* SNII */
   const float mass_snii =
-      normalization * supernovae_ii_get_ejected_mass_fraction_processed(
-                          &sm->snii, log_m_end_step, log_m_beg_step);
+    supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
+      &sm->snii, log_m_avg) * m_avg * number_snii;
 
   sp->feedback_data.mass_ejected = mass_snia + mass_snii;
 
@@ -217,29 +216,33 @@ void stellar_evolution_compute_discrete_feedback_properties(
   /* Get the SNIa yields */
   const float* snia_yields = supernovae_ia_get_yields(&sm->snia);
 
-  /* Compute the SNII yields (without the normalization) */
+  /* Compute the SNII yields */
   float snii_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-  supernovae_ii_get_yields(&sm->snii, log_m_end_step, log_m_beg_step,
-                           snii_yields);
+  supernovae_ii_get_yields_from_raw(&sm->snii, log_m_avg, snii_yields);
 
   /* Compute the mass fraction of non processed elements */
   const float non_processed =
-      normalization * supernovae_ii_get_ejected_mass_fraction(
-                          &sm->snii, log_m_end_step, log_m_beg_step);
+    supernovae_ii_get_ejected_mass_fraction_from_raw(
+        &sm->snii, log_m_avg);
 
   /* Set the yields */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     /* Compute the mass fraction of metals */
     sp->feedback_data.metal_mass_ejected[i] =
-        /* Supernovae Ia yields */
-        snia_yields[i] * number_snia +
-        /* Supernovae II yields */
-        normalization * snii_yields[i] +
-        /* Gas contained in stars initial metallicity */
-        chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
+      /* Supernovae II yields */
+      snii_yields[i] +
+      /* Gas contained in stars initial metallicity */
+      chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
 
     /* Convert it to total mass */
-    sp->feedback_data.metal_mass_ejected[i] *= sp->sf_data.birth_mass;
+    sp->feedback_data.metal_mass_ejected[i] *= m_avg * number_snii;
+
+    /* Supernovae Ia yields */
+    sp->feedback_data.metal_mass_ejected[i] +=
+      snia_yields[i] * number_snia;
+
+    /* Convert everything in code units */
+    sp->feedback_data.metal_mass_ejected[i] *= phys_const->const_solar_mass;
   }
 }
 
@@ -313,8 +316,6 @@ void stellar_evolution_evolve_spart(
   /* Does this star produce a supernovae? */
   if (number_snia_f == 0 && number_snii_f == 0) return;
 
-  sp->feedback_data.number_sn = number_snia_f + number_snii_f;
-
   /* Compute the properties of the feedback (e.g. yields) */
   if (sm->discrete_yields) {
     /* Get the integer number of supernovae */
@@ -325,11 +326,22 @@ void stellar_evolution_evolve_spart(
     const int number_snii = stellar_evolution_compute_integer_number_supernovae(
       sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
 
+    /* Do we have a supernovae? */
+    if (number_snia == 0 && number_snii == 0) return;
+
+    /* Save the number of supernovae */
+    sp->feedback_data.number_sn = number_snia + number_snii;
+
+    /* Compute the yields */
     stellar_evolution_compute_discrete_feedback_properties(
         sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
         m_end_step, m_init, number_snia, number_snii);
 
   } else {
+    /* Save the number of supernovae */
+    sp->feedback_data.number_sn = number_snia_f + number_snii_f;
+
+    /* Compute the yields */
     stellar_evolution_compute_continuous_feedback_properties(
         sp, sm, phys_const, log_m_beg_step, log_m_end_step, m_beg_step,
         m_end_step, m_init, number_snia_f, number_snii_f);
diff --git a/src/feedback/GEAR/stellar_evolution_struct.h b/src/feedback/GEAR/stellar_evolution_struct.h
index cb6280c65c..c80e57711f 100644
--- a/src/feedback/GEAR/stellar_evolution_struct.h
+++ b/src/feedback/GEAR/stellar_evolution_struct.h
@@ -113,15 +113,30 @@ struct supernovae_ia {
  */
 struct supernovae_ii {
 
-  /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
-   */
-  struct interpolation_1d integrated_yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  /*! Yields not integrated */
+  struct {
+    /*! Mass fraction of metals ejected by a supernovae. */
+    struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+    
+    /*! Total mass fraction ejected. */
+    struct interpolation_1d ejected_mass_processed;
 
-  /*! Total mass fraction ejected (integrated over the IMF) */
-  struct interpolation_1d integrated_ejected_mass_processed;
+    /*! Mass fraction ejected and not processed (=> with the star metallicity). */
+    struct interpolation_1d ejected_mass;
+  } raw;
 
-  /*! Mass fraction ejected and not processed (=> with the star metallicity) */
-  struct interpolation_1d integrated_ejected_mass;
+  /*! Yields integrated */
+  struct {
+    /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
+     */
+    struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
+    
+    /*! Total mass fraction ejected (integrated over the IMF) */
+    struct interpolation_1d ejected_mass_processed;
+
+    /*! Mass fraction ejected and not processed (=> with the star metallicity) */
+    struct interpolation_1d ejected_mass;
+  } integrated;
 
   /*! Minimal mass for a SNII */
   float mass_min;
diff --git a/src/feedback/GEAR/supernovae_ii.c b/src/feedback/GEAR/supernovae_ii.c
index 6220e97553..4790a590b4 100644
--- a/src/feedback/GEAR/supernovae_ii.c
+++ b/src/feedback/GEAR/supernovae_ii.c
@@ -96,17 +96,32 @@ float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii, f
  * @param m2 The upper mass in log.
  * @param yields The elements ejected (needs to be allocated).
  */
-void supernovae_ii_get_yields(const struct supernovae_ii *snii, float log_m1,
-                              float log_m2, float *yields) {
+void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, float log_m1,
+					    float log_m2, float *yields) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    float yields_1 = interpolate_1d(&snii->integrated_yields[i], log_m1);
-    float yields_2 = interpolate_1d(&snii->integrated_yields[i], log_m2);
+    float yields_1 = interpolate_1d(&snii->integrated.yields[i], log_m1);
+    float yields_2 = interpolate_1d(&snii->integrated.yields[i], log_m2);
 
     yields[i] = yields_2 - yields_1;
   }
 };
 
+/**
+ * @brief Get the SNII yields per mass.
+ *
+ * @param snii The #supernovae_ii model.
+ * @param log_m The mass in log.
+ * @param yields The elements ejected (needs to be allocated).
+ */
+void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m,
+				       float *yields) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    yields[i] = interpolate_1d(&snii->raw.yields[i], log_m);
+  }
+};
+
 /**
  * @brief Get the ejected mass per mass unit.
  *
@@ -116,15 +131,30 @@ void supernovae_ii_get_yields(const struct supernovae_ii *snii, float log_m1,
  *
  * @return mass_ejected_processed The mass of processsed elements.
  */
-float supernovae_ii_get_ejected_mass_fraction(const struct supernovae_ii *snii,
+float supernovae_ii_get_ejected_mass_fraction_from_integral(const struct supernovae_ii *snii,
                                               float log_m1, float log_m2) {
 
-  float mass_ejected_1 = interpolate_1d(&snii->integrated_ejected_mass, log_m1);
-  float mass_ejected_2 = interpolate_1d(&snii->integrated_ejected_mass, log_m2);
+  float mass_ejected_1 = interpolate_1d(&snii->integrated.ejected_mass, log_m1);
+  float mass_ejected_2 = interpolate_1d(&snii->integrated.ejected_mass, log_m2);
 
   return mass_ejected_2 - mass_ejected_1;
 };
 
+
+/**
+ * @brief Get the ejected mass per mass unit.
+ *
+ * @param snii The #supernovae_ii model.
+ * @param log_m The mass in log.
+ *
+ * @return The mass of processsed elements.
+ */
+float supernovae_ii_get_ejected_mass_fraction_from_raw(const struct supernovae_ii *snii,
+						       float log_m) {
+
+  return interpolate_1d(&snii->raw.ejected_mass, log_m);
+};
+
 /**
  * @brief Get the ejected mass (processed) per mass.
  *
@@ -134,25 +164,48 @@ float supernovae_ii_get_ejected_mass_fraction(const struct supernovae_ii *snii,
  *
  * @return mass_ejected The mass of non processsed elements.
  */
-float supernovae_ii_get_ejected_mass_fraction_processed(
+float supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
     const struct supernovae_ii *snii, float log_m1, float log_m2) {
 
   float mass_ejected_1 =
-      interpolate_1d(&snii->integrated_ejected_mass_processed, log_m1);
+      interpolate_1d(&snii->integrated.ejected_mass_processed, log_m1);
   float mass_ejected_2 =
-      interpolate_1d(&snii->integrated_ejected_mass_processed, log_m2);
+      interpolate_1d(&snii->integrated.ejected_mass_processed, log_m2);
 
   return mass_ejected_2 - mass_ejected_1;
 };
 
+
+/**
+ * @brief Get the ejected mass (processed) per mass.
+ *
+ * @param snii The #supernovae_ii model.
+ * @param log_m The mass in log.
+ *
+ * @return mass_ejected The mass of non processsed elements.
+ */
+float supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
+    const struct supernovae_ii *snii, float log_m) {
+
+  return interpolate_1d(&snii->raw.ejected_mass_processed, log_m);
+};
+
 /**
  * @brief Read an array of SNII yields from the table.
  *
  * @param snii The #supernovae_ii model.
- * @param sm The #stellar_model.
+ * @param interp_raw Interpolation data to initialize (raw).
+ * @param interp_int Interpolation data to initialize (integrated).
+ * @param phys_const The #phys_const.
+ * @param sm * The #stellar_model.
+ * @param group_id The HDF5 group id where to read from.
+ * @param hdf5_dataset_name The dataset name to read.
+ * @param previous_count Number of element in the previous array read.
+ * @param interpolation_size Number of element to keep in the interpolation data.
  */
 void supernovae_ii_read_yields_array(
-    struct supernovae_ii *snii, struct interpolation_1d *interp,
+    struct supernovae_ii *snii, struct interpolation_1d *interp_raw,
+    struct interpolation_1d *interp_int,
     const struct phys_const *phys_const, const struct stellar_model *sm,
     hid_t group_id, const char *hdf5_dataset_name, hsize_t *previous_count,
     int interpolation_size) {
@@ -191,11 +244,16 @@ void supernovae_ii_read_yields_array(
   /* Read the dataset */
   io_read_array_dataset(group_id, hdf5_dataset_name, FLOAT, data, count);
 
+  /* Initialize the raw interpolation */
+  interpolate_1d_init(interp_raw, log10(snii->mass_min), log10(snii->mass_max),
+                      interpolation_size, log_mass_min, step_size, count, data,
+		      boundary_condition_zero);
+
   initial_mass_function_integrate(&sm->imf, data, count, log_mass_min, step_size);
   // TODO: decrease count in order to keep the same distance between points
 
-  /* Initialize the interpolation */
-  interpolate_1d_init(interp, log10(snii->mass_min), log10(snii->mass_max),
+  /* Initialize the integrated interpolation */
+  interpolate_1d_init(interp_int, log10(snii->mass_min), log10(snii->mass_max),
                       interpolation_size, log_mass_min, step_size, count, data,
 		      boundary_condition_const);
 
@@ -235,19 +293,23 @@ void supernovae_ii_read_yields(struct supernovae_ii *snii,
     const char *name = stellar_evolution_get_element_name(sm, i);
 
     /* Read the array */
-    supernovae_ii_read_yields_array(snii, &snii->integrated_yields[i],
+    supernovae_ii_read_yields_array(snii, &snii->raw.yields[i],
+				    &snii->integrated.yields[i],
                                     phys_const, sm, group_id, name,
                                     &previous_count, interpolation_size);
   }
 
   /* Read the mass ejected */
   supernovae_ii_read_yields_array(
-      snii, &snii->integrated_ejected_mass_processed, phys_const, sm, group_id,
+      snii, &snii->raw.ejected_mass_processed,
+      &snii->integrated.ejected_mass_processed,
+      phys_const, sm, group_id,
       "Ej", &previous_count, interpolation_size);
 
   /* Read the mass ejected of non processed gas */
-  supernovae_ii_read_yields_array(snii, &snii->integrated_ejected_mass,
-                                  phys_const, sm, group_id, "Ejnp",
+  supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass, 
+				  &snii->integrated.ejected_mass,
+				  phys_const, sm, group_id, "Ejnp",
                                   &previous_count, interpolation_size);
 
   /* Cleanup everything */
@@ -342,29 +404,64 @@ void supernovae_ii_dump(const struct supernovae_ii *snii, FILE *stream,
 
   /* Dump the yields. */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    if (snii->integrated_yields[i].data == NULL) {
-      continue;
+    const char *element_name = stellar_evolution_get_element_name(sm, i);
+
+    /* Integrated yields */
+    if (snii->integrated.yields[i].data != NULL) {
+      /* Generate name */
+      char name[200];
+      strcpy(name, element_name);
+      strcat(name, "_int");
+
+      /* Write the array */
+      restart_write_blocks((void *)snii->integrated.yields[i].data, sizeof(float),
+			   snii->integrated.yields[i].N, stream, name, name);
+    }
+
+    /* Raw yields */
+    if (snii->raw.yields[i].data != NULL) {
+      /* Generate name */
+      char name[200];
+      strcpy(name, element_name);
+      strcat(name, "_raw");
+
+      /* Write the array */
+      restart_write_blocks((void *)snii->raw.yields[i].data, sizeof(float),
+			   snii->raw.yields[i].N, stream, name, name);
     }
 
-    const char *name = stellar_evolution_get_element_name(sm, i);
-    restart_write_blocks((void *)snii->integrated_yields[i].data, sizeof(float),
-                         snii->integrated_yields[i].N, stream, name, name);
   }
 
-  /*! Dump the processed mass. */
-  if (snii->integrated_ejected_mass_processed.data != NULL) {
-    restart_write_blocks((void *)snii->integrated_ejected_mass_processed.data,
+  /*! Dump the processed mass (integrated). */
+  if (snii->integrated.ejected_mass_processed.data != NULL) {
+    restart_write_blocks((void *)snii->integrated.ejected_mass_processed.data,
                          sizeof(float),
-                         snii->integrated_ejected_mass_processed.N, stream,
-                         "processed_mass", "processed_mass");
+                         snii->integrated.ejected_mass_processed.N, stream,
+                         "processed_mass_int", "processed_mass_int");
   }
 
-  /*! Dump the non processed mass. */
-  if (snii->integrated_ejected_mass.data != NULL) {
-    restart_write_blocks((void *)snii->integrated_ejected_mass.data,
-                         sizeof(float), snii->integrated_ejected_mass.N, stream,
-                         "non_processed_mass", "non_processed_mass");
+  /*! Dump the processed mass (raw). */
+  if (snii->raw.ejected_mass_processed.data != NULL) {
+    restart_write_blocks((void *)snii->raw.ejected_mass_processed.data,
+                         sizeof(float),
+                         snii->raw.ejected_mass_processed.N, stream,
+                         "processed_mass_raw", "processed_mass_raw");
+  }
+
+  /*! Dump the non processed mass (integrated). */
+  if (snii->integrated.ejected_mass.data != NULL) {
+    restart_write_blocks((void *)snii->integrated.ejected_mass.data,
+                         sizeof(float), snii->integrated.ejected_mass.N, stream,
+                         "non_processed_mass_int", "non_processed_mass_int");
   }
+
+  /*! Dump the non processed mass (raw). */
+  if (snii->raw.ejected_mass.data != NULL) {
+    restart_write_blocks((void *)snii->raw.ejected_mass.data,
+                         sizeof(float), snii->raw.ejected_mass.N, stream,
+                         "non_processed_mass_raw", "non_processed_mass_raw");
+  }
+
 }
 
 /**
@@ -382,39 +479,102 @@ void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream,
 
   /* Restore the yields */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    if (snii->integrated_yields[i].data == NULL) {
-      continue;
+    const char *element_name = stellar_evolution_get_element_name(sm, i);
+
+    /* Integrated yields */
+    if (snii->integrated.yields[i].data != NULL) {
+      /* Generate name */
+      char name[200];
+      strcpy(name, element_name);
+      strcat(name, "_int");
+
+      /* Allocate the memory */
+      snii->integrated.yields[i].data =
+        (float *)malloc(sizeof(float) * snii->integrated.yields[i].N);
+      if (snii->integrated.yields[i].data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
+
+      /* Read the data */      
+      restart_read_blocks((void *)snii->integrated.yields[i].data, sizeof(float),
+			  snii->integrated.yields[i].N, stream, NULL, name);
     }
 
-    const char *name = stellar_evolution_get_element_name(sm, i);
-    snii->integrated_yields[i].data =
-        (float *)malloc(sizeof(float) * snii->integrated_yields[i].N);
+    /* Raw yields */
+    if (snii->raw.yields[i].data != NULL) {
+      /* Generate name */
+      char name[200];
+      strcpy(name, element_name);
+      strcat(name, "_raw");
+
+      /* Allocate the memory */
+      snii->raw.yields[i].data =
+        (float *)malloc(sizeof(float) * snii->raw.yields[i].N);
+      if (snii->raw.yields[i].data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
+
+      /* Read the data */      
+      restart_read_blocks((void *)snii->raw.yields[i].data, sizeof(float),
+			  snii->raw.yields[i].N, stream, NULL, name);
+    }
 
-    restart_read_blocks((void *)snii->integrated_yields[i].data, sizeof(float),
-                        snii->integrated_yields[i].N, stream, NULL, name);
   }
 
-  /* Restore the processed mass */
+  /* Restore the processed mass (integrated) */
+  if (snii->integrated.ejected_mass_processed.data != NULL) {
+    snii->integrated.ejected_mass_processed.data = (float *)malloc(
+        sizeof(float) * snii->integrated.ejected_mass_processed.N);
+      if (snii->integrated.ejected_mass_processed.data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
+
+    restart_read_blocks((void *)snii->integrated.ejected_mass_processed.data,
+                        sizeof(float),
+                        snii->integrated.ejected_mass_processed.N, stream, NULL,
+                        "processed_mass_int");
+  }
 
-  if (snii->integrated_ejected_mass_processed.data != NULL) {
-    snii->integrated_ejected_mass_processed.data = (float *)malloc(
-        sizeof(float) * snii->integrated_ejected_mass_processed.N);
+  /* Restore the processed mass (raw) */
+  if (snii->raw.ejected_mass_processed.data != NULL) {
+    snii->raw.ejected_mass_processed.data = (float *)malloc(
+        sizeof(float) * snii->raw.ejected_mass_processed.N);
+      if (snii->raw.ejected_mass_processed.data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
 
-    restart_read_blocks((void *)snii->integrated_ejected_mass_processed.data,
+    restart_read_blocks((void *)snii->raw.ejected_mass_processed.data,
                         sizeof(float),
-                        snii->integrated_ejected_mass_processed.N, stream, NULL,
-                        "processed_mass");
+                        snii->raw.ejected_mass_processed.N, stream, NULL,
+                        "processed_mass_raw");
   }
 
-  /* Restore the non processed mass */
-  if (snii->integrated_ejected_mass.data != NULL) {
-    snii->integrated_ejected_mass.data =
-        (float *)malloc(sizeof(float) * snii->integrated_ejected_mass.N);
+  /* Restore the non processed mass (integrated) */
+  if (snii->integrated.ejected_mass.data != NULL) {
+    snii->integrated.ejected_mass.data =
+        (float *)malloc(sizeof(float) * snii->integrated.ejected_mass.N);
+      if (snii->integrated.ejected_mass.data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
+
+    restart_read_blocks((void *)snii->integrated.ejected_mass.data,
+                        sizeof(float), snii->integrated.ejected_mass.N, stream,
+                        NULL, "non_processed_mass_int");
+  }
 
-    restart_read_blocks((void *)snii->integrated_ejected_mass.data,
-                        sizeof(float), snii->integrated_ejected_mass.N, stream,
-                        NULL, "non_processed_mass");
+  /* Restore the non processed mass (raw) */
+  if (snii->raw.ejected_mass.data != NULL) {
+    snii->raw.ejected_mass.data =
+        (float *)malloc(sizeof(float) * snii->raw.ejected_mass.N);
+      if (snii->raw.ejected_mass.data == NULL) {
+	error("Failed to allocate memory for the yields");
+      }
+
+    restart_read_blocks((void *)snii->raw.ejected_mass.data,
+                        sizeof(float), snii->raw.ejected_mass.N, stream,
+                        NULL, "non_processed_mass_raw");
   }
+
 }
 
 /**
@@ -425,9 +585,12 @@ void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream,
 void supernovae_ii_clean(struct supernovae_ii *snii) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    interpolate_1d_free(&snii->integrated_yields[i]);
+    interpolate_1d_free(&snii->integrated.yields[i]);
+    interpolate_1d_free(&snii->raw.yields[i]);
   }
 
-  interpolate_1d_free(&snii->integrated_ejected_mass_processed);
-  interpolate_1d_free(&snii->integrated_ejected_mass);
+  interpolate_1d_free(&snii->integrated.ejected_mass_processed);
+  interpolate_1d_free(&snii->raw.ejected_mass_processed);
+  interpolate_1d_free(&snii->integrated.ejected_mass);
+  interpolate_1d_free(&snii->raw.ejected_mass);
 }
diff --git a/src/feedback/GEAR/supernovae_ii.h b/src/feedback/GEAR/supernovae_ii.h
index 558cb56b36..3f525df502 100644
--- a/src/feedback/GEAR/supernovae_ii.h
+++ b/src/feedback/GEAR/supernovae_ii.h
@@ -33,15 +33,22 @@ int supernovae_ii_can_explode(const struct supernovae_ii *snii, float m_low,
                               float m_high);
 float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii, float m1,
 					     float m2);
-void supernovae_ii_get_yields(const struct supernovae_ii *snii, float log_m1,
-                              float log_m2, float *yields);
-float supernovae_ii_get_ejected_mass_fraction(const struct supernovae_ii *snii,
-                                              float log_m1, float log_m2);
+void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, float log_m1,
+					    float log_m2, float *yields);
+void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m,
+				       float *yields);
+float supernovae_ii_get_ejected_mass_fraction_from_integral(const struct supernovae_ii *snii,
+							    float log_m1, float log_m2);
+float supernovae_ii_get_ejected_mass_fraction_from_raw(const struct supernovae_ii *snii,
+						       float log_m);
 
-float supernovae_ii_get_ejected_mass_fraction_processed(
+float supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
     const struct supernovae_ii *snii, float log_m1, float log_m2);
+float supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
+    const struct supernovae_ii *snii, float log_m);
 void supernovae_ii_read_yields_array(
-    struct supernovae_ii *snii, struct interpolation_1d *interp,
+    struct supernovae_ii *snii, struct interpolation_1d *interp_raw,
+    struct interpolation_1d *interp_int,
     const struct phys_const *phys_const, const struct stellar_model *sm,
     hid_t group_id, const char *hdf5_dataset_name, hsize_t *previous_count,
     int interpolation_size);
-- 
GitLab