diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index fbaeba529ff1208399d417e7e00419057632cc9b..542d9df5512176037d5602266d177978c37f59d6 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -3,7 +3,8 @@ image: ubuntu:latest
before_script:
- - apt-get update -qq && apt-get install -y -qq build-essential libhdf5-serial-dev gfortran csh git libtool-bin autoconf
+ - export DEBIAN_FRONTEND=noninteractive
+ - apt-get update -qq && apt-get install -y -qq build-essential libhdf5-serial-dev gfortran csh git libtool-bin autoconf openmpi-bin openmpi-common libopenmpi-dev python3 python3-setuptools python3-numpy python3-scipy python3-matplotlib python3-dev python3-h5py
- export LHOME=`pwd`
# Compile grackle
- mkdir grackle
@@ -27,7 +28,7 @@ before_script:
- cd ..
# compile pyswiftsim
- - python setup.py install --user
+ - python3 setup.py install --with-swift $LHOME/swiftsim --with-grackle $GRACKLE_ROOT --user
# compile the docs
- cd docs
diff --git a/docs/RST/examples/stellar.rst b/docs/RST/examples/stellar.rst
index 94e1ac8ee917b03ae0ee1fd50d3886c07f312707..031a8849b978958ec5b183d29be70ae421c16cd0 100644
--- a/docs/RST/examples/stellar.rst
+++ b/docs/RST/examples/stellar.rst
@@ -21,5 +21,14 @@ The lifetime in GEAR is a quadratic in metallicity and mass (see equation 3.32 i
.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/lifetime.png
+Supernovae Rate
+~~~~~~~~~~~~~~~
+
+This example is generated in ``examples/supernovae_rate`` with GEAR's stellar model.
+The supernovae rate is given by the contribution of SNIa (in blue) and SNII (in orange).
+The two peaks in SNIa correspond to the two different type of companion (red giant and main sequence).
+
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/supernovae_rate.png
+
.. [Kroupa2001] `On the variation of the Initial Mass Function <https://arxiv.org/abs/astro-ph/0009005>`_
.. [Poirier2004] `Étude de l'évolution chimique et dynamique d'objets proto-galactiques: Application à l'évolution des galaxies spirales <https://obswww.unige.ch/~revaz/PyChem/_downloads/ff03ba21204a3f36fa4ec885abd545ae/ThesePoirier.pdf>`_
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
index 18bbb26903a25e87667887dc5beeacc8f3feb976..0a3f2ee4c5262c4cdeb3dc22058ae7953af98bd1 100644
--- a/examples/cooling_box/plot_cooling.py
+++ b/examples/cooling_box/plot_cooling.py
@@ -131,7 +131,7 @@ if __name__ == "__main__":
overwrite = {
"GrackleCooling": {
- "WithMetalCooling": with_metals,
+ "with_metal_cooling": with_metals,
}
}
with pyswiftsim.ParameterFile(overwrite) as filename:
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 884f0592da44acc0444e14e41ea63dcba67bf4eb..0fe4d8b42a5c7f99645b9957120f738a36817e32 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -20,6 +20,8 @@
from pyswiftsim.libcooling import coolingRate
import pyswiftsim
+import matplotlib
+
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
@@ -27,8 +29,20 @@ from astropy import units, constants
from time import strftime, gmtime
import sys
+SMALL_SIZE = 13
+MEDIUM_SIZE = 15
+BIGGER_SIZE = 18
+
+plt.rc('font', size=SMALL_SIZE) # controls default text sizes
+plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
+plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
plt.rc('text', usetex=True)
+
pyswiftsim.downloadCoolingTable()
# PARAMETERS
@@ -40,6 +54,12 @@ with_cosmo = False
# include metals?
with_metals = 1
+metallicity = 0.02041
+
+# cooling table
+cooling_tables = [
+ "./CloudyData_UVB=HM2012.h5",
+]
# hydrogen mass fraction
h_mass_frac = 0.76
@@ -50,6 +70,7 @@ if N_rho == 1:
density = np.array([1.])
else:
density = np.logspace(-6, 4, N_rho)
+ plot_cooling_length = False
# temperature in K
N_T = 1000
@@ -113,7 +134,7 @@ def generate_initial_condition(us):
return d
-def plot_solution(rate, data, us):
+def plot_solution(rates, data, us):
print("Plotting solution")
energy = data["energy"]
rho = data["density"]
@@ -131,55 +152,82 @@ def plot_solution(rate, data, us):
energy *= u_v**2
- rate *= u_v**2 / u_t
+ rates = [rate * u_v**2 / u_t for rate in rates]
# lambda cooling
- lam = rate * rho
+ lambdas = [rate * rho for rate in rates]
# do plot
if N_rho == 1:
# plot Lambda vs T
plt.figure()
- plt.loglog(T, np.abs(lam),
- label="SWIFT")
+ for i, lam in enumerate(lambdas):
+ plt.loglog(T, np.abs(lam),
+ label="Table %i" % i)
plt.xlabel("Temperature [K]")
plt.ylabel("$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]")
+ plt.legend()
else:
m_p = constants.m_p.to("g") / units.g
rho /= m_p
-
- shape = [N_rho, N_T]
- cooling_time = energy / rate
- cooling_length = np.sqrt(gamma * (gamma-1.) * energy) * cooling_time
-
- cooling_length = np.log10(np.abs(cooling_length) / units.kpc.to('cm'))
-
- # reshape
- rho = rho.reshape(shape)
- T = T.reshape(shape)
- energy = energy.reshape(shape)
- cooling_length = cooling_length.reshape(shape)
-
- _min = -7
- _max = 7
- N_levels = 100
- levels = np.linspace(_min, _max, N_levels)
- plt.figure()
- plt.contourf(rho, T, cooling_length, levels,
- cmap=plt.get_cmap(colormap))
- plt.xlabel("Density [atom/cm3]")
- plt.ylabel("Temperature [K]")
-
- ax = plt.gca()
- ax.set_xscale("log")
- ax.set_yscale("log")
-
- cbar = plt.colorbar()
- tc = np.arange(_min, _max, 1.5)
- cbar.set_ticks(tc)
- cbar.set_ticklabels(tc)
- cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
+ r = []
+ for rate in rates:
+ shape = [N_rho, N_T]
+ energy = np.reshape(energy, shape)
+ rate = np.reshape(rate, shape)
+
+ cooling_time = energy / rate
+ cooling_length = np.sqrt(gamma * (gamma-1.) * energy)
+ cooling_length *= cooling_time
+
+ cooling_length = np.log10(
+ np.abs(cooling_length) / units.kpc.to('cm'))
+
+ # reshape
+ rho = rho.reshape(shape)
+ T = T.reshape(shape)
+ energy = energy.reshape(shape)
+ cooling_length = cooling_length.reshape(shape)
+
+ levels = 500
+ if plot_cooling_length:
+ _min = -7
+ _max = 7
+ levels = np.linspace(_min, _max, levels)
+ z = cooling_length
+ else:
+ z = np.log10(np.abs(rate))
+
+ plt.figure()
+ r.append(rate)
+ plt.contourf(np.log10(rho), np.log10(T), z,
+ levels,
+ cmap=plt.get_cmap(colormap))
+ plt.xlabel("Log( Density [atom/cm3] )")
+ plt.ylabel("Log( Temperature [K] )")
+
+ cbar = plt.colorbar()
+ if plot_cooling_length:
+ tc = np.arange(_min, _max, 1.5)
+ cbar.set_ticks(tc)
+ cbar.set_ticklabels(tc)
+ cbar.ax.set_ylabel("Log( Cooling Length [kpc] )")
+ else:
+ cbar.ax.set_ylabel("Log( Rate [s / erg] )")
+
+ if len(r) == 2:
+ plt.figure()
+ plt.contourf(rho, T, np.log10(np.abs((r[1] - r[0]) / r[1])), 500,
+ cmap=plt.get_cmap(colormap))
+ plt.xlabel("Density [atom/cm3]")
+ plt.ylabel("Temperature [K]")
+
+ ax = plt.gca()
+ ax.set_xscale("log")
+ ax.set_yscale("log")
+
+ cbar = plt.colorbar()
date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date))
@@ -188,24 +236,35 @@ def plot_solution(rate, data, us):
if __name__ == "__main__":
- overwrite = {
- "GrackleCooling": {
- "WithMetalCooling": with_metals,
+ rates = []
+
+ for cooling_table in cooling_tables:
+ overwrite = {
+ "GrackleCooling": {
+ "with_metal_cooling": with_metals,
+ "redshift": 0,
+ "cloudy_table": cooling_table,
+ },
+ "GEARChemistry": {
+ "initial_metallicity": metallicity,
+ },
}
- }
- with pyswiftsim.ParameterFile(overwrite) as filename:
- parser = pyswiftsim.parseYamlFile(filename)
- us = parser["InternalUnitSystem"]
+ with pyswiftsim.ParameterFile(overwrite) as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ us = parser["InternalUnitSystem"]
+
+ d = generate_initial_condition(us)
- d = generate_initial_condition(us)
+ # du / dt
+ print("Computing cooling...")
- # du / dt
- print("Computing cooling...")
+ rate = coolingRate(filename,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ with_cosmo, d["time_step"])
- rate = coolingRate(filename,
- d["density"].astype(np.float32),
- d["energy"].astype(np.float32),
- with_cosmo, d["time_step"])
+ rates.append(rate)
- plot_solution(rate, d, us)
+ plot_solution(rates, d, us)
diff --git a/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
index ce3d510430474c3986b463878c3edff2ba02d212..c140dca3d370df398e77fa31d871e97e9652c7a6 100644
--- a/examples/cooling_rate/run.sh
+++ b/examples/cooling_rate/run.sh
@@ -1,4 +1,4 @@
#!/bin/bash
./plot_cooling.py 1
-./plot_cooling.py 100
+# ./plot_cooling.py 100
diff --git a/examples/initial_mass_function/plot_initial_mass_function.py b/examples/initial_mass_function/plot_initial_mass_function.py
index bc653edc70a36756b7f65f206fd4e5e02cd56129..873aa4bb6bedaa64f0971476bcbbd2ff5e9fd681 100644
--- a/examples/initial_mass_function/plot_initial_mass_function.py
+++ b/examples/initial_mass_function/plot_initial_mass_function.py
@@ -22,32 +22,45 @@ import pyswiftsim
import numpy as np
import matplotlib.pyplot as plt
+from time import strftime, gmtime
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:
parser = pyswiftsim.parseYamlFile(filename)
- imf_parser = parser["GEARInitialMassFunction"]
- mass_limits = list(imf_parser["mass_limits"])
- mass_min = np.log10(mass_limits[0])
- mass_max = np.log10(mass_limits[-1])
+ # get parser
+ units = parser["InternalUnitSystem"]
+ unit_mass = float(units["UnitMass_in_cgs"])
+ solMass_in_internal = solMass_in_cgs / unit_mass
+
+ 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,
+ endpoint=False)
- # du / dt
print("Computing initial mass function...")
- mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
imf = libstellar.initialMassFunction(filename, mass)
+ imf *= solMass_in_internal
# plot IMF
+ plt.figure()
+ mass /= solMass_in_internal
plt.loglog(mass, imf)
- plt.xlabel("Mass [Msun]")
- plt.ylabel("Mass fraction")
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
+ plt.xlabel("Mass [$M_\mathrm{sun}$]")
+ plt.ylabel("Initial Mass Function [$M_\mathrm{sun}^{-1}$]")
plt.show()
diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py
index 7696311e3cfd048d2dd000bebbbfc4389d448125..3cfcdec817a030b4c977a2f7760e2012db2dd854 100644
--- a/examples/lifetime/plot_lifetime.py
+++ b/examples/lifetime/plot_lifetime.py
@@ -22,25 +22,38 @@ import pyswiftsim
import numpy as np
import matplotlib.pyplot as plt
+from time import strftime, gmtime
plt.rc('text', usetex=True)
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"]
- 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)
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+
+ mass *= solMass_in_cgs / unit_mass
+
# get metallicity
Z = solar_metallicity * np.ones(N, dtype=np.float32)
@@ -49,8 +62,13 @@ if __name__ == "__main__":
lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
- lifetime /= 1e6
+ mass /= solMass_in_cgs / unit_mass
+ mass_from_lifetime /= solMass_in_cgs / unit_mass
+ lifetime /= 1e6 * year_in_cgs / unit_time
+ plt.figure()
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
plt.loglog(mass, lifetime, label="Lifetime from mass")
plt.loglog(mass_from_lifetime, lifetime, "--",
label="Mass from lifetime")
diff --git a/examples/supernovae_rate/plot_rates.py b/examples/supernovae_rate/plot_rates.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac40330d786144d1b0c0379ed9f779b3c3f1f5a7
--- /dev/null
+++ b/examples/supernovae_rate/plot_rates.py
@@ -0,0 +1,84 @@
+#!/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 time import strftime, gmtime
+plt.rc('text', usetex=True)
+
+
+N = 1000
+solar_metallicity = 0.02041
+t_min = 1e6 # in yr
+t_max = 3e10
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
+
+pyswiftsim.downloadYieldsTable()
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+ parser = pyswiftsim.parseYamlFile(filename)
+
+ # Get times
+ times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+
+ times *= year_in_cgs / unit_time
+
+ t1 = times[:-1]
+ t2 = times[1:]
+
+ # get metallicity
+ Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+ print("Computing supernovae rate...")
+
+ rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+ rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+ t1 /= 1e6 * year_in_cgs / unit_time
+
+ unit_rate = year_in_cgs / unit_time
+ unit_rate *= solMass_in_cgs / unit_mass
+
+ rate_snia *= 1e9 * unit_rate
+ rate_snii *= 1e9 * unit_rate
+
+ plt.figure()
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
+ plt.loglog(t1, rate_snia, label="SNIa")
+ plt.loglog(t1, rate_snii, label="SNII")
+
+ plt.xlabel("Time [Myr]")
+ plt.ylabel("Supernovae rate [$M_\odot^{-1} \mathrm{Gyr}^{-1}$]")
+ plt.legend()
+
+ plt.show()
diff --git a/examples/supernovae_rate/run.sh b/examples/supernovae_rate/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..4264cedc3634a1af6e88d23fc110956587a0932b
--- /dev/null
+++ b/examples/supernovae_rate/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python plot_rates.py
diff --git a/examples/supernovae_yields/plot_yields.py b/examples/supernovae_yields/plot_yields.py
new file mode 100644
index 0000000000000000000000000000000000000000..d8f1c94eacdf95735db86ee098bdf86f53093e70
--- /dev/null
+++ b/examples/supernovae_yields/plot_yields.py
@@ -0,0 +1,90 @@
+#!/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 time import strftime, gmtime
+plt.rc('text', usetex=True)
+
+
+N = 1000
+m_min = 0.08 # in solar mass
+m_max = 50
+m_wd = 1.4 # white dwarf mass
+solMass_in_cgs = 1.989e33
+
+pyswiftsim.downloadYieldsTable()
+
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+ parser = pyswiftsim.parseYamlFile(filename)
+
+ # Get masses
+ mass = np.logspace(np.log10(m_min), np.log10(m_max), N,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+
+ mass *= solMass_in_cgs / unit_mass
+
+ m1 = np.ones(N, dtype=np.float32) * 1e-20 # avoid log(0)
+
+ print("Computing supernovae yields...")
+
+ yields_snii = libstellar.supernovaeIIYields(filename, m1, mass)
+ mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+ filename, m1, mass)
+ mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+ filename, m1, mass)
+
+ yields_snia = libstellar.supernovaeIaYields(filename)
+ mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
+
+ elements = libstellar.getElements(filename)
+
+ # units
+ mass /= solMass_in_cgs / unit_mass
+
+ plt.figure()
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
+ for i in range(len(elements)):
+ p = plt.loglog(mass, yields_snii[:, i], label=elements[i])
+ plt.loglog(m_wd, yields_snia[i], "o", color=p[0].get_color())
+
+ p = plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+ plt.loglog(m_wd, mass_ejected_snia, "o", color=p[0].get_color())
+
+ plt.loglog(mass, mass_ejected_np_snii,
+ label="Non Processed Mass Ejected")
+
+ plt.xlabel("$Mass [M_\odot]$")
+ plt.ylabel("Mass fraction")
+ plt.legend()
+
+ plt.show()
diff --git a/examples/supernovae_yields/run.sh b/examples/supernovae_yields/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..358d7c530922339afb1903ef7521946d5f6ebeb2
--- /dev/null
+++ b/examples/supernovae_yields/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python plot_yields.py
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 7852bc713c33d950db8883ca686b71f12e998b4f..de585175385c46080183b649f8d07a63798abf10 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):
@@ -78,7 +98,7 @@ def parseYamlFile(filename):
>>> unit_mass = float(us["UnitMass_in_cgs"])
"""
with open(filename, "r") as stream:
- stream = yaml.load(stream)
+ stream = yaml.load(stream, Loader=yaml.SafeLoader)
return stream
@@ -139,39 +159,28 @@ class ParameterFile:
},
"GrackleCooling": {
- "CloudyTable": "CloudyData_UVB=HM2012.h5",
- "WithUVbackground": 1,
- "Redshift": 0,
- "WithMetalCooling": 1,
- "ProvideVolumetricHeatingRates": 0,
- "ProvideSpecificHeatingRates": 0,
- "SelfShieldingMethod": 0,
- "ConvergenceLimit": 1e-2,
- "MaxSteps": 20000,
- "Omega": 0.8
+ "cloudy_table": "CloudyData_UVB=HM2012.h5",
+ "with_UV_background": 1,
+ "redshift": 0,
+ "with_metal_cooling": 1,
+ "provide_volumetric_heating_rates": 0,
+ "provide_specific_heating_rates": 0,
+ "self_shielding_method": 0,
+ "convergence_limit": 1e-2,
+ "max_steps": 20000,
+ "omega": 0.8,
+ "thermal_time_myr": 5,
},
- "GearChemistry": {
- "InitialMetallicity": 0.01295,
- },
-
- "GEARInitialMassFunction": {
- "number_function_part": 4,
- "exponents": [0.7, -0.8, -1.7, -1.3],
- "mass_limits": [0.05, 0.08, 0.5, 1.0, 50.]
+ "GEARChemistry": {
+ "initial_metallicity": 0.01295,
},
"GEARFeedback": {
- "SupernovaeEnergy_erg": 1e51,
- "ThermalTime_Myr": 5,
- "ChemistryTable": "chemistry.hdf5",
+ "supernovae_energy_erg": 1e51,
+ "yields_table": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
+ "discrete_yields": 0
},
-
- "GEARLifetime": {
- "quadratic": [-40.110, 5.509, 0.7824],
- "linear": [141.929, -15.889, -3.2557],
- "constant": [-261.365, 17.073, 9.8661],
- }
}
"""
Dictionnary containing all the default parameters to write in the
diff --git a/scripts/pyswift_cooling_rate b/scripts/pyswift_cooling_rate
index 262629c6debc640b90821d828f9190ef5b67d872..7c905944b6bda576c1048fcc11828912810c8f9e 100644
--- a/scripts/pyswift_cooling_rate
+++ b/scripts/pyswift_cooling_rate
@@ -25,7 +25,6 @@ import numpy as np
import matplotlib.pyplot as plt
from astropy import units, constants
import yaml
-from time import strftime, gmtime
from argparse import ArgumentParser
import os
@@ -244,8 +243,6 @@ def plot_solution(rate, data, us, opt):
cbar.set_ticklabels(tc)
cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
- date = strftime("%d %b %Y", gmtime())
- plt.title("Date: {}".format(date))
plt.show()
diff --git a/scripts/pyswift_initial_mass_function b/scripts/pyswift_initial_mass_function
index 803e520495fc9d453963975de2e4100d20374bc5..d62cad5d5dd780f808b43c2f3c6bae9e35c3973e 100644
--- a/scripts/pyswift_initial_mass_function
+++ b/scripts/pyswift_initial_mass_function
@@ -26,6 +26,7 @@ from argparse import ArgumentParser
import os
plt.rc('text', usetex=True)
+solMass_in_cgs = 1.98848e33
def parseArguments():
@@ -48,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()
@@ -65,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(
@@ -74,24 +75,22 @@ 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"])
- 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)
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
+ mass /= solMass_in_cgs / unit_mass
+
# plot IMF
plt.loglog(mass, imf)
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
index 0b14c53f589d9a26c21841b4393536a000dc1eec..da0a88dfd0fae15762f91f070ab0639813e767d2 100644
--- a/scripts/pyswift_lifetime
+++ b/scripts/pyswift_lifetime
@@ -26,9 +26,15 @@ from argparse import ArgumentParser
import os
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 initial mass function")
+ parser = ArgumentParser(description="Plot the lifetime of a star")
parser.add_argument(
"--download_table", dest="table",
@@ -47,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(
@@ -69,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(
@@ -78,28 +84,30 @@ 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"])
- 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)
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = float(units["UnitMass_in_cgs"])
+ unit_vel = float(units["UnitVelocity_in_cgs"])
+ unit_length = float(units["UnitLength_in_cgs"])
+ unit_time = unit_length / unit_vel
+
+ mass *= solMass_in_cgs / unit_mass
+
# get metallicity
Z = opt.metallicity * np.ones(opt.N, dtype=np.float32)
print("Computing lifetime...")
lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
- lifetime /= 1e6
+
+ mass /= solMass_in_cgs / unit_mass
+ lifetime /= 1e6 * year_in_cgs / unit_time
plt.loglog(mass, lifetime)
diff --git a/scripts/pyswift_plot_yields b/scripts/pyswift_plot_yields
new file mode 100644
index 0000000000000000000000000000000000000000..469cd26f132089b2d8be50a389a2d2808e92c7d8
--- /dev/null
+++ b/scripts/pyswift_plot_yields
@@ -0,0 +1,142 @@
+#!/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 time import strftime, gmtime
+from argparse import ArgumentParser
+import os
+plt.rc('text', usetex=True)
+
+
+N = 1000
+m_min = 0.08 # in solar mass
+m_max = 50
+m_wd = 1.4 # white dwarf mass
+solMass_in_cgs = 1.989e33
+
+
+def parseArguments():
+ parser = ArgumentParser(description="Plot the supernovae rate")
+
+ 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(
+ "--m_min", dest="m_min",
+ type=float, default=0.08,
+ help="Minimal mass to plot in solar mass")
+
+ parser.add_argument(
+ "--m_max", dest="m_max",
+ type=float, default=50,
+ help="Maximal mass to plot in solar mass")
+
+ parser.add_argument(
+ "--disable-snia", dest="disable_snia",
+ action="store_true",
+ help="Disable the SNIa")
+
+ parser.add_argument(
+ "--disable-snii", dest="disable_snii",
+ action="store_true",
+ help="Disable the SNII")
+
+ args = parser.parse_args()
+
+ return args
+
+
+if __name__ == "__main__":
+ opt = parseArguments()
+
+ if opt.table:
+ pyswiftsim.downloadYieldsTable()
+
+ 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)
+
+ # Get masses
+ mass = np.logspace(np.log10(opt.m_min), np.log10(opt.m_max), N,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = float(units["UnitMass_in_cgs"])
+ unit_vel = float(units["UnitVelocity_in_cgs"])
+ unit_length = float(units["UnitLength_in_cgs"])
+ unit_time = unit_length / unit_vel
+
+ mass *= solMass_in_cgs / unit_mass
+
+ m1 = np.ones(N, dtype=np.float32) * 1e-20 # avoid log(0)
+
+ print("Computing supernovae yields...")
+
+ yields_snii = libstellar.supernovaeIIYields(filename, m1, mass)
+ mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+ filename, m1, mass)
+ mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+ filename, m1, mass)
+
+ yields_snia = libstellar.supernovaeIaYields(filename)
+ mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
+
+ elements = libstellar.getElements(filename)
+
+ # units
+ mass /= solMass_in_cgs / unit_mass
+
+ plt.figure()
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
+ for i in range(len(elements)):
+ p = plt.loglog(mass, yields_snii[:, i], label=elements[i])
+ plt.loglog(m_wd, yields_snia[i], "o", color=p[0].get_color())
+
+ p = plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+ plt.loglog(m_wd, mass_ejected_snia, "o", color=p[0].get_color())
+
+ plt.loglog(mass, mass_ejected_np_snii,
+ label="Non Processed Mass Ejected")
+
+ plt.xlabel("$Mass [M_\odot]$")
+ plt.ylabel("Mass fraction")
+ plt.legend()
+
+ plt.show()
diff --git a/scripts/pyswift_supernovae_rate b/scripts/pyswift_supernovae_rate
new file mode 100644
index 0000000000000000000000000000000000000000..93e461f8b4d6d112839968412acd470401e800a2
--- /dev/null
+++ b/scripts/pyswift_supernovae_rate
@@ -0,0 +1,140 @@
+#!/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)
+
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
+
+
+def parseArguments():
+ parser = ArgumentParser(description="Plot the supernovae rate")
+
+ 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(
+ "--time_min", dest="time_min",
+ type=float, default=1.,
+ help="Minimal time to plot in Myr")
+
+ parser.add_argument(
+ "--time_max", dest="time_max",
+ type=float, default=3e4,
+ help="Maximal time to plot in Myr")
+
+ parser.add_argument(
+ "--metallicity", dest="metallicity",
+ type=float, default=0.02041,
+ help="Metallicity of the stars")
+
+ parser.add_argument(
+ "--disable-snia", dest="disable_snia",
+ action="store_true",
+ help="Disable the SNIa")
+
+ parser.add_argument(
+ "--disable-snii", dest="disable_snii",
+ action="store_true",
+ help="Disable the SNII")
+
+ args = parser.parse_args()
+
+ return args
+
+
+if __name__ == "__main__":
+ opt = parseArguments()
+
+ if opt.table:
+ pyswiftsim.downloadYieldsTable()
+
+ 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)
+
+ # Get times in yr
+ opt.time_min *= 1e6
+ opt.time_max *= 1e6
+
+ times = np.logspace(
+ np.log10(opt.time_min), np.log10(opt.time_max), opt.N + 1,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = float(units["UnitMass_in_cgs"])
+ unit_vel = float(units["UnitVelocity_in_cgs"])
+ unit_length = float(units["UnitLength_in_cgs"])
+ unit_time = unit_length / unit_vel
+
+ times *= year_in_cgs / unit_time
+
+ t1 = times[:-1]
+ t2 = times[1:]
+
+ # get metallicity
+ Z = opt.metallicity * np.ones(opt.N, dtype=np.float32)
+
+ print("Computing supernovae rate...")
+
+ if not opt.disable_snia:
+ rate_snia = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+ if not opt.disable_snii:
+ rate_snii = libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+ t1 /= 1e6 * year_in_cgs / unit_time
+
+ unit_rate = year_in_cgs / unit_time
+ unit_rate *= solMass_in_cgs / unit_mass
+
+ if not opt.disable_snia:
+ rate_snia *= 1e9 * unit_rate
+ plt.loglog(t1, rate_snia, label="SNIa")
+ if not opt.disable_snii:
+ rate_snii *= 1e9 * unit_rate
+ plt.loglog(t1, rate_snii, label="SNII")
+
+ plt.xlabel("Time [Myr]")
+ plt.ylabel("Supernovae rate [$M_\odot^{-1} Gyr^{-1}$]")
+ plt.legend()
+
+ plt.show()
diff --git a/setup.py b/setup.py
index 62a38d6006d0bd59d22d2e1efdd3254612664d7c..c8c1ee94860e164a356784a8dd1a3fbc6e02f3ca 100644
--- a/setup.py
+++ b/setup.py
@@ -80,9 +80,11 @@ if swift_path:
])
# hdf5
+ if "/usr/bin" in hdf5_root:
+ include.append("/usr/include/hdf5/serial")
+
if "bin" in hdf5_root:
hdf5_root = hdf5_root.split("bin")[0]
- include.append(hdf5_root + "/include")
# grackle
if parseCmdLine("--with-grackle", store=True):
@@ -101,6 +103,8 @@ lib_dir = []
if swift_path:
lib_dir.append(swift_path + "/src/.libs")
lib_dir.append(hdf5_root + "/lib")
+ # fix for default install of hdf5
+ lib_dir.append(hdf5_root + "/lib/x86_64-linux-gnu/hdf5")
# src files
c_src = []
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 4fa9731fb44d1a66db652c95c553ab590ddc8dd9..dbe718c9f17cf5c7f5ef01f313b5d694bd60ca0c 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -19,11 +19,6 @@
#include "cooling_wrapper.h"
#include "pyswiftsim_tools.h"
-/* TODO Remove this */
-struct cooling_function_data cooling;
-int initialized;
-
-
/**
* @brief Set the cooling element fractions
*
@@ -79,7 +74,6 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
*/
PyObject* pycooling_rate(PyObject* self, PyObject* args) {
import_array();
-
char *filename;
@@ -125,7 +119,6 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
struct chemistry_global_data chemistry;
chemistry_init_backend(¶ms, &us, &pconst, &chemistry);
-
/* Init entropy floor */
struct entropy_floor_properties floor_props;
entropy_floor_init(&floor_props, &pconst, &us,
@@ -135,7 +128,7 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
struct xpart xp;
hydro_first_init_part(&p, &xp);
-
+
/* return object */
PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
@@ -150,20 +143,24 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
if (fractions != NULL)
- pycooling_set_fractions(&xp, fractions, i);
+ pycooling_set_fractions(&xp, fractions, i);
else {
- /* Need to set the mass in order to estimate the element fractions */
- hydro_set_mass(&p, p.rho);
- p.h = 1;
+ /* Need to set the mass in order to estimate the element fractions */
+ hydro_set_mass(&p, p.rho);
+ p.h = 1;
- cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+ cooling_first_init_part(&pconst, &us, &hydro_props, &cosmo, &cooling, &p, &xp);
}
hydro_set_physical_internal_energy_dt(&p, &cosmo, 0);
+ for(int j = 0; j < GEAR_CHEMISTRY_ELEMENT_COUNT; j++) {
+ p.chemistry_data.smoothed_metal_mass_fraction[j] = p.chemistry_data.metal_mass_fraction[j];
+ }
+
/* compute cooling rate */
cooling_cool_part(&pconst, &us, &cosmo, &hydro_props,
- &floor_props, &cooling, &p, &xp, dt, dt);
+ &floor_props, &cooling, &p, &xp, /* Time */0, dt, dt);
float *tmp = PyArray_GETPTR1(rate, i);
*tmp = hydro_get_physical_internal_energy_dt(&p, &cosmo);
@@ -171,6 +168,7 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
/* Cleanup */
cosmology_clean(&cosmo);
+ cooling_clean(&cooling);
return (PyObject *) rate;
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index e0b97a9c9f70d9033dfaac54ae75bb061f8dfb35..8491ff58ae39a044d0fa2e0f979f56ccb3a96481 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -21,26 +21,12 @@
#include "pyswiftsim_includes.h"
-/* A few global variables */
-extern struct cooling_function_data cooling;
-extern int initialized;
-
-
PyObject* pycooling_rate(PyObject* self, PyObject* args);
-#define PYSWIFTSIM_INIT_COOLING() \
- \
- /* init cooling */ \
- if (initialized == 0) { \
- /* struct cooling_function_data cooling; */ \
- cooling_init_backend(¶ms, &us, &pconst, &cooling); \
- initialized = 1; \
- } \
- else if (initialized == 1) { \
- message("WARNING: not reinitializing the cooling, need to be fixed"); \
- initialized = 2; \
- } \
- \
+#define PYSWIFTSIM_INIT_COOLING() \
+ struct cooling_function_data cooling; \
+ cooling_init_backend(¶ms, &us, &pconst, &hydro_props, \
+ &cooling); \
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/libstellar.c b/src/libstellar.c
index fcd1cb304d5207fd409eb9461337e0bac03ac60e..412bd143057329c9e71ea02ec5db430e5e32234c 100644
--- a/src/libstellar.c
+++ b/src/libstellar.c
@@ -96,7 +96,181 @@ static PyMethodDef libstellar_methods[] = {
">>> mass = libstellar.massFromLifetime(filename, lifetime, Z)\n"
},
+ {"supernovaeIaRate", pysupernovae_ia_rate, METH_VARARGS,
+ "Compute the supernovae Ia rate from the time and the stellar metallicity.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "time1: array\n"
+ "\t Beginning of the time interval for the SNIa in internal units.\n"
+ "time2: array\n"
+ "\t End of the time interval for the SNIa in internal units.\n"
+ "metallicity: array\n"
+ "\t Metallicity of the stars.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "rate: array\n"
+ "\t The rate of SNIa.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> time1 = np.array([1e7, 1e9])\n"
+ ">>> time2 = np.array([2e7, 2e9])\n"
+ ">>> Z = np.array([0.02, 0.01])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> mass = libstellar.supernovaeIaRate(filename, time1, time2, Z)\n"
+ },
+
+ {"supernovaeIIRate", pysupernovae_ii_rate, METH_VARARGS,
+ "Compute the supernovae II rate from the time and the stellar metallicity.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "time1: array\n"
+ "\t Beginning of the time interval for the SNII in internal units.\n"
+ "time2: array\n"
+ "\t End of the time interval for the SNII in internal units.\n"
+ "metallicity: array\n"
+ "\t Metallicity of the stars.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "rate: array\n"
+ "\t The rate of SNII.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> time1 = np.array([1e7, 1e9])\n"
+ ">>> time2 = np.array([2e7, 2e9])\n"
+ ">>> Z = np.array([0.02, 0.01])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> mass = libstellar.supernovaeIIRate(filename, time1, time2, Z)\n"
+ },
+
+ {"supernovaeIIYields", pysupernovae_ii_yields, METH_VARARGS,
+ "Compute the supernovae II yields from the masses.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "mass1: array\n"
+ "\t Beginning of the mass interval for the SNII in internal units.\n"
+ "mass2: array\n"
+ "\t End of the time interval for the SNII in internal units.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "yields: array\n"
+ "\t The yields (mass fraction) of SNII.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass1 = np.array([1.1, 1.3])\n"
+ ">>> mass2 = np.array([1.2, 1.4])\n"
+ ">>> mass = libstellar.supernovaeIIYields('params.yml', mass1, mass2)\n"
+ },
+ {"supernovaeIaYields", pysupernovae_ia_yields, METH_VARARGS,
+ "Compute the supernovae Ia yields from the masses.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "yields: array\n"
+ "\t The yields (mass fraction) of SNIa.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass = libstellar.supernovaeIaYields('params.yml')\n"
+ },
+
+ {"getElements", pystellar_get_elements, METH_VARARGS,
+ "Read the element names..\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "element_names: array\n"
+ "\t The name of all the elements.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass = libstellar.getElements('params.yml')\n"
+ "[Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu, Metals]"
+ },
+
+ {"supernovaeIIMassEjected", pysupernovae_ii_mass_ejected, METH_VARARGS,
+ "Compute the mass ejected by supernovae II from the masses.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "mass1: array\n"
+ "\t Beginning of the mass interval for the SNII in internal units.\n"
+ "mass2: array\n"
+ "\t End of the time interval for the SNII in internal units.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "mass_ejected: array\n"
+ "\t The mass ejected by the supernovae SNII.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass1 = np.array([1.1, 1.3])\n"
+ ">>> mass2 = np.array([1.2, 1.4])\n"
+ ">>> mass = libstellar.supernovaeIIMassEjected('params.yml', mass1, mass2)\n"
+ },
+
+ {"supernovaeIIMassEjectedNonProcessed", pysupernovae_ii_mass_ejected_non_processed, METH_VARARGS,
+ "Compute the (non processed) mass ejected by supernovae II from the masses.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "mass1: array\n"
+ "\t Beginning of the mass interval for the SNII in internal units.\n"
+ "mass2: array\n"
+ "\t End of the time interval for the SNII in internal units.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "mass_ejected: array\n"
+ "\t The non processed mass ejected by the supernovae SNII.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass1 = np.array([1.1, 1.3])\n"
+ ">>> mass2 = np.array([1.2, 1.4])\n"
+ ">>> mass = libstellar.supernovaeIIMassEjectedNonProcessed('params.yml', mass1, mass2)\n"
+ },
+
+ {"supernovaeIaMassEjected", pysupernovae_ia_mass_ejected, METH_VARARGS,
+ "Compute the mass ejected by supernovae Ia.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "mass_ejected: float\n"
+ "\t The mass ejected by the supernovae SNIa.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass = libstellar.supernovaeIaMassEjected('params.yml')\n"
+ },
+
{NULL, NULL, 0, NULL} /* Sentinel */
};
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
index a990b753d2c05ddc620f2115bce366a75fc4ce09..3a23c8d0fd80f6fbff27cc02c2b06e8fc267cf28 100644
--- a/src/stellar_evolution_wrapper.c
+++ b/src/stellar_evolution_wrapper.c
@@ -63,8 +63,18 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
float m = *(float*) PyArray_GETPTR1(mass, i);
float *imf_i = (float*) PyArray_GETPTR1(imf, i);
+ /* change units internal units -> solar mass */
+ m /= pconst.const_solar_mass;
- *imf_i = stellar_evolution_get_imf(&fp.stellar_model.imf, m);
+ if (m > fp.stellar_model.imf.mass_max) {
+ *imf_i = 0.;
+ continue;
+ }
+
+ *imf_i = initial_mass_function_get_imf(&fp.stellar_model.imf, m);
+
+ /* change from internal units to solar mass */
+ *imf_i /= pconst.const_solar_mass;
}
return (PyObject *) imf;
@@ -75,11 +85,13 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
* @brief Compute the lifetime of a star.
*
* args is expecting:
- * TODO
+ * - the name of the parameter file
+ * - a numpy array containing the masses.
+ * - a numpy array containing the metallicity.
*
* @param self calling object
* @param args arguments
- * @return cooling rate
+ * @return the lifetime of the stars
*/
PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
import_array();
@@ -122,9 +134,12 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
float z = *(float*) PyArray_GETPTR1(metallicity, i);
float *tau = (float*) PyArray_GETPTR1(lifetime, i);
+
+ /* change units internal units -> solar mass */
+ m /= pconst.const_solar_mass;
- *tau = stellar_evolution_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
- *tau = pow(10, *tau);
+ *tau = lifetime_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
+ *tau = pow(10, *tau) * pconst.const_year * 1e6;
}
return (PyObject *) lifetime;
@@ -135,11 +150,13 @@ PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
* @brief Compute the mass of a star with a given lifetime.
*
* args is expecting:
- * TODO
+ * - the name of the parameter file
+ * - a numpy array containing the time.
+ * - a numpy array containing the metallicity.
*
* @param self calling object
* @param args arguments
- * @return cooling rate
+ * @return The mass of the stars.
*/
PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
import_array();
@@ -183,10 +200,513 @@ PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
float *m = (float*) PyArray_GETPTR1(mass, i);
- *m = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
- *m = pow(10, *m);
+ /* change units internal units to Myr */
+ tau /= pconst.const_year * 1e6;
+
+ *m = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
+ *m = pow(10, *m) * pconst.const_solar_mass;
}
return (PyObject *) mass;
}
+
+
+
+/**
+ * @brief Compute the supernovae Ia rate.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the lower limit for the time.
+ * - a numpy array containing the upper limit for the time.
+ * - a numpy array containing the metallicity.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return Rate of supernovae Ia
+ */
+PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *time1;
+ PyArrayObject *time2;
+ PyArrayObject *metallicity;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOOO", &filename, &time1, &time2, &metallicity))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(time1, 1, NPY_FLOAT, "Time1") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(time2, 1, NPY_FLOAT, "Time2") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time1, 0))
+ error("Time1 and metallicity should have the same dimension");
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time2, 0))
+ error("Time2 and metallicity should have the same dimension");
+
+
+ size_t N = PyArray_DIM(time1, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(time1), PyArray_DIMS(time1), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float t1 = *(float*) PyArray_GETPTR1(time1, i);
+ float t2 = *(float*) PyArray_GETPTR1(time2, i);
+ float z = *(float*) PyArray_GETPTR1(metallicity, i);
+ float *r = (float*) PyArray_GETPTR1(rate, i);
+
+ /* change units internal units -> Myr */
+ t1 /= pconst.const_year * 1e6;
+ t2 /= pconst.const_year * 1e6;
+
+ float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
+ m1 = pow(10, m1);
+ float m2 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
+ m2 = pow(10, m2);
+
+ /* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
+ *r = supernovae_ia_get_number_per_unit_mass(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
+
+ /* change units 1 / (solMass Myr) -> internal units */
+ *r /= pconst.const_solar_mass * pconst.const_year * 1e6;
+ }
+
+ return (PyObject *) rate;
+
+}
+
+/**
+ * @brief Compute the supernovae II rate.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the lower limit for the time.
+ * - a numpy array containing the upper limit for the time.
+ * - a numpy array containing the metallicity.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return Rate of supernovae II
+ */
+PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *time1;
+ PyArrayObject *time2;
+ PyArrayObject *metallicity;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOOO", &filename, &time1, &time2, &metallicity))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(time1, 1, NPY_FLOAT, "Time1") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(time2, 1, NPY_FLOAT, "Time2") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time1, 0))
+ error("Time1 and metallicity should have the same dimension");
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time2, 0))
+ error("Time2 and metallicity should have the same dimension");
+
+
+ size_t N = PyArray_DIM(time1, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(time1), PyArray_DIMS(time1), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float t1 = *(float*) PyArray_GETPTR1(time1, i);
+ float t2 = *(float*) PyArray_GETPTR1(time2, i);
+ float z = *(float*) PyArray_GETPTR1(metallicity, i);
+ float *r = (float*) PyArray_GETPTR1(rate, i);
+
+ /* change units internal units -> Myr */
+ t1 /= pconst.const_year * 1e6;
+ t2 /= pconst.const_year * 1e6;
+
+ float m1 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t1), z);
+ m1 = pow(10, m1);
+ float m2 = lifetime_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(t2), z);
+ m2 = pow(10, m2);
+
+ /* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
+ *r = supernovae_ii_get_number_per_unit_mass(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
+
+ /* change units 1 / (solMass Myr) -> internal units */
+ *r /= pconst.const_solar_mass * pconst.const_year * 1e6;
+
+ }
+
+ return (PyObject *) rate;
+
+}
+
+/**
+ * @brief Compute the supernovae II yields.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the lower limit for the mass.
+ * - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction of ejected metals
+ */
+PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *mass1;
+ PyArrayObject *mass2;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO", &filename, &mass1, &mass2))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+ error("The inputs should be 1-D arrays");
+
+ if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+ error("Mass1 and Mass2 should have the same dimension");
+
+
+ size_t N = PyArray_DIM(mass1, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ const int ndims = 2;
+ npy_intp dims[2] = {N, GEAR_CHEMISTRY_ELEMENT_COUNT};
+ PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float yields[GEAR_CHEMISTRY_ELEMENT_COUNT] = {0.};
+ const float m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
+ const float log_m1 = log10(m1);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
+ const float log_m2 = log10(m2);
+
+ /* Compute the yields */
+ supernovae_ii_get_yields(&fp.stellar_model.snii, log_m1, log_m2, yields);
+
+ for(int j = 0; j < GEAR_CHEMISTRY_ELEMENT_COUNT; j++) {
+ float *y = (float*) PyArray_GETPTR2(yields_out, i, j);
+ *y = yields[j];
+ }
+ }
+
+ return (PyObject *) yields_out;
+
+}
+
+/**
+ * @brief Compute the mass ejected by a supernovae II.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the lower limit for the mass.
+ * - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction of ejected mass.
+ */
+PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *mass1;
+ PyArrayObject *mass2;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO", &filename, &mass1, &mass2))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+ error("The inputs should be 1-D arrays");
+
+ if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+ error("Mass1 and Mass2 should have the same dimension");
+
+
+ size_t N = PyArray_DIM(mass1, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(
+ PyArray_NDIM(mass1), PyArray_DIMS(mass1), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ const float m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
+ const float log_m1 = log10(m1);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
+ const float log_m2 = log10(m2);
+ float *y = (float*) PyArray_GETPTR1(mass_ejected, i);
+
+ *y = supernovae_ii_get_ejected_mass_fraction_processed(&fp.stellar_model.snii, log_m1, log_m2);
+
+ }
+
+ return (PyObject *) mass_ejected;
+
+}
+
+/**
+ * @brief Compute the non processed mass ejected by a supernovae II.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the lower limit for the mass.
+ * - a numpy array containing the upper limit for the mass.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass ratio of (non processed) mass ejected
+ */
+PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *mass1;
+ PyArrayObject *mass2;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO", &filename, &mass1, &mass2))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(mass1, 1, NPY_FLOAT, "Mass1") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(mass2, 1, NPY_FLOAT, "Mass2") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_NDIM(mass1) != 1 || PyArray_NDIM(mass2) != 1)
+ error("The inputs should be 1-D arrays");
+
+ if (PyArray_DIM(mass1, 0) != PyArray_DIM(mass2, 0))
+ error("Mass1 and Mass2 should have the same dimension");
+
+
+ size_t N = PyArray_DIM(mass1, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(
+ PyArray_NDIM(mass1), PyArray_DIMS(mass1), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ const float m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
+ const float log_m1 = log10(m1);
+ const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
+ const float log_m2 = log10(m2);
+ float *y = (float*) PyArray_GETPTR1(mass_ejected, i);
+
+ *y = supernovae_ii_get_ejected_mass_fraction(&fp.stellar_model.snii, log_m1, log_m2);
+
+ }
+
+ return (PyObject *) mass_ejected;
+
+}
+
+/**
+ * @brief Compute the supernovae Ia yields.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return Mass fraction of the supernovae Ia yields.
+ */
+PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "s", &filename))
+ return NULL;
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ const int ndims = 1;
+ npy_intp dims[1] = {GEAR_CHEMISTRY_ELEMENT_COUNT};
+ PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+ const float *yields = supernovae_ia_get_yields(&fp.stellar_model.snia);
+
+ for(size_t i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+ float *y = (float*) PyArray_GETPTR1(yields_out, i);
+
+ *y = yields[i] / fp.stellar_model.snia.mass_white_dwarf;
+ }
+
+ return (PyObject *) yields_out;
+
+}
+
+/**
+ * @brief Compute the mass ejected by a supernovae Ia.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return mass fraction ejected.
+ */
+PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "s", &filename))
+ return NULL;
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ const int ndims = 1;
+ npy_intp dims[1] = {1};
+ PyArrayObject *mass_ejected = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
+
+ float *y = (float*) PyArray_GETPTR1(mass_ejected, 0);
+
+ *y = supernovae_ia_get_ejected_mass_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
+
+ return (PyObject *) mass_ejected;
+
+}
+
+
+/**
+ * @brief Read the elements names.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pystellar_get_elements(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(args, "s", &filename))
+ return NULL;
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* Save the variables */
+ PyObject *t = PyList_New(GEAR_CHEMISTRY_ELEMENT_COUNT);
+
+ for(int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+ const char *name = stellar_evolution_get_element_name(&fp.stellar_model, i);
+ PyObject *string = PyUnicode_FromString(name);
+ int test = PyList_SetItem(t, i, string);
+ if (test != 0) {
+ error("Failed to add an element");
+ }
+ }
+
+ return t;
+
+}
diff --git a/src/stellar_evolution_wrapper.h b/src/stellar_evolution_wrapper.h
index a76b04d39b1733a8667cc9435859cfeb9fbfbeea..2c3725eeab2837c8d152c016d787e3cd49f26512 100644
--- a/src/stellar_evolution_wrapper.h
+++ b/src/stellar_evolution_wrapper.h
@@ -24,6 +24,15 @@
PyObject* pyinitial_mass_function(PyObject* self, PyObject* args);
PyObject* pylifetime_from_mass(PyObject* self, PyObject* args);
PyObject* pymass_from_lifetime(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_mass_ejected(PyObject* self, PyObject* args);
+PyObject* pysupernovae_ii_mass_ejected_non_processed(PyObject* self, PyObject* args);
+PyObject* pystellar_get_elements(PyObject* self, PyObject* args);
+
#endif // __PYSWIFTSIM_STELLAR_H__
diff --git a/tests/test_initial_mass_function.py b/tests/test_initial_mass_function.py
index ca80c505bb39844035953934cdb84bf221591068..d77b0e24b2b21f831cc5d05b8ec58ee8fe831d11 100644
--- a/tests/test_initial_mass_function.py
+++ b/tests/test_initial_mass_function.py
@@ -26,22 +26,29 @@ 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"]
- mass_min = float(mass_limits[0])
- mass_max = float(mass_limits[1])
mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ # change the units
+ mass *= solMass_in_cgs / unit_mass
+
print("Computing initial mass function...")
imf = libstellar.initialMassFunction(filename, mass)
+ imf *= solMass_in_cgs / unit_mass
+
print("IMF obtained: {}".format(imf))
diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py
index 9a1816e22438b5a857a9945e571961af280cccbc..cc1f75ef9f0ae453ba8065f3d6b20d548fc1b9c7 100644
--- a/tests/test_lifetime.py
+++ b/tests/test_lifetime.py
@@ -27,20 +27,26 @@ plt.rc('text', usetex=True)
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"]
- mass_min = float(mass_limits[0])
- mass_max = float(mass_limits[1])
mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+ # change the units
+ mass *= solMass_in_cgs / unit_mass
+
# get metallicity
Z = solar_metallicity * np.ones(N, dtype=np.float32)
@@ -49,6 +55,10 @@ if __name__ == "__main__":
lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
+ mass /= solMass_in_cgs / unit_mass
+ mass_from_lifetime /= solMass_in_cgs / unit_mass
+ lifetime /= year_in_cgs / unit_time
+
print("Mass used: {}".format(mass))
print("Lifetime obtained: {}".format(lifetime))
diff --git a/tests/test_rates.py b/tests/test_rates.py
new file mode 100644
index 0000000000000000000000000000000000000000..6d879d62712ced21ba6a98a9d83e5d616ab946f1
--- /dev/null
+++ b/tests/test_rates.py
@@ -0,0 +1,73 @@
+#!/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
+t_min = 1e6 # in yr
+t_max = 1e10
+solMass_in_cgs = 1.989e33
+year_in_cgs = 3.15569252e7
+
+pyswiftsim.downloadYieldsTable()
+
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+ parser = pyswiftsim.parseYamlFile(filename)
+
+ # Get times
+ times = np.logspace(np.log10(t_min), np.log10(t_max), N + 1,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+
+ times *= year_in_cgs / unit_time
+
+ t1 = times[:-1]
+ t2 = times[1:]
+
+ # get metallicity
+ Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+ print("Computing supernovae rate...")
+
+ rate = libstellar.supernovaeIaRate(filename, t1, t2, Z)
+ rate += libstellar.supernovaeIIRate(filename, t1, t2, Z)
+
+ rate *= year_in_cgs / unit_time
+ rate *= solMass_in_cgs / unit_mass
+
+ t1 /= year_in_cgs / unit_time
+
+ print("Time used: {}".format(t1))
+
+ print("Rate obtained: {}".format(rate))
diff --git a/tests/test_yields.py b/tests/test_yields.py
new file mode 100644
index 0000000000000000000000000000000000000000..0753da7b998257e7c81145a7b69d5e8f741177e8
--- /dev/null
+++ b/tests/test_yields.py
@@ -0,0 +1,80 @@
+#!/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 = 10
+m_min = 1 # in solar mass
+m_max = 50
+solMass_in_cgs = 1.989e33
+
+pyswiftsim.downloadYieldsTable()
+
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+ parser = pyswiftsim.parseYamlFile(filename)
+
+ # Get masses
+ mass = np.logspace(np.log10(m_min), np.log10(m_max), N + 1,
+ dtype=np.float32)
+
+ # transform into internal units
+ units = parser["InternalUnitSystem"]
+ unit_mass = units["UnitMass_in_cgs"]
+ unit_vel = units["UnitVelocity_in_cgs"]
+ unit_length = units["UnitLength_in_cgs"]
+ unit_time = unit_length / unit_vel
+
+ mass *= solMass_in_cgs / unit_mass
+
+ m1 = mass[:-1]
+ m2 = mass[1:]
+
+ print("Computing supernovae yields...")
+
+ yields_snii = libstellar.supernovaeIIYields(filename, m1, m2)
+ mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+ filename, m1, m2)
+ mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+ filename, m1, m2)
+
+ m1 /= solMass_in_cgs / unit_mass
+
+ yields_snia = libstellar.supernovaeIaYields(filename)
+ mass_ejected_snia = libstellar.supernovaeIaMassEjected(filename)
+
+ print("Mass used: {}".format(m1))
+
+ print("SNII yields obtained: {}".format(yields_snii))
+
+ print("SNIa yields obtained: {}".format(yields_snia))
+
+ print("SNII mass ejected: {}".format(mass_ejected_snii))
+
+ print("SNIa mass ejected: {}".format(mass_ejected_snia))
+
+ print("SNII non processed mass ejected: {}".format(
+ mass_ejected_np_snii))