From eefe0eda02404bc426a0459e207d935d5288750d Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Fri, 7 Jun 2019 15:12:41 +0200
Subject: [PATCH] All feedback parameters are read from tables

---
 .../plot_initial_mass_function.py             | 12 ++--
 examples/lifetime/plot_lifetime.py            |  9 +--
 examples/supernovae_rate/plot_rates.py        |  1 +
 pyswiftsim/__init__.py                        | 71 ++++++++-----------
 scripts/pyswift_initial_mass_function         | 24 ++-----
 scripts/pyswift_lifetime                      | 24 +++----
 scripts/pyswift_supernovae_rate               |  2 +-
 tests/test_initial_mass_function.py           |  7 +-
 tests/test_lifetime.py                        |  7 +-
 tests/test_rates.py                           |  2 +
 10 files changed, 68 insertions(+), 91 deletions(-)

diff --git a/examples/initial_mass_function/plot_initial_mass_function.py b/examples/initial_mass_function/plot_initial_mass_function.py
index bb97239..6957770 100644
--- a/examples/initial_mass_function/plot_initial_mass_function.py
+++ b/examples/initial_mass_function/plot_initial_mass_function.py
@@ -27,6 +27,10 @@ plt.rc('text', usetex=True)
 
 N = 100
 solMass_in_cgs = 1.98848e33
+mass_min = 0.05  # in solar mass
+mass_max = 50  # in solar mass
+
+pyswiftsim.downloadYieldsTable()
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
@@ -34,14 +38,12 @@ if __name__ == "__main__":
         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_msun"])
-        mass_min = np.log10(mass_limits[0] * solMass_in_internal)
-        mass_max = np.log10(mass_limits[-1] * solMass_in_internal)
+        mass_min = np.log10(mass_min * solMass_in_internal)
+        mass_max = np.log10(mass_max * solMass_in_internal)
 
         mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
 
@@ -51,10 +53,10 @@ if __name__ == "__main__":
         imf *= solMass_in_internal
 
         # plot IMF
+        plt.figure()
         mass /= solMass_in_internal
         plt.loglog(mass, imf)
 
-        plt.figure()
         date = strftime("%d %b %Y", gmtime())
         plt.title("Date: {}".format(date))
         plt.xlabel("Mass [$M_\mathrm{sun}$]")
diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py
index 79998ef..3cfcdec 100644
--- a/examples/lifetime/plot_lifetime.py
+++ b/examples/lifetime/plot_lifetime.py
@@ -30,17 +30,18 @@ N = 1000
 solar_metallicity = 0.02041
 solMass_in_cgs = 1.989e33
 year_in_cgs = 3.15569252e7
+mass_min = 0.05
+mass_max = 50
 
+pyswiftsim.downloadYieldsTable()
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
 
         parser = pyswiftsim.parseYamlFile(filename)
         # Get masses
-        imf = parser["GEARInitialMassFunction"]
-        mass_limits = imf["mass_limits_msun"]
-        mass_min = np.log10(mass_limits[0])
-        mass_max = np.log10(mass_limits[-1])
+        mass_min = np.log10(mass_min)
+        mass_max = np.log10(mass_max)
 
         mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
 
diff --git a/examples/supernovae_rate/plot_rates.py b/examples/supernovae_rate/plot_rates.py
index 300c5fc..ac40330 100644
--- a/examples/supernovae_rate/plot_rates.py
+++ b/examples/supernovae_rate/plot_rates.py
@@ -33,6 +33,7 @@ t_max = 3e10
 solMass_in_cgs = 1.989e33
 year_in_cgs = 3.15569252e7
 
+pyswiftsim.downloadYieldsTable()
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 6649bc3..dfdc68f 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -22,6 +22,32 @@ from tempfile import NamedTemporaryFile
 import yaml
 
 
+def download(url, filename):
+    """
+    Download a file from a given url.
+
+    Parameters
+    ==========
+
+    url: str
+        The url to download.
+
+    filename: str
+        The filename to use for saving the file.
+    """
+    progress = [0, 0]
+
+    def report(count, size, total):
+        progress[0] = 100. * count * size / total
+        if progress[0] - progress[1] > 10:
+            progress[1] = progress[0]
+            print("Downloaded {}% ...".format(int(progress[1])))
+
+    if not os.path.isfile(filename):
+        # Download the file from `url` and save it locally under `file_name`:
+        urllib.request.urlretrieve(url, filename, reporthook=report)
+
+
 def downloadCoolingTable():
     """
     Download the cooling table if it does not already exist in the current
@@ -30,10 +56,7 @@ def downloadCoolingTable():
     url = "http://virgodb.cosma.dur.ac.uk/"
     url += "swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5"
     filename = "CloudyData_UVB=HM2012.h5"
-
-    if not os.path.isfile(filename):
-        # Download the file from `url` and save it locally under `file_name`:
-        urllib.request.urlretrieve(url, filename)
+    download(url, filename)
 
 
 def downloadYieldsTable():
@@ -41,13 +64,10 @@ def downloadYieldsTable():
     Download the yields table if it does not already exist in the current
     repository.
     """
-    url = "https://obswww.unige.ch/"
-    url += "lastro/projects/Clastro/PySWIFTsim/"
+    url = "https://obswww.unige.ch/lastro/projects/"
+    url += "Clastro/PySWIFTsim/chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5"
     filename = "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5"
-
-    if not os.path.isfile(filename):
-        # Download the file from `url` and save it locally under `file_name`:
-        urllib.request.urlretrieve(url, filename)
+    download(url, filename)
 
 
 def parseYamlFile(filename):
@@ -155,40 +175,11 @@ class ParameterFile:
             "InitialMetallicity": 0.01295,
         },
 
-        "GEARInitialMassFunction": {
-            "number_function_part": 4,
-            "exponents": [0.7, -0.8, -1.7, -1.3],
-            "mass_limits_msun": [0.05, 0.08, 0.5, 1.0, 50.]
-        },
-
         "GEARFeedback": {
             "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],
+            "YieldsTable": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
         },
-
-        "GEARSupernovaeIa": {
-            "exponent": -0.35,
-            "min_mass_white_dwarf_progenitor": 3.,
-            "max_mass_white_dwarf_progenitor": 8.,
-            "min_mass_red_giant": 0.9,
-            "max_mass_red_giant": 1.5,
-            "coef_red_giant": 0.02,
-            "min_mass_main_sequence": 1.8,
-            "max_mass_main_sequence": 2.5,
-            "coef_main_sequence": 0.05,
-        },
-
-        "GEARSupernovaeII": {
-            "min_mass": 8.,
-            "max_mass": 50.,
-        }
     }
     """
     Dictionnary containing all the default parameters to write in the
diff --git a/scripts/pyswift_initial_mass_function b/scripts/pyswift_initial_mass_function
index fa35ce5..d62cad5 100644
--- a/scripts/pyswift_initial_mass_function
+++ b/scripts/pyswift_initial_mass_function
@@ -49,12 +49,12 @@ def parseArguments():
 
     parser.add_argument(
         "--mass_min", dest="mass_min",
-        type=float, default=None,
+        type=float, default=0.05,
         help="Minimal mass to plot")
 
     parser.add_argument(
         "--mass_max", dest="mass_max",
-        type=float, default=None,
+        type=float, default=50,
         help="Maximal mass to plot")
 
     args = parser.parse_args()
@@ -66,7 +66,7 @@ if __name__ == "__main__":
     opt = parseArguments()
 
     if opt.table:
-        pyswiftsim.downloadChemistryTable()
+        pyswiftsim.downloadYieldsTable()
 
     if opt.params is not None and not os.path.isfile(opt.params):
         raise Exception(
@@ -75,27 +75,17 @@ if __name__ == "__main__":
     with pyswiftsim.ParameterFile(existing=opt.params) as filename:
 
         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_msun"])
-        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_min *= solMass_in_cgs / unit_mass
-        mass_max *= solMass_in_cgs / unit_mass
+        mass_min = np.log10(opt.mass_min)
+        mass_max = np.log10(opt.mass_max)
 
         print("Computing initial mass function...")
 
         mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
+        mass *= solMass_in_cgs / unit_mass
+
         imf = libstellar.initialMassFunction(filename, mass)
 
         imf *= solMass_in_cgs / unit_mass
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
index 9fb7754..b576e86 100644
--- a/scripts/pyswift_lifetime
+++ b/scripts/pyswift_lifetime
@@ -29,6 +29,9 @@ plt.rc('text', usetex=True)
 solMass_in_cgs = 1.989e33
 year_in_cgs = 3.15569252e7
 
+mass_min = 0.05
+mass_max = 50
+
 
 def parseArguments():
     parser = ArgumentParser(description="Plot the lifetime of a star")
@@ -50,12 +53,12 @@ def parseArguments():
 
     parser.add_argument(
         "--mass_min", dest="mass_min",
-        type=float, default=None,
+        type=float, default=0.05,
         help="Minimal mass to plot")
 
     parser.add_argument(
         "--mass_max", dest="mass_max",
-        type=float, default=None,
+        type=float, default=50,
         help="Maximal mass to plot")
 
     parser.add_argument(
@@ -72,7 +75,7 @@ if __name__ == "__main__":
     opt = parseArguments()
 
     if opt.table:
-        pyswiftsim.downloadChemistryTable()
+        pyswiftsim.downloadYieldsTable()
 
     if opt.params is not None and not os.path.isfile(opt.params):
         raise Exception(
@@ -81,18 +84,9 @@ if __name__ == "__main__":
     with pyswiftsim.ParameterFile(existing=opt.params) as filename:
 
         parser = pyswiftsim.parseYamlFile(filename)
-        imf_parser = parser["GEARInitialMassFunction"]
-
-        mass_limits = list(imf_parser["mass_limits_msun"])
-        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_min = np.log10(opt.mass_min)
+        mass_max = np.log10(opt.mass_max)
 
         mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
 
diff --git a/scripts/pyswift_supernovae_rate b/scripts/pyswift_supernovae_rate
index a60baf0..93e461f 100644
--- a/scripts/pyswift_supernovae_rate
+++ b/scripts/pyswift_supernovae_rate
@@ -82,7 +82,7 @@ if __name__ == "__main__":
     opt = parseArguments()
 
     if opt.table:
-        pyswiftsim.downloadChemistryTable()
+        pyswiftsim.downloadYieldsTable()
 
     if opt.params is not None and not os.path.isfile(opt.params):
         raise Exception(
diff --git a/tests/test_initial_mass_function.py b/tests/test_initial_mass_function.py
index 416aa3c..d77b0e2 100644
--- a/tests/test_initial_mass_function.py
+++ b/tests/test_initial_mass_function.py
@@ -27,17 +27,16 @@ plt.rc('text', usetex=True)
 
 N = 100
 solMass_in_cgs = 1.989e33
+mass_min = 0.5
+mass_max = 45
 
+pyswiftsim.downloadYieldsTable()
 
 if __name__ == "__main__":
 
     with pyswiftsim.ParameterFile() as filename:
 
         parser = pyswiftsim.parseYamlFile(filename)
-        imf = parser["GEARInitialMassFunction"]
-        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)
 
diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py
index aac4281..cc1f75e 100644
--- a/tests/test_lifetime.py
+++ b/tests/test_lifetime.py
@@ -29,16 +29,13 @@ N = 100
 solar_metallicity = 0.02041
 solMass_in_cgs = 1.989e33
 year_in_cgs = 3.15569252e7
+mass_min = 0.05  # in solar mass
+mass_max = 50
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
 
         parser = pyswiftsim.parseYamlFile(filename)
-        # Get masses
-        imf = parser["GEARInitialMassFunction"]
-        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)
 
diff --git a/tests/test_rates.py b/tests/test_rates.py
index 7bc38ff..6d879d6 100644
--- a/tests/test_rates.py
+++ b/tests/test_rates.py
@@ -32,6 +32,8 @@ t_max = 1e10
 solMass_in_cgs = 1.989e33
 year_in_cgs = 3.15569252e7
 
+pyswiftsim.downloadYieldsTable()
+
 
 if __name__ == "__main__":
     with pyswiftsim.ParameterFile() as filename:
-- 
GitLab