Page MenuHomec4science
Files F63721568

gsfr
No OneTemporary
Actions

File Metadata

Created
Wed, May 22, 02:06

gsfr
View Options

#!/home/revaz/local/bin/python
'''
Give info on the theoretical star formation rate
Yves Revaz
lun mai 22 15:31:03 CEST 2006
'''
from numarray import *
from Nbody import *
import string
import sys
import os
from libjeans import *
from Nbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import io
from numarray import random_array as RandomArray
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = add_postscript_options(parser)
parser = add_color_options(parser)
parser = add_limits_options(parser)
parser = add_units_options(parser)
parser = add_log_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--cgs",
action="store_true",
dest="cgs",
default = 0,
help="invert into cgs units")
parser.add_option("-r",
action="store",
dest="reduc",
type="int",
default = 1,
help="reduction number",
metavar=" INT")
(options, args) = parser.parse_args()
if options.colors!=None:
exec("options.colors = array([%s])"%(options.colors))
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#############################
# graph
#############################
# get options
files,options = parse_options()
ps = options.ps
col = options.colors
xmin = options.xmin
xmax = options.xmax
ymin = options.ymin
ymax = options.ymax
log = options.log
cgs = options.cgs
ftype = options.ftype
reduc = options.reduc
localsystem = Set_SystemUnits_From_Options(options)
#######################################
# define output system of unit
#######################################
outputunits = UnitSystem('mks',[Unit_kpc, Unit_Ms, Unit_yr , Unit_K, Unit_mol, Unit_C])
#######################################
# LOOP
#######################################
# read files
for file in files:
nb = Nbody(file,ftype=ftype)
nb = nb.select('gas')
mass = nb.mass.astype(Float)
rho = nb.rho.astype(Float)
StarFormationParameter = 0.005
StarFormationTime = 1000.
NStarsFromGas = 1.
dt = 0.038147
# compute sfr1 : bad one
dmdt1 = sum(mass/StarFormationTime)
# compute sfr2
#p = NStarsFromGas*(1. - exp(-dt/StarFormationTime) )
#rand = RandomArray.random(nb.nbody)
#dmdt2 = where(rand<p,mass,0.)
#dmdt2 = sum(dmdt2)/dt
#np = where(rand<p,1,0)
#np = sum(np)
# compute sfr2 : good one
dmdt2 = sum(StarFormationParameter*mass*rho**(1/2.))
# sfr factor
FACT = PhysCte(1,localsystem.UnitDic['kg']/localsystem.UnitDic['s'])
FACT = FACT.into(outputunits)
dmdt1 = dmdt1*FACT
dmdt2 = dmdt2*FACT
print "%8.3f %10.6f %10.6f"%(nb.tnow,dmdt1,dmdt2)

Event Timeline

c4science · Help