From b91850cb7fe4394d0dabf51995c8f5f9b609c379 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Fri, 7 Jun 2019 09:52:02 +0200
Subject: [PATCH] Add units
---
examples/cooling_rate/run.sh | 2 +-
.../plot_initial_mass_function.py | 28 +++++++----
examples/lifetime/plot_lifetime.py | 21 ++++++++-
examples/supernovae_rate/plot_rates.py | 29 ++++++++++--
scripts/pyswift_cooling_rate | 3 --
scripts/pyswift_initial_mass_function | 11 ++++-
scripts/pyswift_lifetime | 18 ++++++-
scripts/pyswift_supernovae_rate | 47 +++++++++++++++----
tests/test_initial_mass_function.py | 10 +++-
tests/test_lifetime.py | 15 +++++-
tests/test_rates.py | 20 +++++++-
11 files changed, 171 insertions(+), 33 deletions(-)
diff --git a/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
index ce3d510..c140dca 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 bc653ed..bb97239 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 7696311..79998ef 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 0e4943b..300c5fc 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 262629c..7c90594 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 803e520..fa35ce5 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 8485f00..9fb7754 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 65bf52b..a60baf0 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 ca80c50..416aa3c 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 3fc7a84..aac4281 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 332f957..7bc38ff 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))
--
GitLab