Add units

examples/cooling_rate/run.sh

 #!/bin/bash
 
 ./plot_cooling.py 1
-./plot_cooling.py 100
+# ./plot_cooling.py 100

examples/initial_mass_function/plot_initial_mass_function.py

@@ -22,32 +22,42 @@ import pyswiftsim
 
 import numpy as np
 import matplotlib.pyplot as plt
+from time import strftime, gmtime
 plt.rc('text', usetex=True)
 
 N = 100
-
+solMass_in_cgs = 1.98848e33
 
 if __name__ == "__main__":
-
     with pyswiftsim.ParameterFile() as filename:
 
         parser = pyswiftsim.parseYamlFile(filename)
+
+        # get parser
         imf_parser = parser["GEARInitialMassFunction"]
+        units = parser["InternalUnitSystem"]
+        unit_mass = float(units["UnitMass_in_cgs"])
+        solMass_in_internal = solMass_in_cgs / unit_mass
 
-        mass_limits = list(imf_parser["mass_limits"])
-        mass_min = np.log10(mass_limits[0])
-        mass_max = np.log10(mass_limits[-1])
+        mass_limits = list(imf_parser["mass_limits_msun"])
+        mass_min = np.log10(mass_limits[0] * solMass_in_internal)
+        mass_max = np.log10(mass_limits[-1] * solMass_in_internal)
+
+        mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
 
-        # du / dt
         print("Computing initial mass function...")
 
-        mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
         imf = libstellar.initialMassFunction(filename, mass)
+        imf *= solMass_in_internal
 
         # plot IMF
+        mass /= solMass_in_internal
         plt.loglog(mass, imf)
 
-        plt.xlabel("Mass [Msun]")
-        plt.ylabel("Mass fraction")
+        plt.figure()
+        date = strftime("%d %b %Y", gmtime())
+        plt.title("Date: {}".format(date))
+        plt.xlabel("Mass [$M_\mathrm{sun}$]")
+        plt.ylabel("Initial Mass Function [$M_\mathrm{sun}^{-1}$]")
 
         plt.show()

examples/lifetime/plot_lifetime.py

@@ -22,11 +22,14 @@ import pyswiftsim
 
 import numpy as np
 import matplotlib.pyplot as plt
+from time import strftime, gmtime
 plt.rc('text', usetex=True)
 
 
 N = 1000
 solar_metallicity = 0.02041
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
 
 
 if __name__ == "__main__":
@@ -35,12 +38,21 @@ 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 = np.log10(mass_limits[0])
         mass_max = np.log10(mass_limits[-1])
 
         mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
 
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        mass *= solMass_in_cgs / unit_mass
+
         # get metallicity
         Z = solar_metallicity * np.ones(N, dtype=np.float32)
 
@@ -49,8 +61,13 @@ if __name__ == "__main__":
         lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
         mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
 
-        lifetime /= 1e6
+        mass /= solMass_in_cgs / unit_mass
+        mass_from_lifetime /= solMass_in_cgs / unit_mass
+        lifetime /= 1e6 * year_in_cgs / unit_time
 
+        plt.figure()
+        date = strftime("%d %b %Y", gmtime())
+        plt.title("Date: {}".format(date))
         plt.loglog(mass, lifetime, label="Lifetime from mass")
         plt.loglog(mass_from_lifetime, lifetime, "--",
                    label="Mass from lifetime")

examples/supernovae_rate/plot_rates.py

@@ -22,6 +22,7 @@ import pyswiftsim
 
 import numpy as np
 import matplotlib.pyplot as plt
+from time import strftime, gmtime
 plt.rc('text', usetex=True)
 
 
@@ -29,14 +30,27 @@ N = 1000
 solar_metallicity = 0.02041
 t_min = 1e6  # in yr
 t_max = 3e10
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
 
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
 
         # Get times
         times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
                             dtype=np.float32)
+
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        times *= year_in_cgs / unit_time
+
         t1 = times[:-1]
         t2 = times[1:]
 
@@ -48,15 +62,22 @@ if __name__ == "__main__":
         rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
         rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
 
-        t1 /= 1e6
-        rate_snia *= 1e9
-        rate_snii *= 1e9
+        t1 /= 1e6 * year_in_cgs / unit_time
+
+        unit_rate = year_in_cgs / unit_time
+        unit_rate *= solMass_in_cgs / unit_mass
+
+        rate_snia *= 1e9 * unit_rate
+        rate_snii *= 1e9 * unit_rate
 
+        plt.figure()
+        date = strftime("%d %b %Y", gmtime())
+        plt.title("Date: {}".format(date))
         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.ylabel("Supernovae rate [$M_\odot^{-1} \mathrm{Gyr}^{-1}$]")
         plt.legend()
 
         plt.show()

scripts/pyswift_cooling_rate

@@ -25,7 +25,6 @@ import numpy as np
 import matplotlib.pyplot as plt
 from astropy import units, constants
 import yaml
-from time import strftime, gmtime
 from argparse import ArgumentParser
 import os
 
@@ -244,8 +243,6 @@ def plot_solution(rate, data, us, opt):
         cbar.set_ticklabels(tc)
         cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
 
-    date = strftime("%d %b %Y", gmtime())
-    plt.title("Date: {}".format(date))
     plt.show()
 
 

scripts/pyswift_initial_mass_function

@@ -26,6 +26,7 @@ from argparse import ArgumentParser
 import os
 
 plt.rc('text', usetex=True)
+solMass_in_cgs = 1.98848e33
 
 
 def parseArguments():
@@ -75,8 +76,10 @@ if __name__ == "__main__":
 
         parser = pyswiftsim.parseYamlFile(filename)
         imf_parser = parser["GEARInitialMassFunction"]
+        units = parser["InternalUnitSystem"]
+        unit_mass = float(units["UnitMass_in_cgs"])
 
-        mass_limits = list(imf_parser["mass_limits"])
+        mass_limits = list(imf_parser["mass_limits_msun"])
         if opt.mass_min is not None:
             mass_min = np.log10(opt.mass_min)
         else:
@@ -87,11 +90,17 @@ if __name__ == "__main__":
         else:
             mass_max = np.log10(mass_limits[-1])
 
+        mass_min *= solMass_in_cgs / unit_mass
+        mass_max *= solMass_in_cgs / unit_mass
+
         print("Computing initial mass function...")
 
         mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
         imf = libstellar.initialMassFunction(filename, mass)
 
+        imf *= solMass_in_cgs / unit_mass
+        mass /= solMass_in_cgs / unit_mass
+
         # plot IMF
         plt.loglog(mass, imf)
 

scripts/pyswift_lifetime

@@ -26,6 +26,9 @@ from argparse import ArgumentParser
 import os
 plt.rc('text', usetex=True)
 
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
+
 
 def parseArguments():
     parser = ArgumentParser(description="Plot the lifetime of a star")
@@ -80,7 +83,7 @@ if __name__ == "__main__":
         parser = pyswiftsim.parseYamlFile(filename)
         imf_parser = parser["GEARInitialMassFunction"]
 
-        mass_limits = list(imf_parser["mass_limits"])
+        mass_limits = list(imf_parser["mass_limits_msun"])
         if opt.mass_min is not None:
             mass_min = np.log10(opt.mass_min)
         else:
@@ -93,13 +96,24 @@ if __name__ == "__main__":
 
         mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
 
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        mass *= solMass_in_cgs / unit_mass
+
         # get metallicity
         Z = opt.metallicity * np.ones(opt.N, dtype=np.float32)
 
         print("Computing lifetime...")
 
         lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
-        lifetime /= 1e6
+
+        mass /= solMass_in_cgs / unit_mass
+        lifetime /= 1e6 * year_in_cgs / unit_time
 
         plt.loglog(mass, lifetime)
 

scripts/pyswift_supernovae_rate

@@ -26,6 +26,9 @@ from argparse import ArgumentParser
 import os
 plt.rc('text', usetex=True)
 
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
+
 
 def parseArguments():
     parser = ArgumentParser(description="Plot the supernovae rate")
@@ -60,6 +63,16 @@ def parseArguments():
         type=float, default=0.02041,
         help="Metallicity of the stars")
 
+    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
@@ -76,14 +89,25 @@ if __name__ == "__main__":
             "The parameter file '{}' does not exist".format(opt.params))
 
     with pyswiftsim.ParameterFile(existing=opt.params) as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
 
-        # Get times
+        # Get times in yr
         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)
+
+        # 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
+
+        times *= year_in_cgs / unit_time
+
         t1 = times[:-1]
         t2 = times[1:]
 
@@ -92,15 +116,22 @@ if __name__ == "__main__":
 
         print("Computing supernovae rate...")
 
-        rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
-        rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+        if not opt.disable_snia:
+            rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+        if not opt.disable_snii:
+            rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+        t1 /= 1e6 * year_in_cgs / unit_time
 
-        t1 /= 1e6
-        rate_snia *= 1e9
-        rate_snii *= 1e9
+        unit_rate = year_in_cgs / unit_time
+        unit_rate *= solMass_in_cgs / unit_mass
 
-        plt.loglog(t1, rate_snia, label="SNIa")
-        plt.loglog(t1, rate_snii, label="SNII")
+        if not opt.disable_snia:
+            rate_snia *= 1e9 * unit_rate
+            plt.loglog(t1, rate_snia, label="SNIa")
+        if not opt.disable_snii:
+            rate_snii *= 1e9 * unit_rate
+            plt.loglog(t1, rate_snii, label="SNII")
 
         plt.xlabel("Time [Myr]")
         plt.ylabel("Supernovae rate [$M_\odot^{-1} Gyr^{-1}$]")

tests/test_initial_mass_function.py

@@ -26,6 +26,7 @@ plt.rc('text', usetex=True)
 
 
 N = 100
+solMass_in_cgs = 1.989e33
 
 
 if __name__ == "__main__":
@@ -34,14 +35,21 @@ if __name__ == "__main__":
 
         parser = pyswiftsim.parseYamlFile(filename)
         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])
 
         mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
 
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        # change the units
+        mass *= solMass_in_cgs / unit_mass
+
         print("Computing initial mass function...")
 
         imf = libstellar.initialMassFunction(filename, mass)
 
+        imf *= solMass_in_cgs / unit_mass
+
         print("IMF obtained: {}".format(imf))

tests/test_lifetime.py

@@ -27,7 +27,8 @@ plt.rc('text', usetex=True)
 
 N = 100
 solar_metallicity = 0.02041
-
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
@@ -41,6 +42,14 @@ if __name__ == "__main__":
 
         mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
 
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+        # change the units
+        mass *= solMass_in_cgs / unit_mass
+
         # get metallicity
         Z = solar_metallicity * np.ones(N, dtype=np.float32)
 
@@ -49,6 +58,10 @@ if __name__ == "__main__":
         lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
         mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
 
+        mass /= solMass_in_cgs / unit_mass
+        mass_from_lifetime /= solMass_in_cgs / unit_mass
+        lifetime /= year_in_cgs / unit_time
+
         print("Mass used: {}".format(mass))
 
         print("Lifetime obtained: {}".format(lifetime))

tests/test_rates.py

@@ -27,16 +27,29 @@ plt.rc('text', usetex=True)
 
 N = 100
 solar_metallicity = 0.02041
-t_min = 1e6 # in yr
+t_min = 1e6  # in yr
 t_max = 1e10
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
 
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
 
         # Get times
         times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
                             dtype=np.float32)
+
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        times *= year_in_cgs / unit_time
+
         t1 = times[:-1]
         t2 = times[1:]
 
@@ -48,6 +61,11 @@ if __name__ == "__main__":
         rate = libstellar.supernovaeIaRate(filename, t1, t2, Z)
         rate += libstellar.supernovaeIIRate(filename, t1, t2, Z)
 
+        rate *= year_in_cgs / unit_time
+        rate *= solMass_in_cgs / unit_mass
+
+        t1 /= year_in_cgs / unit_time
+
         print("Time used: {}".format(t1))
 
         print("Rate obtained: {}".format(rate))