Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F103110238
pplot.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
Thu, Feb 27, 08:56
Size
13 KB
Mime Type
text/x-python
Expires
Sat, Mar 1, 08:56 (2 d)
Engine
blob
Format
Raw Data
Handle
24500988
Attached To
rPNBODY pNbody
pplot.py
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
(
"--y"
,
action
=
"store"
,
dest
=
"y"
,
type
=
"string"
,
default
=
'T'
,
help
=
"y value to plot"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--z"
,
action
=
"store"
,
dest
=
"z"
,
type
=
"string"
,
default
=
None
,
help
=
"z value to plot"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--legend"
,
action
=
"store_true"
,
dest
=
"legend"
,
default
=
False
,
help
=
"add a legend"
)
parser
.
add_option
(
"--colorbar"
,
action
=
"store_true"
,
dest
=
"colorbar"
,
default
=
False
,
help
=
"add a colorbar"
)
parser
.
add_option
(
"--nx"
,
action
=
"store"
,
type
=
"int"
,
dest
=
"nx"
,
default
=
64
,
help
=
"number of bins in x"
)
parser
.
add_option
(
"--ny"
,
action
=
"store"
,
type
=
"int"
,
dest
=
"ny"
,
default
=
64
,
help
=
"number of bins in y"
)
parser
.
add_option
(
"--nxh"
,
action
=
"store"
,
type
=
"int"
,
dest
=
"nxh"
,
default
=
0
,
help
=
"number of bins in x"
)
parser
.
add_option
(
"--nopoints"
,
action
=
"store_true"
,
dest
=
"nopoints"
,
default
=
False
,
help
=
"do not plot points"
)
parser
.
add_option
(
"--map"
,
action
=
"store_true"
,
dest
=
"map"
,
default
=
False
,
help
=
"plot map instead of points"
)
parser
.
add_option
(
"--accumulate"
,
action
=
"store_true"
,
dest
=
"accumulate"
,
default
=
False
,
help
=
"integrate histogramm"
)
parser
.
add_option
(
"--density"
,
action
=
"store_true"
,
dest
=
"density"
,
default
=
False
,
help
=
"compute density"
)
parser
.
add_option
(
"--data"
,
action
=
"store"
,
type
=
"string"
,
dest
=
"data"
,
default
=
None
,
help
=
"data"
)
parser
.
add_option
(
"--mc"
,
action
=
"store"
,
type
=
"int"
,
dest
=
"mc"
,
default
=
0
,
help
=
"number montecarlo point"
)
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]$'
val
=
nb
.
Rxyz
()
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]$'
val
=
nb
.
Rho
()
# unit conversion
fact_g
=
units
.
PhysCte
(
1
,
nb
.
localsystem_of_units
.
UnitDic
[
'kg'
])
fact_g
=
fact_g
.
into
(
units
.
cgs
)
fact_cm
=
units
.
PhysCte
(
1
,
nb
.
localsystem_of_units
.
UnitDic
[
'm'
])
fact_cm
=
fact_cm
.
into
(
units
.
cgs
)
# factor conversion CU->cgs
fact_density
=
fact_g
/
fact_cm
**
3
# g -> H atom
fact_density
=
fact_density
/
ctes
.
PROTONMASS
.
into
(
units
.
cgs
)
val
=
val
*
fact_density
return
val
,
label
elif
val
==
'logrho'
:
label
=
r'$\rm{log\,Density}\,\left[ \rm{atom/cm^3} \right]$'
val
=
nb
.
Rho
()
# unit conversion
fact_g
=
units
.
PhysCte
(
1
,
nb
.
localsystem_of_units
.
UnitDic
[
'kg'
])
fact_g
=
fact_g
.
into
(
units
.
cgs
)
fact_cm
=
units
.
PhysCte
(
1
,
nb
.
localsystem_of_units
.
UnitDic
[
'm'
])
fact_cm
=
fact_cm
.
into
(
units
.
cgs
)
# factor conversion CU->cgs
fact_density
=
fact_g
/
fact_cm
**
3
# g -> H atom
fact_density
=
fact_density
/
ctes
.
PROTONMASS
.
into
(
units
.
cgs
)
val
=
val
*
fact_density
val
=
log10
(
val
)
return
val
,
label
elif
val
==
'Age'
:
label
=
r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
val
=
nb
.
Age
()
return
val
,
label
elif
val
==
'age'
:
label
=
r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units
=
units
.
UnitSystem
(
'local'
,[
units
.
Unit_kpc
,
units
.
Unit_Msol
,
units
.
Unit_Gyr
,
units
.
Unit_K
])
val
=
nb
.
age
()
*
nb
.
localsystem_of_units
.
convertionFactorTo
(
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
#################################
def
ComputeHisto
(
x
,
m
,
v
,
xmin
,
xmax
,
nx
):
#################################
'''
Compute an histogram
'''
shape
=
(
nx
,)
# normalize x
x
=
(
x
-
xmin
)
/
(
xmax
-
xmin
)
x
=
x
.
astype
(
float32
)
m
=
m
.
astype
(
float32
)
v
=
v
.
astype
(
float32
)
# compute the map
h
=
mkmap1d
(
x
,
m
,
v
,
shape
)
.
astype
(
float
)
# compute x
x
=
linspace
(
xmin
,
xmax
,
nx
)
return
x
,
h
#################################
def
ComputeMeanHisto
(
x
,
m
,
v
,
xmin
,
xmax
,
nx
):
#################################
"""
Compute the mean histogram
"""
x
,
h0
=
ComputeHisto
(
x
,
m
,
ones
(
len
(
m
)),
xmin
,
xmax
,
nx
)
x
,
h1
=
ComputeHisto
(
x
,
m
,
v
,
xmin
,
xmax
,
nx
)
h1
=
where
(
h0
==
0
,
0
,
h1
)
h0
=
where
(
h0
==
0
,
1
,
h0
)
h
=
h1
/
h0
return
x
,
h
def
get_histo
(
nb
,
val1
,
val2
,
val1min
,
val1max
,
n
):
if
(
val2
==
'm'
)
or
(
val2
==
'Lum'
):
mode
=
'extensive'
else
:
mode
=
'intensive'
print
val2
sys
.
exit
()
val1
,
l1
=
get_value_and_label
(
nb
,
val1
)
val2
,
l2
=
get_value_and_label
(
nb
,
val2
)
if
(
mode
==
'extensive'
):
m
=
ones
(
nb
.
nbody
)
x
,
vx
=
ComputeHisto
(
val1
,
m
,
val2
,
val1min
,
val1max
,
n
)
else
:
m
=
nb
.
mass
x
,
vx
=
ComputeMeanHisto
(
val1
,
m
,
val2
,
val1min
,
val1max
,
n
)
xlabel
=
l1
ylabel
=
l2
return
x
,
vx
,
xlabel
,
ylabel
#######################################
# MakePlot
#######################################
def
MakePlot
(
files
,
opt
):
# some inits
colors
=
pt
.
Colors
(
n
=
len
(
files
))
datas
=
[]
# define local units
unit_params
=
pt
.
do_units_options
(
opt
)
# read files
for
file
in
files
:
nb
=
Nbody
(
file
,
ftype
=
opt
.
ftype
)
################
# units
################
# define local units
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
"---------------------------------------------------------"
################
# get values
################
if
not
opt
.
nopoints
:
x
,
xlabel
=
get_value_and_label
(
nb
,
opt
.
x
)
y
,
ylabel
=
get_value_and_label
(
nb
,
opt
.
y
)
if
opt
.
z
!=
None
:
z
,
zlabel
=
get_value_and_label
(
nb
,
opt
.
z
)
data
=
pt
.
DataPoints
(
x
,
y
,
color
=
z
,
label
=
file
,
tpe
=
'points'
)
else
:
data
=
pt
.
DataPoints
(
x
,
y
,
color
=
colors
.
get
(),
label
=
file
,
tpe
=
'points'
)
datas
.
append
(
data
)
################
# get mean values
################
if
opt
.
nxh
>
0
:
v
,
h
,
xlabel
,
ylabel
=
get_histo
(
nb
,
opt
.
x
,
opt
.
y
,
opt
.
xmin
,
opt
.
xmax
,
opt
.
nxh
)
if
not
opt
.
accumulate
:
x
=
v
y
=
h
else
:
x
=
v
y
=
add
.
accumulate
(
h
)
if
opt
.
density
:
V
=
4.
/
3.
*
pi
*
(
v
[
1
:])
**
3
-
4.
/
3.
*
pi
*
(
v
[:
-
1
])
**
3
y
=
y
[:
-
1
]
/
V
x
=
x
[:
-
1
]
data
=
pt
.
DataPoints
(
x
,
y
,
color
=
colors
.
get
(),
label
=
file
,
tpe
=
'both'
)
datas
.
append
(
data
)
norm
=
None
# now, plot
if
not
opt
.
map
:
for
d
in
datas
:
if
d
.
tpe
==
'points'
or
d
.
tpe
==
'both'
:
if
opt
.
log
!=
None
:
if
string
.
find
(
opt
.
log
,
'z'
)
!=-
1
:
#norm = pt.colors.LogNorm(opt.zmin,opt.zmax)
norm
=
pt
.
colors
.
LogNorm
()
else
:
norm
=
None
cmap
=
pt
.
GetColormap
(
'rainbow4'
,
revesed
=
True
)
pt
.
scatter
(
d
.
x
,
d
.
y
,
c
=
d
.
color
,
s
=
5
,
linewidths
=
0
,
marker
=
'o'
,
vmin
=
opt
.
zmin
,
vmax
=
opt
.
zmax
,
norm
=
norm
,
cmap
=
cmap
)
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
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
)
if
opt
.
map
:
# now, plot
for
d
in
datas
:
x
=
d
.
x
y
=
d
.
y
z
=
zeros
(
len
(
x
))
# set log
#if log!=None:
# if string.find(opt.log,'x') != -1:
# x = log10(x)
# if string.find(opt.log,'y') != -1:
# y = log10(y)
c
=
isreal
(
x
)
*
isreal
(
y
)
x
=
compress
(
c
,
x
)
y
=
compress
(
c
,
y
)
# re-scale between 0 and 1
xs
=
1
-
(
x
-
xmax
)
/
(
xmin
-
xmax
)
ys
=
1
-
(
y
-
ymax
)
/
(
ymin
-
ymax
)
pos
=
transpose
(
array
([
xs
,
ys
,
z
],
float32
))
m
=
ones
(
len
(
x
),
float32
)
data
=
mapping
.
mkmap2d
(
pos
,
m
,
m
,(
opt
.
nx
,
opt
.
ny
))
data
=
transpose
(
data
)
/
sum
(
ravel
(
data
))
#########################
# compute and draw level
#########################
rdata
=
ravel
(
data
)
rdata
=
rdata
/
sum
(
rdata
)
zmin
=
min
(
rdata
)
zmax
=
max
(
rdata
)
def
levfct
(
x
,
flev
):
da
=
where
(
rdata
<
x
,
0.
,
rdata
)
return
sum
(
da
)
-
flev
# here, we cut at 90% of the flux
levelmax
=
zmax
levelmin
=
optimize
.
bisect
(
levfct
,
a
=
zmin
,
b
=
zmax
,
args
=
(
0.90
,),
xtol
=
1e-10
,
maxiter
=
100
)
cmap
=
pt
.
GetColormap
(
'heat'
,
revesed
=
True
)
#im = pt.imshow(data, interpolation='bilinear', origin='lower',cmap=cmap, extent=(xmin,xmax,ymin,ymax),aspect='auto')
im
=
pt
.
imshow
(
data
,
interpolation
=
'bilinear'
,
origin
=
'lower'
,
cmap
=
cmap
,
extent
=
(
xmin
,
xmax
,
ymin
,
ymax
),
aspect
=
'auto'
,
vmin
=
levelmin
,
vmax
=
levelmax
)
#cset = pt.contour(data,array([levelmin]),origin='lower',linewidths=1,extent=(xmin,xmax,ymin,ymax))
#########################
# compute monte carlo
#########################
if
opt
.
mc
>
0
:
from
pNbody
import
montecarlolib
datamc
=
transpose
(
data
)
pos
=
montecarlolib
.
mc2d
(
datamc
.
astype
(
float32
),
opt
.
mc
,
random
.
random
()
*
1000
)
# transform into physical coord
x
=
pos
[:,
0
]
y
=
pos
[:,
1
]
x
=
x
*
(
xmax
-
xmin
)
+
xmin
y
=
y
*
(
ymax
-
ymin
)
+
ymin
#pt.scatter(x,y,marker='o',s=1)
xe
=
0.1
*
ones
(
len
(
y
))
ye
=
0.2
*
ones
(
len
(
y
))
pt
.
errorbar
(
x
,
y
,
xerr
=
xe
,
yerr
=
ye
,
fmt
=
None
,
ecolor
=
'g'
,
lw
=
0.1
)
# 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
)
if
opt
.
colorbar
:
#lev = array([-2,0,15])
#l_f = pt.LogFormatter(10, labelOnlyBase=False)
#cb = pt.colorbar(ticks=lev, format = l_f)
cb
=
pt
.
colorbar
()
cb
.
set_label
(
zlabel
)
########################################################################
# 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