Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72770550
plot_rotation_curve.py
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
Tue, Jul 16, 23:23
Size
5 KB
Mime Type
text/x-python
Expires
Thu, Jul 18, 23:23 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19101768
Attached To
rPNBODY pNbody
plot_rotation_curve.py
View Options
#!/usr/bin/env python
'''
Plot value computed in a spherical grid, as a function of the radius
Yves Revaz
Wed Jan 9 09:29:46 CET 2008
'''
import
Ptools
as
pt
from
pNbody
import
*
from
pNbody
import
myNumeric
from
pNbody
import
libdisk
from
pNbody
import
libgrid
from
pNbody
import
profiles
from
pNbody
import
cosmo
import
string
import
sys
import
os
from
numpy
import
*
from
optparse
import
OptionParser
from
scipy.optimize
import
leastsq
from
scipy.optimize
import
bisection
#from mpfit import *
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
)
parser
=
pt
.
add_legend_options
(
parser
)
parser
.
add_option
(
"--Rmax"
,
action
=
"store"
,
dest
=
"Rmax"
,
type
=
"float"
,
default
=
50.
,
help
=
"max radius of bins"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--nR"
,
action
=
"store"
,
dest
=
"nR"
,
type
=
"int"
,
default
=
32
,
help
=
"number of bins in r"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--nt"
,
action
=
"store"
,
dest
=
"nt"
,
type
=
"int"
,
default
=
32
,
help
=
"number of bins in r"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--ErrTolTheta"
,
action
=
"store"
,
dest
=
"ErrTolTheta"
,
type
=
"float"
,
default
=
0.5
,
help
=
"Error tolerance theta"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--fct"
,
action
=
"store"
,
dest
=
"fct"
,
type
=
"string"
,
default
=
None
,
help
=
"transformation function ex : 'r**3' should be of the form M(r)"
,
metavar
=
" STR"
)
parser
.
add_option
(
"--fctm"
,
action
=
"store"
,
dest
=
"fctm"
,
type
=
"string"
,
default
=
None
,
help
=
"inverse transformation function ex : 'r**(1/3.)'"
,
metavar
=
" STR"
)
parser
.
add_option
(
"--AdaptativeSoftenning"
,
action
=
"store_true"
,
dest
=
"AdaptativeSoftenning"
,
default
=
False
,
help
=
"AdaptativeSoftenning"
)
parser
.
add_option
(
"--eps"
,
action
=
"store"
,
dest
=
"eps"
,
type
=
"float"
,
default
=
0.28
,
help
=
"smoothing length"
,
metavar
=
" FLOAT"
)
(
options
,
args
)
=
parser
.
parse_args
()
if
len
(
args
)
==
0
:
print
"you must specify a filename"
sys
.
exit
(
0
)
files
=
args
return
files
,
options
#######################################
# units
#######################################
# gear system of units
params
=
{}
params
[
'UnitVelocity_in_cm_per_s'
]
=
20725573.785998672
params
[
'UnitMass_in_g'
]
=
1.989e+43
params
[
'UnitLength_in_cm'
]
=
3.085e+21
gear_units
=
units
.
Set_SystemUnits_From_Params
(
params
)
out_units
=
units
.
UnitSystem
(
'local'
,[
units
.
Unit_km
,
units
.
Unit_Msol
,
units
.
Unit_s
,
units
.
Unit_K
])
tokms
=
gear_units
.
convertionFactorTo
(
out_units
.
UnitVelocity
)
#######################################
# MakePlot
#######################################
def
MakePlot
(
files
,
opt
):
# some inits
colors
=
pt
.
Colors
(
n
=
len
(
files
))
datas
=
[]
velocity_factor
=
207.
# read files
for
i
,
file
in
enumerate
(
files
):
nb
=
Nbody
(
file
,
ftype
=
opt
.
ftype
)
# 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
)
#######################
# set the grid scaling
#######################
if
opt
.
fct
!=
None
:
print
"g =lambda r:(
%s
)"
%
opt
.
fct
exec
(
"g =lambda r:(
%s
)"
%
opt
.
fct
)
else
:
g
=
None
if
opt
.
fctm
!=
None
:
print
"gm =lambda r:(
%s
)"
%
opt
.
fctm
exec
(
"gm =lambda r:(
%s
)"
%
opt
.
fctm
)
else
:
gm
=
None
##############################
# compute the tree
##############################
nb
.
getTree
(
force_computation
=
True
,
ErrTolTheta
=
opt
.
ErrTolTheta
)
##############################
# create cylindrical rt grid
##############################
Grt
=
libgrid
.
Cylindrical_2drt_Grid
(
rmin
=
0
,
rmax
=
opt
.
Rmax
,
nr
=
opt
.
nR
,
nt
=
opt
.
nt
,
z
=
0
,
g
=
g
,
gm
=
gm
)
# radial acceleration (selection)
Accx
,
Accy
,
Accz
=
Grt
.
get_AccelerationMap
(
nb
,
eps
=
opt
.
eps
,
UseTree
=
True
,
AdaptativeSoftenning
=
opt
.
AdaptativeSoftenning
)
Ar
=
sqrt
(
Accx
**
2
+
Accy
**
2
)
Ar
=
sum
(
Ar
,
1
)
/
opt
.
nt
dPhi
=
Ar
# radius
R
,
t
=
Grt
.
get_rt
()
# velocity
Vc
=
libdisk
.
Vcirc
(
R
,
dPhi
)
*
tokms
datas
.
append
(
pt
.
DataPoints
(
R
,
Vc
,
color
=
colors
.
next
(),
label
=
i
))
# plot
for
d
in
datas
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
# set limits and draw axis
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
opt
.
log
)
xlabel
=
r"$\rm{Radius}\,[\rm{kpc}]$"
ylabel
=
r"$\rm{Velocity}\,[\rm{km/s}]$"
# plot axis
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
opt
.
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
pt
.
grid
(
False
)
if
opt
.
legend
:
pt
.
LegendFromDataPoints
(
datas
)
########################################################################
# 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