Page MenuHomec4science
Files F117332279

pychem_sample_IMF
No OneTemporary
Actions

File Metadata

Created
Fri, Jun 13, 01:59

pychem_sample_IMF
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
print "number of stars",N
print "total mass ",sum(ms)*chemistry.CodeUnits_to_SolarMass_Factor()
#print min(ms)*UnitMass_in_g/SOLAR_MASS
#print max(ms)*UnitMass_in_g/SOLAR_MASS
# nombre de SNs
ms_ns = compress(ms>chemistry.get_SNII_Mmin(),ms)
print "number of SNII ",len(ms_ns)
##############################################
# plot the imf
##############################################
# to ms
ms = ms*UnitMass_in_g/SOLAR_MASS
mmin = mmin *UnitMass_in_g/SOLAR_MASS
mmax = mmax *UnitMass_in_g/SOLAR_MASS
#ms = compress(ms>8,ms)
#print len(ms)
n = 10000.
db = (mmax-mmin)/n
bins = arange(mmin,mmax+db,db)
# histogram
n = searchsorted(sort(ms),bins)
n = concatenate([n,[len(ms)]])
hx = n[1:]-n[:-1]
hx = hx*bins
print "total mass (histogram)",sum(hx)
pt.plot(bins,hx)
pt.plot(bins_th,hx_th)
################
# get values
################
x = ms
y = chemistry.get_imf(x)
pt.loglog()
pt.xlabel(r'$\rm{Mass}\,[\rm{M}_\odot]$')
pt.ylabel(r'$\rm{number\,of\,stars}\,,\,\rm{mass\,fraction}$')
pt.show()

Event Timeline

c4science · Help