From f1be5083b2937a50fec3b66b64f97d2f9dda2d11 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 3 Jun 2019 17:09:49 +0200
Subject: [PATCH] Lifetime working
---
examples/lifetime/plot_lifetime.py | 57 ++++++++++++++
examples/lifetime/run.sh | 3 +
pyswiftsim/__init__.py | 6 ++
scripts/pyswift_initial_mass_function | 3 +-
scripts/pyswift_lifetime | 109 ++++++++++++++++++++++++++
tests/test_initial_mass_function.py | 24 +++++-
tests/test_lifetime.py | 51 ++++++++++++
7 files changed, 248 insertions(+), 5 deletions(-)
create mode 100644 examples/lifetime/plot_lifetime.py
create mode 100644 examples/lifetime/run.sh
create mode 100644 scripts/pyswift_lifetime
create mode 100644 tests/test_lifetime.py
diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py
new file mode 100644
index 0000000..a1e1e9f
--- /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 0000000..a180bae
--- /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 a306875..0025114 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 0e7a674..803e520 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 0000000..0b14c53
--- /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 152937f..ca80c50 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 0000000..6fb9f09
--- /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))
--
GitLab