Change units

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...")
 

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;