Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F102365716
old.pRotCurve
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
Wed, Feb 19, 23:17
Size
6 KB
Mime Type
text/x-python
Expires
Fri, Feb 21, 23:17 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
24339525
Attached To
rGTOOLS Gtools
old.pRotCurve
View Options
#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
import string
from scipy import optimize
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parse = pt.add_units_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'r',
help="x value to plot",
metavar=" STRING")
parser.add_option("--nx",
action="store",
type="int",
dest="nx",
default = 14,
help="number of bins")
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
def get_value_and_label(nb,val):
if val == 'R':
label = r'$\rm{Radius}\,\left[ kpc \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Rxyz(units=out_units.UnitLength)
return val,label
elif val == 'T':
label = r'$\rm{Temperature}\,\left[ \rm{K} \right]$'
val = nb.T()
return val,label
elif val == 'logT':
label = r'$\rm{log\,Temperature}\,\left[ \rm{K} \right]$'
val = nb.T()
val = log10(val)
return val,label
elif val == 'rho':
label = r'$\rm{Density}\,\left[ \rm{atom/cm^3} \right]$'
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
Unit_atom.set_symbol('atom')
out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
val = nb.Rho(units=out_units.UnitDensity)
return val,label
elif val == 'logrho':
label = r'$\rm{log\,Density}\,\left[ \rm{atom/cm^3} \right]$'
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
Unit_atom.set_symbol('atom')
out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
val = nb.Rho(units=out_units.UnitDensity)
val = log10(val)
return val,label
elif val == 'Tcool':
label = r'$\rm{Cooling\,Time}\,\left[ \rm{Myr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Tcool(units=out_units.UnitTime)
return val,label
# old implementation
#
# val = nb.Tcool()
#
## convert in Myr
#outputunits = units.UnitSystem('mks',[units.Unit_kpc, units.Unit_Ms, units.Unit_Myr , units.Unit_K, units.Unit_mol, units.Unit_C])
#codeunits = nb.localsystem_of_units
#
#FACT = units.PhysCte(1,codeunits.UnitDic['s'])
#FACT = FACT.into(outputunits)
#
#val = val*FACT /1e3 # in Gyr
#return val,label
elif val == 'Age':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
val = nb.StellarAge(units=out_units.UnitTime)
return val,label
elif val == 'aFe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.metals[:,0]
return val,label
elif val == 'aMgFe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.metals[:,1]/(nb.metals[:,0]+1e-20)
return val,label
elif val == 'Fe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.Fe()
return val,label
elif val == 'Mg':
label = r'$\left[ \rm{Mg}/\rm{H} \right]$'
val = nb.Mg()
return val,label
elif val == 'MgFe':
label = r'$\left[ \rm{Mg}/\rm{Fe} \right]$'
val = nb.MgFe()
return val,label
else:
label = r'%s'%val
print "val = %s"%val
exec("val = %s"%val)
return val, label
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
# read files
for file in files:
nb = Nbody(file,ftype=opt.ftype)
################
# units
################
# define local units
unit_params = pt.do_units_options(opt)
nb.set_local_system_of_units(params=unit_params)
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
################
# apply options
################
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_select_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
################
# some info
################
print "---------------------------------------------------------"
nb.localsystem_of_units.info()
nb.ComovingIntegrationInfo()
print "---------------------------------------------------------"
################################
# compute physical quantities
################################
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(r'$\rm{Stellar\,\,Fraction}$',fontsize=pt.labelfont)
pt.grid(False)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment