Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_sample_IMF_and_properties
b3a44a0e17d6master

Last Change
Raw File

pychem_sample_IMF_and_properties
View Options

#!/usr/bin/env python
import Ptools as pt
from optparse import OptionParser
from pNbody import *
from pNbody import units
import string
from scipy import optimize
from PyChem import chemistry
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_postscript_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'Fe',
help="x value to plot",
metavar=" STRING")
parser.add_option("--y",
action="store",
dest="y",
type="string",
default = 'Mg',
help="y value to plot",
metavar=" STRING")
parser.add_option("--dt",
action="store",
dest="dt",
type="float",
default = 0.1,
help="dt",
metavar=" FLOAT")
parser.add_option("--mstar",
action="store",
dest="mstar",
type="float",
default = 1e5,
help="initial mass of the SSP in solar mass",
metavar=" FLOAT")
parser.add_option("--tstar",
action="store",
dest="tstar",
type="float",
default = 0,
help="formation time of the SSP",
metavar=" FLOAT")
parser.add_option("-o",
action="store",
dest="obs",
type="string",
default = 'Y',
help="observable to plot",
metavar=" STRING")
parser.add_option("--timeunit",
action="store",
dest="timeunit",
type="string",
default = None,
help="unit of time",
metavar=" STRING")
parser.add_option("--NumberOfTables",
action="store",
dest="NumberOfTables",
type="int",
default = 1,
help="NumberOfTables",
metavar=" INT")
parser.add_option("--DefaultTable",
action="store",
dest="DefaultTable",
type="int",
default = 0,
help="DefaultTable",
metavar=" INT")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# M A I N
########################################################################
SOLAR_MASS = 1.989e33
UnitLength_in_cm = 3.085e+21
UnitMass_in_g = 1.989e+43
UnitVelocity_in_cm_per_s = 20725573.785998672
UnitTime_in_s = 148849920000000.0
files,opt = parse_options()
file = files[0]
# some parameters
M0 = opt.mstar *SOLAR_MASS/UnitMass_in_g # gas mass of the SSP (in code unit)
# init
chemistry.init_chimie(file,opt.NumberOfTables,opt.DefaultTable)
mmax = chemistry.get_Mmax()
mmin = chemistry.get_Mmin()
##################
# theoretical imf
##################
'''
here, the imf is normed in order that
the integral of it, i.e, the mass fraction
of stars is equal to 1.
'''
bins_th = arange(mmin,mmax,1e-12).astype(float32)
hx_th = chemistry.get_imf(bins_th).astype(float32)
bins_th = bins_th.astype(float)*UnitMass_in_g/SOLAR_MASS
##################
# sampled imf
##################
'''
here, the imf is normed, in order that the total mass of
the stars = M0
==> the normalisation is different from the theoretical imf.
'''
# compute the number of stars per mass between m1 and m2 (dep on M0)
# N this is thus the number of stars in a particle of mass M0
N = chemistry.get_imf_N(array([mmin]),array([mmax]))[0]*M0
# compute the masses
ms = chemistry.imf_sampling(int(N),1)
# in Msol
Mtot = sum(ms)*chemistry.CodeUnits_to_SolarMass_Factor()
print "number of stars",N
print "total mass %g [Msol] "%Mtot
#print min(ms)*UnitMass_in_g/SOLAR_MASS
#print max(ms)*UnitMass_in_g/SOLAR_MASS
# nombre de SNII
ms_ns = compress(ms>chemistry.get_SNII_Mmin(),ms)
NSNII = len(ms_ns)
MSNII = sum(ms_ns)*chemistry.CodeUnits_to_SolarMass_Factor()
MSNII_MTOT = MSNII/Mtot
print "number of SNII ",NSNII
print "enregy ejected in SNII %g ergs"%(1e51*NSNII)
print "gas mass fraction ejected in SNII %g [Msol]"%MSNII_MTOT
c4science · Help