diff --git a/examples/initial_mass_function/plot_initial_mass_function.py b/examples/initial_mass_function/plot_initial_mass_function.py
index 6957770991a3998a304179bd5a005fc84b4a7a9e..873aa4bb6bedaa64f0971476bcbbd2ff5e9fd681 100644
--- a/examples/initial_mass_function/plot_initial_mass_function.py
+++ b/examples/initial_mass_function/plot_initial_mass_function.py
@@ -45,7 +45,8 @@ if __name__ == "__main__":
mass_min = np.log10(mass_min * solMass_in_internal)
mass_max = np.log10(mass_max * solMass_in_internal)
- mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
+ mass = np.logspace(mass_min, mass_max, N, dtype=np.float32,
+ endpoint=False)
print("Computing initial mass function...")
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
index d38bbe63fda08fc2dc6633a901f57a9a87b61b7c..007e369a8a5550e95d3d15cdfbb0ab2fb882d22d 100644
--- a/src/stellar_evolution_wrapper.c
+++ b/src/stellar_evolution_wrapper.c
@@ -63,6 +63,13 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
float m = *(float*) PyArray_GETPTR1(mass, i);
float *imf_i = (float*) PyArray_GETPTR1(imf, i);
+ /* change units internal units -> solar mass */
+ m /= pconst.const_solar_mass;
+
+ if (m > fp.stellar_model.imf.mass_max) {
+ *imf_i = 0.;
+ continue;
+ }
*imf_i = initial_mass_function_get_imf(&fp.stellar_model.imf, m);
}
@@ -124,9 +131,12 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
float z = *(float*) PyArray_GETPTR1(metallicity, i);
float *tau = (float*) PyArray_GETPTR1(lifetime, i);
+
+ /* change units internal units -> solar mass */
+ m /= pconst.const_solar_mass;
*tau = lifetime_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
- *tau = pow(10, *tau);
+ *tau = pow(10, *tau) * pconst.const_year * 1e6;
}
return (PyObject *) lifetime;
@@ -187,8 +197,11 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
float *m = (float*) PyArray_GETPTR1(mass, i);
+ /* change units internal units to Myr */
+ tau /= pconst.const_year * 1e6;
+
*m = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
- *m = pow(10, *m);
+ *m = pow(10, *m) * pconst.const_solar_mass;
}
return (PyObject *) mass;
@@ -257,7 +270,9 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
float z = *(float*) PyArray_GETPTR1(metallicity, i);
float *r = (float*) PyArray_GETPTR1(rate, i);
-
+ /* change units internal units -> Myr */
+ t1 /= pconst.const_year * 1e6;
+ t2 /= pconst.const_year * 1e6;
float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
m1 = pow(10, m1);
@@ -266,6 +281,9 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
/* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
*r = supernovae_ia_get_number(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
+
+ /* change units 1 / (solMass Myr) -> internal units */
+ *r /= pconst.const_solar_mass * pconst.const_year * 1e6;
}
return (PyObject *) rate;
@@ -332,7 +350,9 @@ PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
float z = *(float*) PyArray_GETPTR1(metallicity, i);
float *r = (float*) PyArray_GETPTR1(rate, i);
-
+ /* change units internal units -> Myr */
+ t1 /= pconst.const_year * 1e6;
+ t2 /= pconst.const_year * 1e6;
float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
m1 = pow(10, m1);
@@ -341,6 +361,10 @@ PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
/* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
*r = supernovae_ii_get_number(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
+
+ /* change units 1 / (solMass Myr) -> internal units */
+ *r /= pconst.const_solar_mass * pconst.const_year * 1e6;
+
}
return (PyObject *) rate;
@@ -402,9 +426,9 @@ PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
for(size_t i = 0; i < N; i++) {
float yields[CHEMISTRY_ELEMENT_COUNT] = {0.};
- const float m1 = *(float*) PyArray_GETPTR1(mass1, i);
+ const float m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
const float log_m1 = log10(m1);
- const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
const float log_m2 = log10(m2);
/* Compute the yields */
@@ -473,13 +497,13 @@ PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args) {
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 m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
const float log_m1 = log10(m1);
- const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
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);
+ *y = supernovae_ii_get_ejected_mass_fraction_processed(&fp.stellar_model.snii, log_m1, log_m2);
}
@@ -540,13 +564,13 @@ PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* a
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 m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
const float log_m1 = log10(m1);
- const float m2 = *(float*) PyArray_GETPTR1(mass2, i);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
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);
+ *y = supernovae_ii_get_ejected_mass_fraction(&fp.stellar_model.snii, log_m1, log_m2);
}
@@ -634,7 +658,7 @@ PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args) {
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;
+ *y = supernovae_ia_get_ejected_mass_fraction_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
return (PyObject *) mass_ejected;