diff --git a/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
index ce3d510430474c3986b463878c3edff2ba02d212..c140dca3d370df398e77fa31d871e97e9652c7a6 100644
--- a/examples/cooling_rate/run.sh
+++ b/examples/cooling_rate/run.sh
@@ -1,4 +1,4 @@
#!/bin/bash
./plot_cooling.py 1
-./plot_cooling.py 100
+# ./plot_cooling.py 100
diff --git a/examples/initial_mass_function/plot_initial_mass_function.py b/examples/initial_mass_function/plot_initial_mass_function.py
index bc653edc70a36756b7f65f206fd4e5e02cd56129..bb972399eee6e1a756321695e5e63326c3e0b6d5 100644
--- a/examples/initial_mass_function/plot_initial_mass_function.py
+++ b/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()
diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py
index 7696311e3cfd048d2dd000bebbbfc4389d448125..79998ef4f98eed12a4c78efd94809c0861208011 100644
--- a/examples/lifetime/plot_lifetime.py
+++ b/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")
diff --git a/examples/supernovae_rate/plot_rates.py b/examples/supernovae_rate/plot_rates.py
index 0e4943b3b96324764ee5eb19af812cbd2e357e54..300c5fcb314d79796d7fb54e67ea49070dfe2ea1 100644
--- a/examples/supernovae_rate/plot_rates.py
+++ b/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()
diff --git a/scripts/pyswift_cooling_rate b/scripts/pyswift_cooling_rate
index 262629c6debc640b90821d828f9190ef5b67d872..7c905944b6bda576c1048fcc11828912810c8f9e 100644
--- a/scripts/pyswift_cooling_rate
+++ b/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()
diff --git a/scripts/pyswift_initial_mass_function b/scripts/pyswift_initial_mass_function
index 803e520495fc9d453963975de2e4100d20374bc5..fa35ce5a8d50329fe5b9fc0bbba2ffcc932ae667 100644
--- a/scripts/pyswift_initial_mass_function
+++ b/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)
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
index 8485f006cd4e36deba35ceda434c6f7bd9785178..9fb7754db9c4ceadde0fedaf4007c375c2c1a7ee 100644
--- a/scripts/pyswift_lifetime
+++ b/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)
diff --git a/scripts/pyswift_supernovae_rate b/scripts/pyswift_supernovae_rate
index 65bf52b12a6c9c6d6dbf588038ad368ac2073455..a60baf0f6d493e59964cd441df2c9e2aff914be5 100644
--- a/scripts/pyswift_supernovae_rate
+++ b/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}$]")
diff --git a/tests/test_initial_mass_function.py b/tests/test_initial_mass_function.py
index ca80c505bb39844035953934cdb84bf221591068..416aa3c4f7905b419200358cd9632bc1af4513db 100644
--- a/tests/test_initial_mass_function.py
+++ b/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))
diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py
index 3fc7a84d9a9a91eedaf7ccc033058d88bb6e93cf..aac4281d47c93cf5cb7ac0c1262e119736d25b8c 100644
--- a/tests/test_lifetime.py
+++ b/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))
diff --git a/tests/test_rates.py b/tests/test_rates.py
index 332f957001b460a1feb1913e1ffbbf7fcec4cc0f..7bc38ff4f37d96e98fc822225fb4c594cc367ab5 100644
--- a/tests/test_rates.py
+++ b/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))