Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77990208
pychem_sample_IMF_optimal
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
Sat, Aug 17, 17:11
Size
6 KB
Mime Type
text/x-python
Expires
Mon, Aug 19, 17:11 (2 d)
Engine
blob
Format
Raw Data
Handle
19962462
Attached To
rGEAR Gear
pychem_sample_IMF_optimal
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
##################
# optimal imf sampling
##################
Mecl = opt.mstar # cluster mass, or mass of the IM
# init
m_max_star = chemistry.optimal_sampling_init_norm(Mecl)
m = m_max_star
M_diced = 0
ms = []
i = 1
while (M_diced <= Mecl):
M_diced += m
#print m
ms.append(m)
#print i,m
if chemistry.optimal_sampling_stop_loop(m):
break;
m = chemistry.optimal_sampling_get_next_mass(m)
i=i+1
ms_opt = array(ms)
print "m_max = ",m_max_star
##################
# 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 = 100.
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
# histogram opt
n = searchsorted(sort(ms_opt),bins)
n = concatenate([n,[len(ms_opt)]])
hx_opt = n[1:]-n[:-1]
hx_opt = hx_opt*bins
print "total mass (histogram random)",sum(hx)
print "total mass (histogram opt) ",sum(hx_opt)
print "max mass (random) ",ms.max()
print "max mass (opt) ",ms_opt.max()
pt.plot(bins,hx,label='random')
pt.plot(bins,hx_opt,label='optimal')
#pt.plot(bins_th,hx_th,label='theoretical')
pt.legend()
################
# 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()
######################
#
pt.figure()
ms.sort()
x = ms
y = add.accumulate(ms)
pt.plot(x,y,label='random')
ms_opt.sort()
x = ms_opt
y = add.accumulate(ms_opt)
pt.plot(x,y,label='optimal')
pt.legend(loc='lower right')
#pt.semilogy()
pt.xlabel(r'$\rm{Star\,Mass}\,[\rm{M}_\odot]$')
pt.ylabel(r'$\rm{Total\,Mass}$')
pt.show()
#
######################
Event Timeline
Log In to Comment