diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py new file mode 100644 index 0000000000000000000000000000000000000000..a1e1e9f2fec9e23a1d8bef49cd95361046617a9d --- /dev/null +++ b/examples/lifetime/plot_lifetime.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +# ######################################################################## +# This file is part of PYSWIFT. +# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# ######################################################################## + +from pyswiftsim import libstellar +import pyswiftsim + +import numpy as np +import matplotlib.pyplot as plt +plt.rc('text', usetex=True) + + +N = 1000 +solar_metallicity = 0.02041 + + +if __name__ == "__main__": + with pyswiftsim.ParameterFile() as filename: + + parser = pyswiftsim.parseYamlFile(filename) + # Get masses + imf = parser["GEARInitialMassFunction"] + mass_limits = imf["mass_limits"] + 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) + + # get metallicity + Z = solar_metallicity * np.ones(N, dtype=np.float32) + + print("Computing lifetime...") + + lifetime = libstellar.lifetimeFromMass(filename, mass, Z) + lifetime /= 1e6 + + plt.loglog(mass, lifetime) + + plt.xlabel("Mass [Msun]") + plt.ylabel("Lifetime [Myr]") + + plt.show() diff --git a/examples/lifetime/run.sh b/examples/lifetime/run.sh new file mode 100644 index 0000000000000000000000000000000000000000..a180bae803170bc4e66db1432a2b367dffcdfe36 --- /dev/null +++ b/examples/lifetime/run.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +./plot_lifetime.py diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py index a30687558761c727d98cf4d34cc96708525cdd34..002511420fac8e158ac031bd79c41fbeea6bf6ef 100644 --- a/pyswiftsim/__init__.py +++ b/pyswiftsim/__init__.py @@ -97,6 +97,12 @@ class ParameterFile: "SupernovaeEnergy_erg": 1e51, "ThermalTime_Myr": 5, "ChemistryTable": "chemistry.hdf5", + }, + + "GEARLifetime": { + "quadratic": [-40.110, 5.509, 0.7824], + "linear": [141.929, -15.889, -3.2557], + "constant": [-261.365, 17.073, 9.8661], } } diff --git a/scripts/pyswift_initial_mass_function b/scripts/pyswift_initial_mass_function index 0e7a674597d03977d2546ebda4bc3bb1c5eee970..803e520495fc9d453963975de2e4100d20374bc5 100644 --- a/scripts/pyswift_initial_mass_function +++ b/scripts/pyswift_initial_mass_function @@ -87,8 +87,7 @@ if __name__ == "__main__": else: mass_max = np.log10(mass_limits[-1]) - # du / dt - print("Computing cooling...") + print("Computing initial mass function...") mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32) imf = libstellar.initialMassFunction(filename, mass) diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime new file mode 100644 index 0000000000000000000000000000000000000000..0b14c53f589d9a26c21841b4393536a000dc1eec --- /dev/null +++ b/scripts/pyswift_lifetime @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +# ######################################################################## +# This file is part of PYSWIFT. +# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# ######################################################################## + +from pyswiftsim import libstellar +import pyswiftsim + +import numpy as np +import matplotlib.pyplot as plt +from argparse import ArgumentParser +import os +plt.rc('text', usetex=True) + + +def parseArguments(): + parser = ArgumentParser(description="Plot the initial mass function") + + parser.add_argument( + "--download_table", dest="table", + action="store_true", + help="Download the default chemistry table") + + parser.add_argument( + "--params", dest="params", + type=str, default=None, + help="SWIFT parameters file") + + parser.add_argument( + "--N", dest="N", + type=int, default=1000, + help="Number of points") + + parser.add_argument( + "--mass_min", dest="mass_min", + type=float, default=None, + help="Minimal mass to plot") + + parser.add_argument( + "--mass_max", dest="mass_max", + type=float, default=None, + help="Maximal mass to plot") + + parser.add_argument( + "--metallicity", dest="metallicity", + type=float, default=0.02041, + help="Metallicity of the stars") + + args = parser.parse_args() + + return args + + +if __name__ == "__main__": + opt = parseArguments() + + if opt.table: + pyswiftsim.downloadChemistryTable() + + if opt.params is not None and not os.path.isfile(opt.params): + raise Exception( + "The parameter file '{}' does not exist".format(opt.params)) + + with pyswiftsim.ParameterFile(existing=opt.params) as filename: + + parser = pyswiftsim.parseYamlFile(filename) + imf_parser = parser["GEARInitialMassFunction"] + + mass_limits = list(imf_parser["mass_limits"]) + if opt.mass_min is not None: + mass_min = np.log10(opt.mass_min) + else: + mass_min = np.log10(mass_limits[0]) + + if opt.mass_max is not None: + mass_max = np.log10(opt.mass_max) + else: + mass_max = np.log10(mass_limits[-1]) + + mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32) + + # get metallicity + Z = opt.metallicity * np.ones(opt.N, dtype=np.float32) + + print("Computing lifetime...") + + lifetime = libstellar.lifetimeFromMass(filename, mass, Z) + lifetime /= 1e6 + + plt.loglog(mass, lifetime) + + plt.xlabel("Mass [Msun]") + plt.ylabel("Lifetime [Myr]") + + plt.show() diff --git a/tests/test_initial_mass_function.py b/tests/test_initial_mass_function.py index 152937f2c0293d0dfd54db8e4697d792c756e9eb..ca80c505bb39844035953934cdb84bf221591068 100644 --- a/tests/test_initial_mass_function.py +++ b/tests/test_initial_mass_function.py @@ -1,4 +1,21 @@ #!/usr/bin/env python3 +# ######################################################################## +# This file is part of PYSWIFT. +# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# ######################################################################## from pyswiftsim import libstellar import pyswiftsim @@ -16,9 +33,10 @@ if __name__ == "__main__": with pyswiftsim.ParameterFile() as filename: parser = pyswiftsim.parseYamlFile(filename) - # imf = parser["GEARInitialMassFunction"] - mass_min = float(1.) - mass_max = float(50.) + imf = parser["GEARInitialMassFunction"] + mass_limits = imf["mass_limits"] + mass_min = float(mass_limits[0]) + mass_max = float(mass_limits[1]) mass = np.linspace(mass_min, mass_max, N, dtype=np.float32) diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py new file mode 100644 index 0000000000000000000000000000000000000000..6fb9f092a99d17c0a937fec8a424ee8df2187c11 --- /dev/null +++ b/tests/test_lifetime.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +# ######################################################################## +# This file is part of PYSWIFT. +# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# ######################################################################## + +from pyswiftsim import libstellar +import pyswiftsim + +import numpy as np +import matplotlib.pyplot as plt +plt.rc('text', usetex=True) + + +N = 100 +solar_metallicity = 0.02041 + + +if __name__ == "__main__": + with pyswiftsim.ParameterFile() as filename: + + parser = pyswiftsim.parseYamlFile(filename) + # Get masses + imf = parser["GEARInitialMassFunction"] + mass_limits = imf["mass_limits"] + mass_min = float(mass_limits[0]) + mass_max = float(mass_limits[1]) + + mass = np.linspace(mass_min, mass_max, N, dtype=np.float32) + + # get metallicity + Z = solar_metallicity * np.ones(N, dtype=np.float32) + + print("Computing lifetime...") + + lifetime = libstellar.lifetimeFromMass(filename, mass, Z) + + print("Lifetime obtained: {}".format(lifetime))