Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F72775238
gextractall.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
Wed, Jul 17, 00:02
Size
11 KB
Mime Type
text/x-python
Expires
Fri, Jul 19, 00:02 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19101663
Attached To
rPNBODY pNbody
gextractall.py
View Options
#!/usr/bin/env python
'''
Extract and plot energy and mass values contained in the
output Gadget file called by default "energy.txt".
Yves Revaz
ven jun 9 10:43:59 CEST 2006
'''
from
numpy
import
*
from
pNbody
import
*
import
string
import
sys
import
os
import
copy
as
docopy
from
pNbody.libutil
import
histogram
from
pNbody
import
libgrid
from
optparse
import
OptionParser
from
Gtools
import
*
from
Gtools
import
io
import
Ptools
as
pt
def
parse_options
():
usage
=
"usage: %prog [options] file"
parser
=
OptionParser
(
usage
=
usage
)
parser
=
pt
.
add_postscript_options
(
parser
)
parser
=
pt
.
add_limits_options
(
parser
)
parser
=
pt
.
add_log_options
(
parser
)
parser
=
pt
.
add_ftype_options
(
parser
)
parser
=
pt
.
add_reduc_options
(
parser
)
#parser = pt.add_center_options(parser)
parser
=
pt
.
add_display_options
(
parser
)
#parser = pt.add_select_options(parser)
parser
=
pt
.
add_cmd_options
(
parser
)
parser
=
pt
.
add_info_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
=
64
,
help
=
"number of bins in t"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--eps"
,
action
=
"store"
,
dest
=
"eps"
,
type
=
"float"
,
default
=
0.28
,
help
=
"smoothing length"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--nmin"
,
action
=
"store"
,
dest
=
"nmin"
,
type
=
"float"
,
default
=
32
,
help
=
"min number of particles in a cell to accept value"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--G"
,
action
=
"store"
,
dest
=
"G"
,
type
=
"float"
,
default
=
1.
,
help
=
"Gravitational constant value"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--xmode"
,
action
=
"store"
,
dest
=
"xmode"
,
type
=
"float"
,
default
=
2.
,
help
=
"mode for X-Toomre parameter"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--ErrTolTheta"
,
action
=
"store"
,
dest
=
"ErrTolTheta"
,
type
=
"float"
,
default
=
0.5
,
help
=
"Error tolerance theta"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--AdaptativeSoftenning"
,
action
=
"store_true"
,
dest
=
"AdaptativeSoftenning"
,
default
=
False
,
help
=
"AdaptativeSoftenning"
)
parser
.
add_option
(
"--statfile"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"statfile"
,
default
=
'stat.h5'
,
help
=
"stat output file"
)
parser
.
add_option
(
"--disk"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"disk"
,
default
=
"('gas','disk','stars')"
,
help
=
"stat output file"
)
parser
.
add_option
(
"--components"
,
action
=
"store"
,
dest
=
"components"
,
type
=
"string"
,
default
=
"('gas','halo','disk','bulge','stars')"
,
help
=
"list of components"
,
metavar
=
" TUPLE"
)
parser
.
add_option
(
"--x"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"x"
,
default
=
'R'
,
help
=
"x"
)
parser
.
add_option
(
"--y"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"y"
,
default
=
'vct'
,
help
=
"y"
)
parser
.
add_option
(
"--mode"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"mode"
,
default
=
'all'
,
help
=
"mode"
)
(
options
,
args
)
=
parser
.
parse_args
()
if
len
(
args
)
==
0
:
print
"you must specify a filename"
sys
.
exit
(
0
)
files
=
args
return
files
,
options
#######################################
# MakePlot
#######################################
def
MakePlot
(
dirs
,
opt
):
# some inits
palette
=
pt
.
GetPalette
()
colors
=
pt
.
SetColorsForFiles
(
files
,
palette
)
labels
=
[]
#######################################
# deal with mode
#######################################
if
opt
.
mode
==
'all'
:
opt
.
ComputePotential
=
True
opt
.
ComputeSurfaceDensity
=
True
opt
.
ComputeDispertions
=
True
opt
.
ComputeStability
=
True
elif
opt
.
mode
==
'Sden'
:
opt
.
ComputePotential
=
False
opt
.
ComputeSurfaceDensity
=
True
opt
.
ComputeDispertions
=
False
opt
.
ComputeStability
=
False
elif
opt
.
mode
==
'velocities'
:
opt
.
ComputePotential
=
False
opt
.
ComputeSurfaceDensity
=
False
opt
.
ComputeDispertions
=
True
opt
.
ComputeStability
=
False
elif
opt
.
mode
==
'stability'
:
opt
.
ComputePotential
=
True
opt
.
ComputeSurfaceDensity
=
True
opt
.
ComputeDispertions
=
True
opt
.
ComputeStability
=
True
#######################################
# LOOP
#######################################
# read files
for
file
in
files
:
######################################
# open file and apply option
######################################
nb
=
Nbody
(
file
,
ftype
=
opt
.
ftype
)
nb
=
pt
.
do_reduc_options
(
nb
,
opt
)
nb
=
pt
.
do_cmd_options
(
nb
,
opt
)
###############
# create total
###############
nbt
=
docopy
.
deepcopy
(
nb
)
###############
# create disk
###############
nbd
=
docopy
.
deepcopy
(
nb
)
if
opt
.
disk
!=
None
:
opt
.
disk
=
eval
(
opt
.
disk
,
dict
(
__builtins__
=
None
))
print
"disk = "
,
opt
.
disk
nbd
=
nbd
.
select
(
opt
.
disk
)
######################################################
# here, we should loop over component
if
opt
.
components
!=
None
:
opt
.
components
=
eval
(
opt
.
components
,
dict
(
__builtins__
=
None
))
print
"components = "
,
opt
.
components
for
component
in
opt
.
components
:
print
"-----------------------------"
print
"-- component = "
,
component
print
"-----------------------------"
nb
=
nbt
.
select
(
component
)
opt
.
statfile
=
"
%s
.h5"
%
(
component
)
if
nb
.
nbody
>
0
:
###############
# other info
###############
#nb = pt.do_center_options(nb,opt)
nb
=
pt
.
do_info_options
(
nb
,
opt
)
nb
=
pt
.
do_display_options
(
nb
,
opt
)
######################################
# computes values
######################################
# create cylindrical rt grid
rc
=
1.
g
=
lambda
r
:
log
(
r
/
rc
+
1.
)
gm
=
lambda
r
:
rc
*
(
exp
(
r
)
-
1.
)
g
=
None
;
gm
=
None
Grt
=
libgrid
.
Cylindrical_2drt_Grid
(
rmin
=
0
,
rmax
=
opt
.
Rmax
,
nr
=
opt
.
nR
,
nt
=
opt
.
nt
,
z
=
0
,
g
=
g
,
gm
=
gm
)
# build the tree
if
opt
.
ComputePotential
:
print
"ComputeTree"
nbt
.
getTree
(
force_computation
=
True
,
ErrTolTheta
=
opt
.
ErrTolTheta
)
nb
.
getTree
(
force_computation
=
True
,
ErrTolTheta
=
opt
.
ErrTolTheta
)
# radius
R
,
t
=
Grt
.
get_rt
()
stats
=
{}
stats
[
'nR'
]
=
opt
.
nR
stats
[
'nt'
]
=
opt
.
nt
stats
[
'Rmax'
]
=
opt
.
Rmax
stats
[
'eps'
]
=
opt
.
eps
stats
[
'ErrTolTheta'
]
=
opt
.
ErrTolTheta
stats
[
'AdaptativeSoftenning'
]
=
opt
.
AdaptativeSoftenning
stats
[
'xmode'
]
=
opt
.
xmode
stats
[
'nmin'
]
=
opt
.
nmin
stats
[
'G'
]
=
opt
.
G
stats
[
'R'
]
=
R
stats
[
't'
]
=
t
###################################
# Surface density
###################################
if
opt
.
ComputeSurfaceDensity
:
print
"ComputeSurfaceDensity"
Sden
=
Grt
.
get_SurfaceDensityMap
(
nb
)
Sden
=
sum
(
Sden
,
1
)
/
opt
.
nt
Sdend
=
Grt
.
get_SurfaceDensityMap
(
nbd
)
Sdend
=
sum
(
Sdend
,
1
)
/
opt
.
nt
stats
[
'Sden'
]
=
Sden
stats
[
'Sdend'
]
=
Sdend
###################################
# rotation curve and frequencies
###################################
if
opt
.
ComputePotential
:
print
"ComputePotential"
# radial acceleration (total)
Accx
,
Accy
,
Accz
=
Grt
.
get_AccelerationMap
(
nbt
,
eps
=
opt
.
eps
,
UseTree
=
True
,
AdaptativeSoftenning
=
opt
.
AdaptativeSoftenning
)
Ar
=
sqrt
(
Accx
**
2
+
Accy
**
2
)
Ar
=
sum
(
Ar
,
1
)
/
opt
.
nt
dPhit
=
Ar
d2Phit
=
libgrid
.
get_First_Derivative
(
dPhit
,
R
)
kappat
=
libdisk
.
Kappa
(
R
,
dPhit
,
d2Phit
)
omegat
=
libdisk
.
Omega
(
R
,
dPhit
)
vct
=
libdisk
.
Vcirc
(
R
,
dPhit
)
#nut = libdisk.Nu(z,Phit)
# 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
d2Phi
=
libgrid
.
get_First_Derivative
(
dPhi
,
R
)
kappa
=
libdisk
.
Kappa
(
R
,
dPhi
,
d2Phi
)
omega
=
libdisk
.
Omega
(
R
,
dPhi
)
vc
=
libdisk
.
Vcirc
(
R
,
dPhi
)
#nu = libdisk.Nu(z,Phi)
#stats['Phi'] = Phi
stats
[
'dPhi'
]
=
dPhi
stats
[
'd2Phi'
]
=
d2Phi
stats
[
'kappa'
]
=
kappa
stats
[
'omega'
]
=
omega
#stats['nu'] = nu
#stats['Phit'] = Phit
stats
[
'dPhit'
]
=
dPhit
stats
[
'd2Phit'
]
=
d2Phit
stats
[
'kappat'
]
=
kappat
stats
[
'omegat'
]
=
omegat
#stats['nut'] = nut
stats
[
'vct'
]
=
vct
stats
[
'vc'
]
=
vc
if
opt
.
ComputeDispertions
:
print
"ComputeDispertions"
###################################
# number of points per cell
###################################
nn
=
Grt
.
get_NumberMap
(
nb
)
nm
=
sum
(
nn
,
1
)
/
opt
.
nt
###################################
# mean velocities
###################################
vr
=
Grt
.
get_MeanValMap
(
nb
,
nb
.
Vr
())
vr
=
get1dMeanFrom2dMap
(
vr
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
vm
=
Grt
.
get_MeanValMap
(
nb
,
nb
.
Vt
())
vm
=
get1dMeanFrom2dMap
(
vm
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
vz
=
Grt
.
get_MeanValMap
(
nb
,
nb
.
Vz
())
vz
=
get1dMeanFrom2dMap
(
vz
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
###################################
# velocity dispertions
###################################
sr
=
Grt
.
get_SigmaValMap
(
nb
,
nb
.
Vr
())
sr
=
get1dMeanFrom2dMap
(
sr
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
sp
=
Grt
.
get_SigmaValMap
(
nb
,
nb
.
Vt
())
sp
=
get1dMeanFrom2dMap
(
sp
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
sz
=
Grt
.
get_SigmaValMap
(
nb
,
nb
.
Vz
())
sz
=
get1dMeanFrom2dMap
(
sz
,
nn
,
nmin
=
opt
.
nmin
,
axis
=
1
)
stats
[
'nn'
]
=
nn
stats
[
'nm'
]
=
nm
stats
[
'vr'
]
=
vr
stats
[
'vm'
]
=
vm
stats
[
'vz'
]
=
vz
stats
[
'sr'
]
=
sr
stats
[
'sp'
]
=
sp
stats
[
'sz'
]
=
sz
###################################
# Q,X,A
###################################
if
opt
.
ComputeStability
:
print
"ComputeStability"
Q
=
libdisk
.
QToomre
(
opt
.
G
,
R
,
sr
,
kappat
,
Sdend
)
X
=
libdisk
.
XToomre
(
opt
.
G
,
R
,
kappat
,
Sdend
,
opt
.
xmode
)
A
=
libdisk
.
AAraki
(
sr
,
sz
)
stats
[
'Q'
]
=
Q
stats
[
'A'
]
=
A
stats
[
'X'
]
=
X
###################################
# save hdf
##################s#################
pt
.
io
.
write_hdf5
(
opt
.
statfile
,
stats
)
###################################
# now, plot
###################################
x
=
stats
[
opt
.
x
]
y
=
stats
[
opt
.
y
]
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimits
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
x
,
y
,
opt
.
log
)
pt
.
plot
(
x
,
y
)
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
opt
.
log
)
pt
.
show
()
if
__name__
==
'__main__'
:
files
,
opt
=
parse_options
()
pt
.
InitPlot
(
files
,
opt
)
#pt.figure(figsize=(8*2,6*2))
#pt.figure(dpi=10)
pt
.
pcolors
#fig = pt.gcf()
#fig.subplots_adjust(left=0.1)
#fig.subplots_adjust(right=1)
#fig.subplots_adjust(bottom=0.12)
#fig.subplots_adjust(top=0.95)
#fig.subplots_adjust(wspace=0.25)
#fig.subplots_adjust(hspace=0.02)
MakePlot
(
files
,
opt
)
pt
.
EndPlot
(
files
,
opt
)
Event Timeline
Log In to Comment