Yields correctly implemented

examples/supernovae_yields/plot_yields.py

@@ -63,6 +63,7 @@ if __name__ == "__main__":
             filename, m1, mass)
 
         yields_snia = libstellar.supernovaeIaYields(filename)
+        mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
 
         elements = libstellar.getElements(filename)
 
@@ -76,7 +77,8 @@ if __name__ == "__main__":
             p = plt.loglog(mass, yields_snii[:, i], label=elements[i])
             plt.loglog(m_wd, yields_snia[i], "o", color=p[0].get_color())
 
-        plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+        p = plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+        plt.loglog(m_wd, mass_ejected_snia, "o", color=p[0].get_color())
 
         plt.loglog(mass, mass_ejected_np_snii,
                    label="Non Processed Mass Ejected")

pyswiftsim/__init__.py

@@ -168,7 +168,8 @@ class ParameterFile:
             "SelfShieldingMethod": 0,
             "ConvergenceLimit": 1e-2,
             "MaxSteps": 20000,
-            "Omega": 0.8
+            "Omega": 0.8,
+            "ThermalTime_Myr": 5,
         },
 
         "GearChemistry": {
@@ -177,7 +178,6 @@ class ParameterFile:
 
         "GEARFeedback": {
             "SupernovaeEnergy_erg": 1e51,
-            "ThermalTime_Myr": 5,
             "YieldsTable": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
         },
     }

scripts/pyswift_plot_yields

+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+from time import strftime, gmtime
+from argparse import ArgumentParser
+import os
+plt.rc('text', usetex=True)
+
+
+N = 1000
+m_min = 0.08  # in solar mass
+m_max = 50
+m_wd = 1.4  # white dwarf mass
+solMass_in_cgs = 1.989e33
+
+
+def parseArguments():
+    parser = ArgumentParser(description="Plot the supernovae rate")
+
+    parser.add_argument(
+        "--download_table", dest="table",
+        action="store_true",
+        help="Download the default chemistry table")
+
+    parser.add_argument(
+        "--params", dest="params",
+        type=str, default=None,
+        help="SWIFT parameters file")
+
+    parser.add_argument(
+        "--N", dest="N",
+        type=int, default=1000,
+        help="Number of points")
+
+    parser.add_argument(
+        "--m_min", dest="m_min",
+        type=float, default=0.08,
+        help="Minimal mass to plot in solar mass")
+
+    parser.add_argument(
+        "--m_max", dest="m_max",
+        type=float, default=50,
+        help="Maximal mass to plot in solar mass")
+
+    parser.add_argument(
+        "--disable-snia", dest="disable_snia",
+        action="store_true",
+        help="Disable the SNIa")
+
+    parser.add_argument(
+        "--disable-snii", dest="disable_snii",
+        action="store_true",
+        help="Disable the SNII")
+
+    args = parser.parse_args()
+
+    return args
+
+
+if __name__ == "__main__":
+    opt = parseArguments()
+
+    if opt.table:
+        pyswiftsim.downloadYieldsTable()
+
+    if opt.params is not None and not os.path.isfile(opt.params):
+        raise Exception(
+            "The parameter file '{}' does not exist".format(opt.params))
+
+    with pyswiftsim.ParameterFile(existing=opt.params) as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
+
+        # Get masses
+        mass = np.logspace(np.log10(opt.m_min), np.log10(opt.m_max), N,
+                           dtype=np.float32)
+
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = float(units["UnitMass_in_cgs"])
+        unit_vel = float(units["UnitVelocity_in_cgs"])
+        unit_length = float(units["UnitLength_in_cgs"])
+        unit_time = unit_length / unit_vel
+
+        mass *= solMass_in_cgs / unit_mass
+
+        m1 = np.ones(N, dtype=np.float32) * 1e-20  # avoid log(0)
+
+        print("Computing supernovae yields...")
+
+        yields_snii = libstellar.supernovaeIIYields(filename, m1, mass)
+        mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+            filename, m1, mass)
+        mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+            filename, m1, mass)
+
+        yields_snia = libstellar.supernovaeIaYields(filename)
+        mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
+
+        elements = libstellar.getElements(filename)
+
+        # units
+        mass /= solMass_in_cgs / unit_mass
+
+        plt.figure()
+        date = strftime("%d %b %Y", gmtime())
+        plt.title("Date: {}".format(date))
+        for i in range(len(elements)):
+            p = plt.loglog(mass, yields_snii[:, i], label=elements[i])
+            plt.loglog(m_wd, yields_snia[i], "o", color=p[0].get_color())
+
+        p = plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+        plt.loglog(m_wd, mass_ejected_snia, "o", color=p[0].get_color())
+
+        plt.loglog(mass, mass_ejected_np_snii,
+                   label="Non Processed Mass Ejected")
+
+        plt.xlabel("$Mass [M_\odot]$")
+        plt.ylabel("Mass fraction")
+        plt.legend()
+
+        plt.show()

src/cooling_wrapper.c

@@ -163,7 +163,7 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
 
       /* compute cooling rate */
       cooling_cool_part(&pconst, &us, &cosmo, &hydro_props,
-			&floor_props, &cooling, &p, &xp, dt, dt);
+			&floor_props, &cooling, &p, &xp, /* Time */0, dt, dt);
       float *tmp = PyArray_GETPTR1(rate, i);
       *tmp = hydro_get_physical_internal_energy_dt(&p, &cosmo);
 

src/libstellar.c

@@ -130,9 +130,9 @@ static PyMethodDef libstellar_methods[] = {
    "filename: str\n"
    "\t Parameter file\n"
    "time1: array\n"
-   "\t Beginning of the time interval for the SNIa in internal units.\n"
+   "\t Beginning of the time interval for the SNII in internal units.\n"
    "time2: array\n"
-   "\t End of the time interval for the SNIa in internal units.\n"
+   "\t End of the time interval for the SNII in internal units.\n"
    "metallicity: array\n"
    "\t Metallicity of the stars.\n"
    "\n"
@@ -150,7 +150,127 @@ static PyMethodDef libstellar_methods[] = {
    ">>> mass = libstellar.supernovaeIIRate(filename, time1, time2, Z)\n"
   },
 
+  {"supernovaeIIYields", pysupernovae_ii_yields, METH_VARARGS,
+   "Compute the supernovae II yields from the masses.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "mass1: array\n"
+   "\t Beginning of the mass interval for the SNII in internal units.\n"
+   "mass2: array\n"
+   "\t End of the time interval for the SNII in internal units.\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "yields: array\n"
+   "\t The yields (mass fraction) of SNII.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass1 = np.array([1.1, 1.3])\n"
+   ">>> mass2 = np.array([1.2, 1.4])\n"
+   ">>> mass = libstellar.supernovaeIIYields('params.yml', mass1, mass2)\n"
+  },
+
+  {"supernovaeIaYields", pysupernovae_ia_yields, METH_VARARGS,
+   "Compute the supernovae Ia yields from the masses.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "yields: array\n"
+   "\t The yields (mass fraction) of SNIa.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass = libstellar.supernovaeIaYields('params.yml')\n"
+  },
+
+  {"getElements", pystellar_get_elements, METH_VARARGS,
+   "Read the element names..\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "element_names: array\n"
+   "\t The name of all the elements.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass = libstellar.getElements('params.yml')\n"
+   "[Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu, Metals]"
+  },
 
+  {"supernovaeIIMassEjected", pysupernovae_ii_mass_ejected, METH_VARARGS,
+   "Compute the mass ejected by supernovae II from the masses.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "mass1: array\n"
+   "\t Beginning of the mass interval for the SNII in internal units.\n"
+   "mass2: array\n"
+   "\t End of the time interval for the SNII in internal units.\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "mass_ejected: array\n"
+   "\t The mass ejected by the supernovae SNII.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass1 = np.array([1.1, 1.3])\n"
+   ">>> mass2 = np.array([1.2, 1.4])\n"
+   ">>> mass = libstellar.supernovaeIIMassEjected('params.yml', mass1, mass2)\n"
+  },
+
+  {"supernovaeIIMassEjectedNonProcessed", pysupernovae_ii_mass_ejected_non_processed, METH_VARARGS,
+   "Compute the (non processed) mass ejected by supernovae II from the masses.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "mass1: array\n"
+   "\t Beginning of the mass interval for the SNII in internal units.\n"
+   "mass2: array\n"
+   "\t End of the time interval for the SNII in internal units.\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "mass_ejected: array\n"
+   "\t The non processed mass ejected by the supernovae SNII.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass1 = np.array([1.1, 1.3])\n"
+   ">>> mass2 = np.array([1.2, 1.4])\n"
+   ">>> mass = libstellar.supernovaeIIMassEjectedNonProcessed('params.yml', mass1, mass2)\n"
+  },
+
+  {"supernovaeIaMassEjected", pysupernovae_ia_mass_ejected, METH_VARARGS,
+   "Compute the mass ejected by supernovae Ia.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "mass_ejected: float\n"
+   "\t The mass ejected by the supernovae SNIa.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> mass = libstellar.supernovaeIaMassEjected('params.yml')\n"
+  },
+  
   {NULL, NULL, 0, NULL}        /* Sentinel */
 };      
       

src/stellar_evolution_wrapper.c

@@ -64,7 +64,7 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
     float *imf_i = (float*) PyArray_GETPTR1(imf, i);
 
 
-    *imf_i = stellar_evolution_get_imf(&fp.stellar_model.imf, m);
+    *imf_i = initial_mass_function_get_imf(&fp.stellar_model.imf, m);
   }
 
   return (PyObject *) imf;
@@ -125,7 +125,7 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
     float *tau = (float*) PyArray_GETPTR1(lifetime, i);
 
     
-    *tau = stellar_evolution_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
+    *tau = lifetime_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
     *tau = pow(10, *tau);
   }
 
@@ -187,7 +187,7 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
     float *m = (float*) PyArray_GETPTR1(mass, i);
 
     
-    *m = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
+    *m = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
     *m = pow(10, *m);
   }
 
@@ -208,7 +208,7 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
  *
  * @param self calling object
  * @param args arguments
- * @return cooling rate
+ * @return Rate of supernovae Ia
  */
 PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
   import_array();
@@ -259,13 +259,13 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
 
     
     
-    float m1 = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
+    float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
     m1 = pow(10, m1);
-    float m2 = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
+    float m2 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
     m2 = pow(10, m2);
 
     /* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
-    *r = stellar_evolution_get_number_supernovae_ia(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
+    *r = supernovae_ia_get_number(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
   }
 
   return (PyObject *) rate;
@@ -283,7 +283,7 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
  *
  * @param self calling object
  * @param args arguments
- * @return cooling rate
+ * @return Rate of supernovae II
  */
 PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
   import_array();
@@ -334,15 +334,352 @@ PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
 
     
     
-    float m1 = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
+    float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
     m1 = pow(10, m1);
-    float m2 = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
+    float m2 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
     m2 = pow(10, m2);
 
     /* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
-    *r = stellar_evolution_get_number_supernovae_ii(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
+    *r = supernovae_ii_get_number(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
   }
 
   return (PyObject *) rate;
   
 }
+
+/**
+ * @brief Compute the supernovae II yields.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *  - a numpy array containing the lower limit for the mass.
+ *  - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction of ejected metals
+ */
+PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  PyArrayObject *mass1;
+  PyArrayObject *mass2;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "sOO", &filename, &mass1, &mass2))
+    return NULL;
+
+  /* check numpy array */
+  if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+    return NULL;
+
+  if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+    error("The inputs should be 1-D arrays");
+
+  if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+    error("Mass1 and Mass2 should have the same dimension");
+
+
+  size_t N = PyArray_DIM(mass1, 0);
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* return object */
+  const int ndims = 2;
+  npy_intp dims[2] = {N, CHEMISTRY_ELEMENT_COUNT};
+  PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+  for(size_t i = 0; i < N; i++) {
+    float yields[CHEMISTRY_ELEMENT_COUNT] = {0.};
+    const float m1 = *(float*) PyArray_GETPTR1(mass1, i);
+    const float log_m1 = log10(m1);
+    const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+    const float log_m2 = log10(m2);
+
+    /* Compute the yields */
+    supernovae_ii_get_yields(&fp.stellar_model.snii, log_m1, log_m2, yields);
+
+    for(int j = 0; j < CHEMISTRY_ELEMENT_COUNT; j++) {
+      float *y = (float*) PyArray_GETPTR2(yields_out, i, j);
+      *y = yields[j];
+    }
+  }
+
+  return (PyObject *) yields_out;
+  
+}
+
+/**
+ * @brief Compute the mass ejected by a supernovae II.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *  - a numpy array containing the lower limit for the mass.
+ *  - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction of ejected mass.
+ */
+PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  PyArrayObject *mass1;
+  PyArrayObject *mass2;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "sOO", &filename, &mass1, &mass2))
+    return NULL;
+
+  /* check numpy array */
+  if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+    return NULL;
+
+  if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+    error("The inputs should be 1-D arrays");
+
+  if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+    error("Mass1 and Mass2 should have the same dimension");
+
+
+  size_t N = PyArray_DIM(mass1, 0);
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* return object */
+  PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(
+      PyArray_NDIM(mass1), PyArray_DIMS(mass1), NPY_FLOAT);
+
+  for(size_t i = 0; i < N; i++) {
+    const float m1 = *(float*) PyArray_GETPTR1(mass1, i);
+    const float log_m1 = log10(m1);
+    const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+    const float log_m2 = log10(m2);
+    float *y = (float*) PyArray_GETPTR1(mass_ejected, i);
+
+    *y = supernovae_ii_get_ejected_mass_processed(&fp.stellar_model.snii, log_m1, log_m2);
+
+  }
+
+  return (PyObject *) mass_ejected;
+  
+}
+
+/**
+ * @brief Compute the non processed mass ejected by a supernovae II.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *  - a numpy array containing the lower limit for the mass.
+ *  - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass ratio of (non processed) mass ejected
+ */
+PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  PyArrayObject *mass1;
+  PyArrayObject *mass2;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "sOO", &filename, &mass1, &mass2))
+    return NULL;
+
+  /* check numpy array */
+  if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+    return NULL;
+
+  if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+    error("The inputs should be 1-D arrays");
+
+  if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+    error("Mass1 and Mass2 should have the same dimension");
+
+
+  size_t N = PyArray_DIM(mass1, 0);
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* return object */
+  PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(
+      PyArray_NDIM(mass1), PyArray_DIMS(mass1), NPY_FLOAT);
+
+  for(size_t i = 0; i < N; i++) {
+    const float m1 = *(float*) PyArray_GETPTR1(mass1, i);
+    const float log_m1 = log10(m1);
+    const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+    const float log_m2 = log10(m2);
+    float *y = (float*) PyArray_GETPTR1(mass_ejected, i);
+
+    *y = supernovae_ii_get_ejected_mass(&fp.stellar_model.snii, log_m1, log_m2);
+    
+  }
+
+  return (PyObject *) mass_ejected;
+  
+}
+
+/**
+ * @brief Compute the supernovae Ia yields.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return Mass fraction of the supernovae Ia yields.
+ */
+PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "s", &filename))
+    return NULL;
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* return object */
+  const int ndims = 1;
+  npy_intp dims[1] = {CHEMISTRY_ELEMENT_COUNT};
+  PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+  const float *yields = supernovae_ia_get_yields(&fp.stellar_model.snia);
+
+  for(size_t i = 0; i < CHEMISTRY_ELEMENT_COUNT; i++) {
+    float *y = (float*) PyArray_GETPTR1(yields_out, i);
+
+    *y = yields[i] / fp.stellar_model.snia.mass_white_dwarf;
+  }
+
+  return (PyObject *) yields_out;
+  
+}
+
+/**
+ * @brief Compute the mass ejected by a supernovae Ia.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction ejected.
+ */
+PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "s", &filename))
+    return NULL;
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* return object */
+  const int ndims = 1;
+  npy_intp dims[1] = {1};
+  PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+  float *y = (float*) PyArray_GETPTR1(mass_ejected, 0);
+
+  *y = supernovae_ia_get_ejected_mass_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
+
+  return (PyObject *) mass_ejected;
+  
+}
+
+
+/**
+ * @brief Read the elements names.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pystellar_get_elements(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(args, "s", &filename))
+    return NULL;
+
+  int with_cosmo = 0;
+
+  PYSWIFTSIM_INIT_STRUCTS();
+  
+  /* Init stellar physics */
+  struct feedback_props fp;
+  feedback_props_init(&fp, &pconst, &us, &params, &hydro_props, &cosmo);
+
+  /* Save the variables */
+  PyObject *t = PyList_New(CHEMISTRY_ELEMENT_COUNT);
+
+  for(int i = 0; i < CHEMISTRY_ELEMENT_COUNT; i++) {
+    const char *name = stellar_evolution_get_element_name(&fp.stellar_model, i);
+    PyObject *string = PyUnicode_FromString(name);
+    int test = PyList_SetItem(t, i, string);
+    if (test != 0) {
+      error("Failed to add an element");
+    }
+  }
+
+  return t;
+  
+}

src/stellar_evolution_wrapper.h

@@ -26,6 +26,13 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args);
 PyObject* pymass_from_lifetime(PyObject* self, PyObject* args);
 PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args);
 PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* args);
+PyObject* pystellar_get_elements(PyObject* self, PyObject* args);
+
 
 #endif // __PYSWIFTSIM_STELLAR_H__
 

tests/test_yields.py

@@ -64,6 +64,7 @@ if __name__ == "__main__":
         m1 /= solMass_in_cgs / unit_mass
 
         yields_snia = libstellar.supernovaeIaYields(filename)
+        mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
 
         print("Mass used: {}".format(m1))
 
@@ -73,5 +74,7 @@ if __name__ == "__main__":
 
         print("SNII mass ejected: {}".format(mass_ejected_snii))
 
+        print("SNIa mass ejected: {}".format(mass_ejected_snia))
+
         print("SNII non processed mass ejected: {}".format(
             mass_ejected_np_snii))