Lifetime working

examples/lifetime/plot_lifetime.py

+#!/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()

examples/lifetime/run.sh

+#!/bin/bash
+
+./plot_lifetime.py

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],
         }
     }
 

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)

scripts/pyswift_lifetime

+#!/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()

tests/test_initial_mass_function.py

 #!/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)
 

tests/test_lifetime.py

+#!/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))