Skip to content
Snippets Groups Projects
Commit eefe0eda authored by Loic Hausammann's avatar Loic Hausammann
Browse files

All feedback parameters are read from tables

parent b91850cb
No related branches found
No related tags found
1 merge request!10Update SWIFT calls
...@@ -27,6 +27,10 @@ plt.rc('text', usetex=True) ...@@ -27,6 +27,10 @@ plt.rc('text', usetex=True)
N = 100 N = 100
solMass_in_cgs = 1.98848e33 solMass_in_cgs = 1.98848e33
mass_min = 0.05 # in solar mass
mass_max = 50 # in solar mass
pyswiftsim.downloadYieldsTable()
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
...@@ -34,14 +38,12 @@ if __name__ == "__main__": ...@@ -34,14 +38,12 @@ if __name__ == "__main__":
parser = pyswiftsim.parseYamlFile(filename) parser = pyswiftsim.parseYamlFile(filename)
# get parser # get parser
imf_parser = parser["GEARInitialMassFunction"]
units = parser["InternalUnitSystem"] units = parser["InternalUnitSystem"]
unit_mass = float(units["UnitMass_in_cgs"]) unit_mass = float(units["UnitMass_in_cgs"])
solMass_in_internal = solMass_in_cgs / unit_mass solMass_in_internal = solMass_in_cgs / unit_mass
mass_limits = list(imf_parser["mass_limits_msun"]) mass_min = np.log10(mass_min * solMass_in_internal)
mass_min = np.log10(mass_limits[0] * solMass_in_internal) mass_max = np.log10(mass_max * solMass_in_internal)
mass_max = np.log10(mass_limits[-1] * solMass_in_internal)
mass = np.logspace(mass_min, mass_max, N, dtype=np.float32) mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
...@@ -51,10 +53,10 @@ if __name__ == "__main__": ...@@ -51,10 +53,10 @@ if __name__ == "__main__":
imf *= solMass_in_internal imf *= solMass_in_internal
# plot IMF # plot IMF
plt.figure()
mass /= solMass_in_internal mass /= solMass_in_internal
plt.loglog(mass, imf) plt.loglog(mass, imf)
plt.figure()
date = strftime("%d %b %Y", gmtime()) date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date)) plt.title("Date: {}".format(date))
plt.xlabel("Mass [$M_\mathrm{sun}$]") plt.xlabel("Mass [$M_\mathrm{sun}$]")
......
...@@ -30,17 +30,18 @@ N = 1000 ...@@ -30,17 +30,18 @@ N = 1000
solar_metallicity = 0.02041 solar_metallicity = 0.02041
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
year_in_cgs = 3.15569252e7 year_in_cgs = 3.15569252e7
mass_min = 0.05
mass_max = 50
pyswiftsim.downloadYieldsTable()
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
parser = pyswiftsim.parseYamlFile(filename) parser = pyswiftsim.parseYamlFile(filename)
# Get masses # Get masses
imf = parser["GEARInitialMassFunction"] mass_min = np.log10(mass_min)
mass_limits = imf["mass_limits_msun"] mass_max = np.log10(mass_max)
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) mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
......
...@@ -33,6 +33,7 @@ t_max = 3e10 ...@@ -33,6 +33,7 @@ t_max = 3e10
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
year_in_cgs = 3.15569252e7 year_in_cgs = 3.15569252e7
pyswiftsim.downloadYieldsTable()
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
......
...@@ -22,6 +22,32 @@ from tempfile import NamedTemporaryFile ...@@ -22,6 +22,32 @@ from tempfile import NamedTemporaryFile
import yaml 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(): def downloadCoolingTable():
""" """
Download the cooling table if it does not already exist in the current Download the cooling table if it does not already exist in the current
...@@ -30,10 +56,7 @@ def downloadCoolingTable(): ...@@ -30,10 +56,7 @@ def downloadCoolingTable():
url = "http://virgodb.cosma.dur.ac.uk/" url = "http://virgodb.cosma.dur.ac.uk/"
url += "swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5" url += "swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5"
filename = "CloudyData_UVB=HM2012.h5" filename = "CloudyData_UVB=HM2012.h5"
download(url, filename)
if not os.path.isfile(filename):
# Download the file from `url` and save it locally under `file_name`:
urllib.request.urlretrieve(url, filename)
def downloadYieldsTable(): def downloadYieldsTable():
...@@ -41,13 +64,10 @@ def downloadYieldsTable(): ...@@ -41,13 +64,10 @@ def downloadYieldsTable():
Download the yields table if it does not already exist in the current Download the yields table if it does not already exist in the current
repository. repository.
""" """
url = "https://obswww.unige.ch/" url = "https://obswww.unige.ch/lastro/projects/"
url += "lastro/projects/Clastro/PySWIFTsim/" url += "Clastro/PySWIFTsim/chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5"
filename = "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5" filename = "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5"
download(url, filename)
if not os.path.isfile(filename):
# Download the file from `url` and save it locally under `file_name`:
urllib.request.urlretrieve(url, filename)
def parseYamlFile(filename): def parseYamlFile(filename):
...@@ -155,40 +175,11 @@ class ParameterFile: ...@@ -155,40 +175,11 @@ class ParameterFile:
"InitialMetallicity": 0.01295, "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": { "GEARFeedback": {
"SupernovaeEnergy_erg": 1e51, "SupernovaeEnergy_erg": 1e51,
"ThermalTime_Myr": 5, "ThermalTime_Myr": 5,
"ChemistryTable": "chemistry.hdf5", "YieldsTable": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
},
"GEARLifetime": {
"quadratic": [-40.110, 5.509, 0.7824],
"linear": [141.929, -15.889, -3.2557],
"constant": [-261.365, 17.073, 9.8661],
}, },
"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 Dictionnary containing all the default parameters to write in the
......
...@@ -49,12 +49,12 @@ def parseArguments(): ...@@ -49,12 +49,12 @@ def parseArguments():
parser.add_argument( parser.add_argument(
"--mass_min", dest="mass_min", "--mass_min", dest="mass_min",
type=float, default=None, type=float, default=0.05,
help="Minimal mass to plot") help="Minimal mass to plot")
parser.add_argument( parser.add_argument(
"--mass_max", dest="mass_max", "--mass_max", dest="mass_max",
type=float, default=None, type=float, default=50,
help="Maximal mass to plot") help="Maximal mass to plot")
args = parser.parse_args() args = parser.parse_args()
...@@ -66,7 +66,7 @@ if __name__ == "__main__": ...@@ -66,7 +66,7 @@ if __name__ == "__main__":
opt = parseArguments() opt = parseArguments()
if opt.table: if opt.table:
pyswiftsim.downloadChemistryTable() pyswiftsim.downloadYieldsTable()
if opt.params is not None and not os.path.isfile(opt.params): if opt.params is not None and not os.path.isfile(opt.params):
raise Exception( raise Exception(
...@@ -75,27 +75,17 @@ if __name__ == "__main__": ...@@ -75,27 +75,17 @@ if __name__ == "__main__":
with pyswiftsim.ParameterFile(existing=opt.params) as filename: with pyswiftsim.ParameterFile(existing=opt.params) as filename:
parser = pyswiftsim.parseYamlFile(filename) parser = pyswiftsim.parseYamlFile(filename)
imf_parser = parser["GEARInitialMassFunction"]
units = parser["InternalUnitSystem"] units = parser["InternalUnitSystem"]
unit_mass = float(units["UnitMass_in_cgs"]) unit_mass = float(units["UnitMass_in_cgs"])
mass_limits = list(imf_parser["mass_limits_msun"]) mass_min = np.log10(opt.mass_min)
if opt.mass_min is not None: mass_max = np.log10(opt.mass_max)
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
print("Computing initial mass function...") print("Computing initial mass function...")
mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32) mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
mass *= solMass_in_cgs / unit_mass
imf = libstellar.initialMassFunction(filename, mass) imf = libstellar.initialMassFunction(filename, mass)
imf *= solMass_in_cgs / unit_mass imf *= solMass_in_cgs / unit_mass
......
...@@ -29,6 +29,9 @@ plt.rc('text', usetex=True) ...@@ -29,6 +29,9 @@ plt.rc('text', usetex=True)
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
year_in_cgs = 3.15569252e7 year_in_cgs = 3.15569252e7
mass_min = 0.05
mass_max = 50
def parseArguments(): def parseArguments():
parser = ArgumentParser(description="Plot the lifetime of a star") parser = ArgumentParser(description="Plot the lifetime of a star")
...@@ -50,12 +53,12 @@ def parseArguments(): ...@@ -50,12 +53,12 @@ def parseArguments():
parser.add_argument( parser.add_argument(
"--mass_min", dest="mass_min", "--mass_min", dest="mass_min",
type=float, default=None, type=float, default=0.05,
help="Minimal mass to plot") help="Minimal mass to plot")
parser.add_argument( parser.add_argument(
"--mass_max", dest="mass_max", "--mass_max", dest="mass_max",
type=float, default=None, type=float, default=50,
help="Maximal mass to plot") help="Maximal mass to plot")
parser.add_argument( parser.add_argument(
...@@ -72,7 +75,7 @@ if __name__ == "__main__": ...@@ -72,7 +75,7 @@ if __name__ == "__main__":
opt = parseArguments() opt = parseArguments()
if opt.table: if opt.table:
pyswiftsim.downloadChemistryTable() pyswiftsim.downloadYieldsTable()
if opt.params is not None and not os.path.isfile(opt.params): if opt.params is not None and not os.path.isfile(opt.params):
raise Exception( raise Exception(
...@@ -81,18 +84,9 @@ if __name__ == "__main__": ...@@ -81,18 +84,9 @@ if __name__ == "__main__":
with pyswiftsim.ParameterFile(existing=opt.params) as filename: with pyswiftsim.ParameterFile(existing=opt.params) as filename:
parser = pyswiftsim.parseYamlFile(filename) parser = pyswiftsim.parseYamlFile(filename)
imf_parser = parser["GEARInitialMassFunction"]
mass_min = np.log10(opt.mass_min)
mass_limits = list(imf_parser["mass_limits_msun"]) mass_max = np.log10(opt.mass_max)
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) mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
......
...@@ -82,7 +82,7 @@ if __name__ == "__main__": ...@@ -82,7 +82,7 @@ if __name__ == "__main__":
opt = parseArguments() opt = parseArguments()
if opt.table: if opt.table:
pyswiftsim.downloadChemistryTable() pyswiftsim.downloadYieldsTable()
if opt.params is not None and not os.path.isfile(opt.params): if opt.params is not None and not os.path.isfile(opt.params):
raise Exception( raise Exception(
......
...@@ -27,17 +27,16 @@ plt.rc('text', usetex=True) ...@@ -27,17 +27,16 @@ plt.rc('text', usetex=True)
N = 100 N = 100
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
mass_min = 0.5
mass_max = 45
pyswiftsim.downloadYieldsTable()
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
parser = pyswiftsim.parseYamlFile(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) mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
......
...@@ -29,16 +29,13 @@ N = 100 ...@@ -29,16 +29,13 @@ N = 100
solar_metallicity = 0.02041 solar_metallicity = 0.02041
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
year_in_cgs = 3.15569252e7 year_in_cgs = 3.15569252e7
mass_min = 0.05 # in solar mass
mass_max = 50
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
parser = pyswiftsim.parseYamlFile(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) mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
......
...@@ -32,6 +32,8 @@ t_max = 1e10 ...@@ -32,6 +32,8 @@ t_max = 1e10
solMass_in_cgs = 1.989e33 solMass_in_cgs = 1.989e33
year_in_cgs = 3.15569252e7 year_in_cgs = 3.15569252e7
pyswiftsim.downloadYieldsTable()
if __name__ == "__main__": if __name__ == "__main__":
with pyswiftsim.ParameterFile() as filename: with pyswiftsim.ParameterFile() as filename:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment