diff --git a/docs/RST/examples/stellar.rst b/docs/RST/examples/stellar.rst
index 94e1ac8ee917b03ae0ee1fd50d3886c07f312707..031a8849b978958ec5b183d29be70ae421c16cd0 100644
--- a/docs/RST/examples/stellar.rst
+++ b/docs/RST/examples/stellar.rst
@@ -21,5 +21,14 @@ The lifetime in GEAR is a quadratic in metallicity and mass (see equation 3.32 i
 .. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/lifetime.png
 
 
+Supernovae Rate
+~~~~~~~~~~~~~~~
+
+This example is generated in ``examples/supernovae_rate`` with GEAR's stellar model.
+The supernovae rate is given by the contribution of SNIa (in blue) and SNII (in orange).
+The two peaks in SNIa correspond to the two different type of companion (red giant and main sequence).
+
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/supernovae_rate.png
+
 .. [Kroupa2001] `On the variation of the Initial Mass Function <https://arxiv.org/abs/astro-ph/0009005>`_
 .. [Poirier2004] `Étude de l'évolution chimique et dynamique d'objets proto-galactiques: Application à l'évolution des galaxies spirales <https://obswww.unige.ch/~revaz/PyChem/_downloads/ff03ba21204a3f36fa4ec885abd545ae/ThesePoirier.pdf>`_
diff --git a/examples/supernovae_rate/plot_rates.py b/examples/supernovae_rate/plot_rates.py
new file mode 100644
index 0000000000000000000000000000000000000000..0e4943b3b96324764ee5eb19af812cbd2e357e54
--- /dev/null
+++ b/examples/supernovae_rate/plot_rates.py
@@ -0,0 +1,62 @@
+#!/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
+plt.rc('text', usetex=True)
+
+
+N = 1000
+solar_metallicity = 0.02041
+t_min = 1e6  # in yr
+t_max = 3e10
+
+
+if __name__ == "__main__":
+    with pyswiftsim.ParameterFile() as filename:
+
+        # Get times
+        times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
+                            dtype=np.float32)
+        t1 = times[:-1]
+        t2 = times[1:]
+
+        # get metallicity
+        Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+        print("Computing supernovae rate...")
+
+        rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+        rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+        t1 /= 1e6
+        rate_snia *= 1e9
+        rate_snii *= 1e9
+
+        plt.loglog(t1, rate_snia, label="SNIa")
+        plt.loglog(t1, rate_snii, label="SNII")
+
+        plt.xlabel("Time [Myr]")
+        plt.ylabel("Supernovae rate [$M_\odot^{-1} Gyr^{-1}$]")
+        plt.legend()
+
+        plt.show()
diff --git a/examples/supernovae_rate/run.sh b/examples/supernovae_rate/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..4264cedc3634a1af6e88d23fc110956587a0932b
--- /dev/null
+++ b/examples/supernovae_rate/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python plot_rates.py
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 7852bc713c33d950db8883ca686b71f12e998b4f..6649bc3091ac65aec3aefa91e5d7dbfe3512a0fe 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -78,7 +78,7 @@ def parseYamlFile(filename):
     >>> unit_mass = float(us["UnitMass_in_cgs"])
     """
     with open(filename, "r") as stream:
-        stream = yaml.load(stream)
+        stream = yaml.load(stream, Loader=yaml.SafeLoader)
         return stream
 
 
@@ -158,7 +158,7 @@ class ParameterFile:
         "GEARInitialMassFunction": {
             "number_function_part": 4,
             "exponents": [0.7, -0.8, -1.7, -1.3],
-            "mass_limits": [0.05, 0.08, 0.5, 1.0, 50.]
+            "mass_limits_msun": [0.05, 0.08, 0.5, 1.0, 50.]
         },
 
         "GEARFeedback": {
@@ -171,6 +171,23 @@ class ParameterFile:
             "quadratic": [-40.110, 5.509, 0.7824],
             "linear": [141.929, -15.889, -3.2557],
             "constant": [-261.365, 17.073, 9.8661],
+        },
+
+        "GEARSupernovaeIa": {
+            "exponent": -0.35,
+            "min_mass_white_dwarf_progenitor": 3.,
+            "max_mass_white_dwarf_progenitor": 8.,
+            "min_mass_red_giant": 0.9,
+            "max_mass_red_giant": 1.5,
+            "coef_red_giant": 0.02,
+            "min_mass_main_sequence": 1.8,
+            "max_mass_main_sequence": 2.5,
+            "coef_main_sequence": 0.05,
+        },
+
+        "GEARSupernovaeII": {
+            "min_mass": 8.,
+            "max_mass": 50.,
         }
     }
     """
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
index 0b14c53f589d9a26c21841b4393536a000dc1eec..8485f006cd4e36deba35ceda434c6f7bd9785178 100644
--- a/scripts/pyswift_lifetime
+++ b/scripts/pyswift_lifetime
@@ -28,7 +28,7 @@ plt.rc('text', usetex=True)
 
 
 def parseArguments():
-    parser = ArgumentParser(description="Plot the initial mass function")
+    parser = ArgumentParser(description="Plot the lifetime of a star")
 
     parser.add_argument(
         "--download_table", dest="table",
diff --git a/scripts/pyswift_supernovae_rate b/scripts/pyswift_supernovae_rate
new file mode 100644
index 0000000000000000000000000000000000000000..65bf52b12a6c9c6d6dbf588038ad368ac2073455
--- /dev/null
+++ b/scripts/pyswift_supernovae_rate
@@ -0,0 +1,109 @@
+#!/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 argparse import ArgumentParser
+import os
+plt.rc('text', usetex=True)
+
+
+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(
+        "--time_min", dest="time_min",
+        type=float, default=1.,
+        help="Minimal time to plot in Myr")
+
+    parser.add_argument(
+        "--time_max", dest="time_max",
+        type=float, default=3e4,
+        help="Maximal time to plot in Myr")
+
+    parser.add_argument(
+        "--metallicity", dest="metallicity",
+        type=float, default=0.02041,
+        help="Metallicity of the stars")
+
+    args = parser.parse_args()
+
+    return args
+
+
+if __name__ == "__main__":
+    opt = parseArguments()
+
+    if opt.table:
+        pyswiftsim.downloadChemistryTable()
+
+    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:
+
+        # Get times
+        opt.time_min *= 1e6
+        opt.time_max *= 1e6
+
+        times = np.logspace(
+            np.log10(opt.time_min), np.log10(opt.time_max), opt.N + 1,
+            dtype=np.float32)
+        t1 = times[:-1]
+        t2 = times[1:]
+
+        # get metallicity
+        Z = opt.metallicity * np.ones(opt.N, dtype=np.float32)
+
+        print("Computing supernovae rate...")
+
+        rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+        rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+        t1 /= 1e6
+        rate_snia *= 1e9
+        rate_snii *= 1e9
+
+        plt.loglog(t1, rate_snia, label="SNIa")
+        plt.loglog(t1, rate_snii, label="SNII")
+
+        plt.xlabel("Time [Myr]")
+        plt.ylabel("Supernovae rate [$M_\odot^{-1} Gyr^{-1}$]")
+        plt.legend()
+
+        plt.show()
diff --git a/src/libstellar.c b/src/libstellar.c
index fcd1cb304d5207fd409eb9461337e0bac03ac60e..f9957d08ac626ee5ba5e57402906d2c3ebdf6c93 100644
--- a/src/libstellar.c
+++ b/src/libstellar.c
@@ -96,6 +96,60 @@ static PyMethodDef libstellar_methods[] = {
    ">>> mass = libstellar.massFromLifetime(filename, lifetime, Z)\n"
   },
 
+  {"supernovaeIaRate", pysupernovae_ia_rate, METH_VARARGS,
+   "Compute the supernovae Ia rate from the time and the stellar metallicity.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "time1: array\n"
+   "\t Beginning of the time interval for the SNIa in internal units.\n"
+   "time2: array\n"
+   "\t End of the time interval for the SNIa in internal units.\n"
+   "metallicity: array\n"
+   "\t Metallicity of the stars.\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "rate: array\n"
+   "\t The rate of SNIa.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> time1 = np.array([1e7, 1e9])\n"
+   ">>> time2 = np.array([2e7, 2e9])\n"
+   ">>> Z = np.array([0.02, 0.01])\n"
+   ">>> filename = 'params.yml'\n"
+   ">>> mass = libstellar.supernovaeIaRate(filename, time1, time2, Z)\n"
+  },
+
+  {"supernovaeIIRate", pysupernovae_ii_rate, METH_VARARGS,
+   "Compute the supernovae II rate from the time and the stellar metallicity.\n\n"
+   "Parameters\n"
+   "----------\n\n"
+   "filename: str\n"
+   "\t Parameter file\n"
+   "time1: array\n"
+   "\t Beginning of the time interval for the SNIa in internal units.\n"
+   "time2: array\n"
+   "\t End of the time interval for the SNIa in internal units.\n"
+   "metallicity: array\n"
+   "\t Metallicity of the stars.\n"
+   "\n"
+   "Returns\n"
+   "-------\n"
+   "rate: array\n"
+   "\t The rate of SNII.\n"
+   "\n"
+   "Examples\n"
+   "--------\n"
+   ">>> time1 = np.array([1e7, 1e9])\n"
+   ">>> time2 = np.array([2e7, 2e9])\n"
+   ">>> Z = np.array([0.02, 0.01])\n"
+   ">>> filename = 'params.yml'\n"
+   ">>> mass = libstellar.supernovaeIIRate(filename, time1, time2, Z)\n"
+  },
+
 
   {NULL, NULL, 0, NULL}        /* Sentinel */
 };      
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
index a990b753d2c05ddc620f2115bce366a75fc4ce09..bbd31d41bd5a557ae7677f2e503ea94abffc9a38 100644
--- a/src/stellar_evolution_wrapper.c
+++ b/src/stellar_evolution_wrapper.c
@@ -75,11 +75,13 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
  * @brief Compute the lifetime of a star.
  *
  * args is expecting:
- * TODO
+ *  - the name of the parameter file
+ *  - a numpy array containing the masses.
+ *  - a numpy array containing the metallicity.
  *
  * @param self calling object
  * @param args arguments
- * @return cooling rate
+ * @return the lifetime of the stars
  */
 PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
   import_array();
@@ -135,11 +137,13 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
  * @brief Compute the mass of a star with a given lifetime.
  *
  * args is expecting:
- * TODO
+ *  - the name of the parameter file
+ *  - a numpy array containing the time.
+ *  - a numpy array containing the metallicity.
  *
  * @param self calling object
  * @param args arguments
- * @return cooling rate
+ * @return The mass of the stars.
  */
 PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
   import_array();
@@ -190,3 +194,155 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
   return (PyObject *) mass;
   
 }
+
+
+
+/**
+ * @brief Compute the supernovae Ia rate.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *  - a numpy array containing the lower limit for the time.
+ *  - a numpy array containing the upper limit for the time.
+ *  - a numpy array containing the metallicity.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  PyArrayObject *time1;
+  PyArrayObject *time2;
+  PyArrayObject *metallicity;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "sOOO", &filename, &time1, &time2, &metallicity))
+    return NULL;
+
+  /* check numpy array */
+  if (pytools_check_array(time1, 1, NPY_FLOAT, "Time1") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(time2, 1, NPY_FLOAT, "Time2") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+    return NULL;
+
+  if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time1, 0))
+    error("Time1 and metallicity should have the same dimension");
+  if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time2, 0))
+    error("Time2 and metallicity should have the same dimension");
+
+
+  size_t N = PyArray_DIM(time1, 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 *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(time1), PyArray_DIMS(time1), NPY_FLOAT);
+
+  for(size_t i = 0; i < N; i++) {
+    float t1 = *(float*) PyArray_GETPTR1(time1, i);
+    float t2 = *(float*) PyArray_GETPTR1(time2, i);
+    float z = *(float*) PyArray_GETPTR1(metallicity, i);
+    float *r = (float*) PyArray_GETPTR1(rate, i);
+
+    
+    
+    float m1 = stellar_evolution_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);
+    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);
+  }
+
+  return (PyObject *) rate;
+  
+}
+
+/**
+ * @brief Compute the supernovae II rate.
+ *
+ * args is expecting:
+ *  - the name of the parameter file
+ *  - a numpy array containing the lower limit for the time.
+ *  - a numpy array containing the upper limit for the time.
+ *  - a numpy array containing the metallicity.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
+  import_array();
+  
+  char *filename;
+
+  PyArrayObject *time1;
+  PyArrayObject *time2;
+  PyArrayObject *metallicity;
+
+  /* parse argument */
+  if (!PyArg_ParseTuple(
+      args, "sOOO", &filename, &time1, &time2, &metallicity))
+    return NULL;
+
+  /* check numpy array */
+  if (pytools_check_array(time1, 1, NPY_FLOAT, "Time1") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(time2, 1, NPY_FLOAT, "Time2") != SWIFT_SUCCESS)
+    return NULL;
+  if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+    return NULL;
+
+  if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time1, 0))
+    error("Time1 and metallicity should have the same dimension");
+  if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time2, 0))
+    error("Time2 and metallicity should have the same dimension");
+
+
+  size_t N = PyArray_DIM(time1, 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 *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(time1), PyArray_DIMS(time1), NPY_FLOAT);
+
+  for(size_t i = 0; i < N; i++) {
+    float t1 = *(float*) PyArray_GETPTR1(time1, i);
+    float t2 = *(float*) PyArray_GETPTR1(time2, i);
+    float z = *(float*) PyArray_GETPTR1(metallicity, i);
+    float *r = (float*) PyArray_GETPTR1(rate, i);
+
+    
+    
+    float m1 = stellar_evolution_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);
+    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);
+  }
+
+  return (PyObject *) rate;
+  
+}
diff --git a/src/stellar_evolution_wrapper.h b/src/stellar_evolution_wrapper.h
index a76b04d39b1733a8667cc9435859cfeb9fbfbeea..8cce237edc04f4ba1f20662ec2197a36411d69d8 100644
--- a/src/stellar_evolution_wrapper.h
+++ b/src/stellar_evolution_wrapper.h
@@ -24,6 +24,8 @@
 PyObject* pyinitial_mass_function(PyObject* self, PyObject* args);
 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);
 
 #endif // __PYSWIFTSIM_STELLAR_H__
 
diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py
index 9a1816e22438b5a857a9945e571961af280cccbc..3fc7a84d9a9a91eedaf7ccc033058d88bb6e93cf 100644
--- a/tests/test_lifetime.py
+++ b/tests/test_lifetime.py
@@ -35,7 +35,7 @@ if __name__ == "__main__":
         parser = pyswiftsim.parseYamlFile(filename)
         # Get masses
         imf = parser["GEARInitialMassFunction"]
-        mass_limits = imf["mass_limits"]
+        mass_limits = imf["mass_limits_msun"]
         mass_min = float(mass_limits[0])
         mass_max = float(mass_limits[1])
 
diff --git a/tests/test_rates.py b/tests/test_rates.py
new file mode 100644
index 0000000000000000000000000000000000000000..332f957001b460a1feb1913e1ffbbf7fcec4cc0f
--- /dev/null
+++ b/tests/test_rates.py
@@ -0,0 +1,53 @@
+#!/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
+plt.rc('text', usetex=True)
+
+
+N = 100
+solar_metallicity = 0.02041
+t_min = 1e6 # in yr
+t_max = 1e10
+
+
+if __name__ == "__main__":
+    with pyswiftsim.ParameterFile() as filename:
+
+        # Get times
+        times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
+                            dtype=np.float32)
+        t1 = times[:-1]
+        t2 = times[1:]
+
+        # get metallicity
+        Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+        print("Computing supernovae rate...")
+
+        rate = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+        rate += libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+        print("Time used: {}".format(t1))
+
+        print("Rate obtained: {}".format(rate))