plot_initial_mass_function.py

#!/usr/bin/env python3
# ########################################################################
# This file is part of PYSWIFT.
# Copyright (c) 2018 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 = 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)

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

        print("Computing initial mass function...")

        imf = libstellar.initialMassFunction(filename, mass)
        imf *= solMass_in_internal

        # plot IMF
        plt.figure()
        mass /= solMass_in_internal
        plt.loglog(mass, imf)

        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()