Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73081180
pychem_sample_IMF
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jul 18, 09:43
Size
5 KB
Mime Type
text/x-python
Expires
Sat, Jul 20, 09:43 (2 d)
Engine
blob
Format
Raw Data
Handle
19138741
Attached To
rGEAR Gear
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)*UnitMass_in_g/SOLAR_MASS
#print min(ms)*UnitMass_in_g/SOLAR_MASS
#print max(ms)*UnitMass_in_g/SOLAR_MASS
##############################################
# 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
Log In to Comment