Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F105423782
pychem_SSP_discret_evolution2
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
Mon, Mar 17, 03:05
Size
27 KB
Mime Type
text/x-python
Expires
Wed, Mar 19, 03:05 (2 d)
Engine
blob
Format
Raw Data
Handle
24963967
Attached To
rGEAR Gear
pychem_SSP_discret_evolution2
View Options
#!/usr/bin/env python
from
optparse
import
OptionParser
from
PyChem
import
chemistry
from
numpy
import
*
import
sys
import
string
import
Ptools
as
pt
def
parse_options
():
usage
=
"usage: %prog [options] file"
parser
=
OptionParser
(
usage
=
usage
)
parser
=
pt
.
add_limits_options
(
parser
)
parser
=
pt
.
add_log_options
(
parser
)
parser
=
pt
.
add_postscript_options
(
parser
)
parser
.
add_option
(
"--x"
,
action
=
"store"
,
dest
=
"x"
,
type
=
"string"
,
default
=
'Fe'
,
help
=
"x value to plot"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--y"
,
action
=
"store"
,
dest
=
"y"
,
type
=
"string"
,
default
=
'Mg'
,
help
=
"y value to plot"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--dt"
,
action
=
"store"
,
dest
=
"dt"
,
type
=
"float"
,
default
=
0.1
,
help
=
"dt"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--mstar"
,
action
=
"store"
,
dest
=
"mstar"
,
type
=
"float"
,
default
=
1e5
,
help
=
"initial mass of the SSP in solar mass"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--tstar"
,
action
=
"store"
,
dest
=
"tstar"
,
type
=
"float"
,
default
=
0
,
help
=
"formation time of the SSP"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"-o"
,
action
=
"store"
,
dest
=
"obs"
,
type
=
"string"
,
default
=
'Y'
,
help
=
"observable to plot"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--timeunit"
,
action
=
"store"
,
dest
=
"timeunit"
,
type
=
"string"
,
default
=
None
,
help
=
"unit of time"
,
metavar
=
" STRING"
)
parser
.
add_option
(
"--NumberOfTables"
,
action
=
"store"
,
dest
=
"NumberOfTables"
,
type
=
"int"
,
default
=
1
,
help
=
"NumberOfTables"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--DefaultTable"
,
action
=
"store"
,
dest
=
"DefaultTable"
,
type
=
"int"
,
default
=
0
,
help
=
"DefaultTable"
,
metavar
=
" INT"
)
parser
.
add_option
(
"--discsn"
,
action
=
"store_true"
,
dest
=
"discsn"
,
default
=
False
,
help
=
"discretize sn"
,
metavar
=
" FLOAT"
)
parser
.
add_option
(
"--seed"
,
action
=
"store"
,
dest
=
"seed"
,
type
=
"int"
,
default
=
None
,
help
=
"initial random seed"
,
metavar
=
" INT"
)
(
options
,
args
)
=
parser
.
parse_args
()
pt
.
check_files_number
(
args
)
files
=
args
return
files
,
options
########################################################################
# M A I N
########################################################################
SOLAR_MASS
=
1.989e33
UnitLength_in_cm
=
3.085e+21
UnitMass_in_g
=
1.989e+43
UnitVelocity_in_cm_per_s
=
20725573.785998672
UnitTime_in_s
=
148849920000000.0
feedback_energy
=
10
**
51
/
UnitMass_in_g
/
(
UnitVelocity_in_cm_per_s
)
**
2
# parse options
files
,
opt
=
parse_options
()
if
opt
.
seed
!=
None
:
random
.
seed
(
opt
.
seed
)
# units
opt
.
ux
=
1.0
if
opt
.
timeunit
==
None
:
opt
.
ux
=
1.0
elif
opt
.
timeunit
==
'Myr'
:
opt
.
ux
=
UnitTime_in_s
/
(
3600
*
24
*
365
)
/
1e6
elif
opt
.
timeunit
==
'Gyr'
:
opt
.
ux
=
UnitTime_in_s
/
(
3600
*
24
*
365
)
/
1e9
# to plot
opt
.
obs
=
string
.
split
(
opt
.
obs
,
','
)
# init chimie
chemistry
.
init_chimie
(
files
[
0
],
opt
.
NumberOfTables
,
opt
.
DefaultTable
)
elts
=
chemistry
.
get_elts_labels
()
nelts
=
chemistry
.
get_nelts
()
print
elts
# some parameters
dt
=
opt
.
dt
tmax
=
opt
.
xmax
/
opt
.
ux
t_star
=
opt
.
tstar
# formation time of the SSP
M0
=
opt
.
mstar
*
SOLAR_MASS
/
UnitMass_in_g
# gas mass of the SSP (in code unit)
Mstacont
=
M0
# star mass
Mstaj
=
M0
Mstajcont
=
M0
Mgascont
=
0.015
*
M0
# gas mass
Mgasj
=
0.015
*
M0
# gas mass
Mgasjcont
=
0.015
*
M0
# gas mass
# stellar initial element abondance
metals
=
zeros
(
chemistry
.
get_nelts
(),
float
)
#metals[0] = 0.001771 *1
#metals[1] = 0.00091245 *1
#metals[2] = 0.0 *1
#metals[3] = 0.02 *1 # metal abondance (=metallicity z)
# gas initial element abondance
metalgcont
=
zeros
(
chemistry
.
get_nelts
(),
float
)
metalgj
=
zeros
(
chemistry
.
get_nelts
(),
float
)
metalgjcont
=
zeros
(
chemistry
.
get_nelts
(),
float
)
# minimum live time of a star
minlivetime
=
chemistry
.
star_lifetime
(
metals
[
-
1
],
chemistry
.
get_Mmax
())
maxlivetime
=
chemistry
.
star_lifetime
(
metals
[
-
1
],
chemistry
.
get_Mmin
())
print
minlivetime
print
maxlivetime
# times
times
=
arange
(
dt
,
tmax
,
dt
)
.
astype
(
float
)
nt
=
len
(
times
)
# output arrays
TIMES
=
zeros
(
nt
,
float
)
# time
NSNIas
=
zeros
(
nt
,
float
)
# number of SNIa
NSNIIs
=
zeros
(
nt
,
float
)
# number of SNII
NSNIaconts
=
zeros
(
nt
,
float
)
# number of SNIa (continuous value)
NSNIIconts
=
zeros
(
nt
,
float
)
# number of SNII (continuous value)
NDYINs
=
zeros
(
nt
,
float
)
# number of dying stars
NDYINconts
=
zeros
(
nt
,
float
)
# number of dying stars
Es
=
zeros
(
nt
,
float
)
# injected energy
Mecont
=
zeros
(
nt
,
float
)
# ejected mass (continuous)
Mej
=
zeros
(
nt
,
float
)
# ejected mass (computed from NSN)
Mejcont
=
zeros
(
nt
,
float
)
# ejected mass (computed from NSN,continuous)
Mstaconts
=
zeros
(
nt
,
float
)
# star mass (continuous)
Mstajs
=
zeros
(
nt
,
float
)
# star mass
Mstajconts
=
zeros
(
nt
,
float
)
# star mass
Mgasconts
=
zeros
(
nt
,
float
)
# gas mass
Mgasjs
=
zeros
(
nt
,
float
)
# gas mass
Mgasjconts
=
zeros
(
nt
,
float
)
# gas mass
# elements to follow
# ejected
dataecont
=
{}
dataej
=
{}
dataejcont
=
{}
for
elt
in
elts
:
dataecont
[
elt
]
=
zeros
(
nt
,
float
)
dataej
[
elt
]
=
zeros
(
nt
,
float
)
dataejcont
[
elt
]
=
zeros
(
nt
,
float
)
# gas
datagcont
=
{}
datagj
=
{}
datagjcont
=
{}
for
elt
in
elts
:
datagcont
[
elt
]
=
zeros
(
nt
,
float
)
datagj
[
elt
]
=
zeros
(
nt
,
float
)
datagjcont
[
elt
]
=
zeros
(
nt
,
float
)
###############################################
# main loop
###############################################
for
ti
,
time
in
enumerate
(
times
):
tpdt
=
time
+
dt
if
(
tpdt
-
t_star
>=
minlivetime
)
and
(
tpdt
-
t_star
<
maxlivetime
):
m1
=
chemistry
.
star_mass_from_age
(
metals
[
-
1
],
tpdt
-
t_star
)
if
(
time
-
t_star
>=
minlivetime
):
m2
=
chemistry
.
star_mass_from_age
(
metals
[
-
1
],
time
-
t_star
)
else
:
m2
=
chemistry
.
get_Mmax
()
if
(
m1
>
m2
):
m1
=
m2
raise
"m1>m2 !!!"
# number of SNIa
NSNIa
=
chemistry
.
SNIa_rate
(
m1
,
m2
)
*
M0
# number of SNII
NSNII
=
chemistry
.
SNII_rate
(
m1
,
m2
)
*
M0
# number of dying stars
NDYIN
=
chemistry
.
DYIN_rate
(
m1
,
m2
)
*
M0
#########################
# discrete stuff
NSNIIcont
=
NSNII
NSNIacont
=
NSNIa
NDYINcont
=
NDYIN
if
opt
.
discsn
:
fNSNIa
=
NSNIa
-
floor
(
NSNIa
)
# fraction
NSNIa
=
floor
(
NSNIa
)
# integer
if
random
.
random
()
<
fNSNIa
:
NSNIa
=
NSNIa
+
1
if
opt
.
discsn
:
fNSNII
=
NSNII
-
floor
(
NSNII
)
# fraction
NSNII
=
floor
(
NSNII
)
# integer
if
random
.
random
()
<
fNSNII
:
NSNII
=
NSNII
+
1
if
opt
.
discsn
:
fNDYIN
=
NDYIN
-
floor
(
NDYIN
)
# fraction
NDYIN
=
floor
(
NDYIN
)
# integer
if
random
.
random
()
<
fNDYIN
:
NDYIN
=
NDYIN
+
1
#
#########################
##################################################
# mass and yields ejection
##################################################
# standard way (continuous explosion)
EjectedMass
=
chemistry
.
Total_mass_ejection
(
m1
,
m2
,
M0
,
metals
)
TotalEjectedEltMasscont
=
EjectedMass
[
2
:]
TotalEjectedGasMasscont
=
EjectedMass
[
0
]
# new way (discrete explosion)
EjectedMassj
=
chemistry
.
Total_single_mass_ejection
(
0.5
*
(
m1
+
m2
),
metals
,
float
(
NSNII
),
float
(
NSNIa
),
float
(
NDYIN
))
TotalEjectedEltMassj
=
EjectedMassj
[
2
:]
TotalEjectedGasMassj
=
EjectedMassj
[
0
]
# new way (continuous explosion)
EjectedMassjcont
=
chemistry
.
Total_single_mass_ejection
(
0.5
*
(
m1
+
m2
),
metals
,
NSNIIcont
,
NSNIacont
,
NDYINcont
)
TotalEjectedEltMassjcont
=
EjectedMassjcont
[
2
:]
TotalEjectedGasMassjcont
=
EjectedMassjcont
[
0
]
# other check
#SNII_EjectedMass = chemistry.SNII_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNII
#SNIa_EjectedMass = chemistry.SNIa_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNIa
#TotalEjectedEltMassjcont = SNII_EjectedMass[2:] + SNIa_EjectedMass[2:]
#TotalEjectedGasMassjcont = SNII_EjectedMass[0] + SNIa_EjectedMass[0]
'''
# all dying stars (= SNII + lighter masses)
DYIN_TotalEjectedGasMassj = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIa # !!! need to add SNIa
DYIN_TotalEjectedGasMassjcont = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIacont # !!! need to add SNIa
'''
##################################################
# energy
##################################################
# energy injected
TotalInjectedEnergy
=
feedback_energy
*
(
NSNIa
+
NSNII
)
# mass fraction of ejected elements
emetalcont
=
TotalEjectedEltMasscont
/
TotalEjectedGasMasscont
emetalj
=
TotalEjectedEltMassj
/
TotalEjectedGasMassj
emetaljcont
=
TotalEjectedEltMassjcont
/
TotalEjectedGasMassjcont
# metal enrichment
metalgcont
=
((
metalgcont
*
Mgascont
)
+
TotalEjectedEltMasscont
)
metalgj
=
((
metalgj
*
Mgasj
)
+
TotalEjectedEltMassj
)
metalgjcont
=
((
metalgjcont
*
Mgasjcont
)
+
TotalEjectedEltMassjcont
)
# Total Mass
Mstacont
=
Mstacont
-
TotalEjectedGasMasscont
Mstaj
=
Mstaj
-
TotalEjectedGasMassj
Mstajcont
=
Mstajcont
-
TotalEjectedGasMassjcont
# Gas Mass
Mgascont
=
Mgascont
+
TotalEjectedGasMasscont
Mgasj
=
Mgasj
+
TotalEjectedGasMassj
Mgasjcont
=
Mgasjcont
+
TotalEjectedGasMassjcont
# metal enrichment
metalgcont
=
metalgcont
/
Mgascont
metalgj
=
metalgj
/
Mgasj
metalgjcont
=
metalgjcont
/
Mgasjcont
# record output
TIMES
[
ti
]
=
time
NSNIas
[
ti
]
=
NSNIa
NSNIIs
[
ti
]
=
NSNII
NSNIaconts
[
ti
]
=
NSNIacont
NSNIIconts
[
ti
]
=
NSNIIcont
NDYINs
[
ti
]
=
NDYIN
NDYINconts
[
ti
]
=
NDYINcont
Es
[
ti
]
=
TotalInjectedEnergy
Mecont
[
ti
]
=
TotalEjectedGasMasscont
Mej
[
ti
]
=
TotalEjectedGasMassj
Mejcont
[
ti
]
=
TotalEjectedGasMassjcont
Mstaconts
[
ti
]
=
Mstacont
# stellar mass
Mstajs
[
ti
]
=
Mstaj
# stellar mass
Mstajconts
[
ti
]
=
Mstajcont
# stellar mass
Mgasconts
[
ti
]
=
Mgascont
Mgasjs
[
ti
]
=
Mgasj
Mgasjconts
[
ti
]
=
Mgasjcont
for
j
,
elt
in
enumerate
(
elts
):
dataecont
[
elt
][
ti
]
=
emetalcont
[
j
]
datagcont
[
elt
][
ti
]
=
metalgcont
[
j
]
dataej
[
elt
][
ti
]
=
emetalj
[
j
]
datagj
[
elt
][
ti
]
=
metalgj
[
j
]
dataejcont
[
elt
][
ti
]
=
emetaljcont
[
j
]
datagjcont
[
elt
][
ti
]
=
metalgjcont
[
j
]
# remove time=0 values
c
=
(
TIMES
!=
0
)
TIMES
=
compress
(
c
,
TIMES
)
*
opt
.
ux
NSNIas
=
compress
(
c
,
NSNIas
)
NSNIIs
=
compress
(
c
,
NSNIIs
)
NSNIaconts
=
compress
(
c
,
NSNIaconts
)
NSNIIconts
=
compress
(
c
,
NSNIIconts
)
NDYINs
=
compress
(
c
,
NDYINs
)
NDYINconts
=
compress
(
c
,
NDYINconts
)
Es
=
compress
(
c
,
Es
)
Mecont
=
compress
(
c
,
Mecont
)
Mej
=
compress
(
c
,
Mej
)
Mejcont
=
compress
(
c
,
Mejcont
)
Mstaconts
=
compress
(
c
,
Mstaconts
)
Mstajs
=
compress
(
c
,
Mstajs
)
Mstajconts
=
compress
(
c
,
Mstajconts
)
Mgasconts
=
compress
(
c
,
Mgasconts
)
Mgasjs
=
compress
(
c
,
Mgasjs
)
Mgasjconts
=
compress
(
c
,
Mgasjconts
)
for
elt
in
elts
:
dataecont
[
elt
]
=
compress
(
c
,
dataecont
[
elt
])
datagcont
[
elt
]
=
compress
(
c
,
datagcont
[
elt
])
dataej
[
elt
]
=
compress
(
c
,
dataej
[
elt
])
datagj
[
elt
]
=
compress
(
c
,
datagj
[
elt
])
dataejcont
[
elt
]
=
compress
(
c
,
dataejcont
[
elt
])
datagjcont
[
elt
]
=
compress
(
c
,
datagjcont
[
elt
])
# get Solar Abundances
SolarAbun
=
chemistry
.
get_elts_SolarMassAbundances
()
Xsconts
=
{}
Xsjs
=
{}
Xsjconts
=
{}
Xgconts
=
{}
Xgjs
=
{}
Xgjconts
=
{}
for
elt
in
elts
:
Xsconts
[
elt
]
=
log10
(
dataecont
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xgconts
[
elt
]
=
log10
(
datagcont
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xsjs
[
elt
]
=
log10
(
dataej
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xgjs
[
elt
]
=
log10
(
datagj
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xsjconts
[
elt
]
=
log10
(
dataejcont
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xgjconts
[
elt
]
=
log10
(
datagjcont
[
elt
]
/
SolarAbun
[
elt
]
+
1.0e-10
)
Xscont
=
Xsconts
[
opt
.
x
]
Yscont
=
Xsconts
[
opt
.
y
]
Xgcont
=
Xgconts
[
opt
.
x
]
Ygcont
=
Xgconts
[
opt
.
y
]
Xsj
=
Xsjs
[
opt
.
x
]
Ysj
=
Xsjs
[
opt
.
y
]
Xgj
=
Xgjs
[
opt
.
x
]
Ygj
=
Xgjs
[
opt
.
y
]
Xsjcont
=
Xsjconts
[
opt
.
x
]
Ysjcont
=
Xsjconts
[
opt
.
y
]
Xgjcont
=
Xgjconts
[
opt
.
x
]
Ygjcont
=
Xgjconts
[
opt
.
y
]
#########################################################
# plot
#########################################################
pt
.
InitPlot
([],
opt
)
pt
.
pcolors
fig
=
pt
.
gcf
()
#fig.subplots_adjust(left=0.02)
#fig.subplots_adjust(right=0.99)
#fig.subplots_adjust(bottom=0.1)
#fig.subplots_adjust(top=0.99)
#fig.subplots_adjust(wspace=0.0)
fig
.
subplots_adjust
(
hspace
=
0.25
)
np
=
len
(
opt
.
obs
)
def
PlotYields
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{[
%s
/H],[
%s
/H],[
%s
/
%s
]}$"
%
(
opt
.
x
,
opt
.
y
,
opt
.
y
,
opt
.
x
)
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Xscont
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Yscont
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Yscont
-
Xscont
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Xsj
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ysj
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ysj
-
Xsj
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'-'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Xsjcont
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
)
#datas.append(data)
data
=
pt
.
DataPoints
(
TIMES
,
Ysjcont
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
)
#datas.append(data)
data
=
pt
.
DataPoints
(
TIMES
,
Ysjcont
-
Xsjcont
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'--'
)
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
,
ls
=
d
.
linestyle
)
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
opt
.
log
)
ymin
=
-
2
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
opt
.
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotGasElts
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{gas\,[
%s
/H],[
%s
/H],[
%s
/
%s
]}$"
%
(
opt
.
x
,
opt
.
y
,
opt
.
y
,
opt
.
x
)
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Xgcont
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ygcont
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ygcont
-
Xgcont
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Xgj
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ygj
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
,
linestyle
=
'--'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ygj
-
Xgj
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'-'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Xgjcont
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'b'
)
#datas.append(data)
data
=
pt
.
DataPoints
(
TIMES
,
Ygjcont
,
tpe
=
'line'
,
label
=
'Ys'
,
color
=
'g'
)
#datas.append(data)
data
=
pt
.
DataPoints
(
TIMES
,
Ygjcont
-
Xgjcont
,
tpe
=
'line'
,
label
=
'Ys-Xs'
,
color
=
'r'
,
linestyle
=
'--'
)
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
,
ls
=
d
.
linestyle
)
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
opt
.
log
)
ymin
=
-
2
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
opt
.
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotEnergy
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Supernova\,Energy}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Es
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
if
opt
.
log
==
None
:
log
=
'y'
else
:
log
=
opt
.
log
+
'y'
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotEjectedMass
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Ejected\,Mass}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Mecont
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Mej
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Mejcont
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'g'
)
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
if
opt
.
log
==
None
:
log
=
opt
.
log
else
:
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotCumEjectedMass
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Cumulative\,Ejected\,Mass}$"
pt
.
subplot
(
np
,
1
,
i
)
cumcont
=
add
.
accumulate
(
Mecont
)
cumj
=
add
.
accumulate
(
Mej
)
cumjcont
=
add
.
accumulate
(
Mejcont
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
cumcont
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
cumj
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
cumjcont
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'g'
)
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
if
opt
.
log
==
None
:
log
=
opt
.
log
else
:
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotStellarMass
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Stellar\,Mass}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Mstaconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Mstajs
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Mstajconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'g'
)
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
if
opt
.
log
==
None
:
log
=
opt
.
log
else
:
log
=
opt
.
log
opt
.
ymin
=
0
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotNDYIN
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Number\,dying\,stars}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
NDYINs
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
NDYINconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotCumNDYIN
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Cum\,Number\,of\,dying\,stars}$"
pt
.
subplot
(
np
,
1
,
i
)
cum
=
add
.
accumulate
(
NDYINs
)
cumc
=
add
.
accumulate
(
NDYINconts
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
cum
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
cumc
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotNSNII
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Number\,of\,SNII}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
NSNIIs
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
NSNIIconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotCumNSNII
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Cum\,Number\,of\,SNII}$"
pt
.
subplot
(
np
,
1
,
i
)
cum
=
add
.
accumulate
(
NSNIIs
)
cumc
=
add
.
accumulate
(
NSNIIconts
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
cum
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
cumc
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotNSNIa
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Number\,of\,SNIa}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
NSNIas
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
NSNIaconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotCumNSNIa
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Cum\,Number\,of\,SNIa}$"
pt
.
subplot
(
np
,
1
,
i
)
cum
=
add
.
accumulate
(
NSNIas
)
cumc
=
add
.
accumulate
(
NSNIaconts
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
cum
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
cumc
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotNSNIIDYIN
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Number\,of\,SNII\,and\,dying}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
NSNIIs
+
NDYINs
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
NSNIIconts
+
NDYINconts
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotFevsTime
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Fe}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Ygj
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'k'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
TIMES
,
Ygjcont
,
tpe
=
'line'
,
label
=
'Es'
,
color
=
'r'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotDiff
(
i
):
xlabel
=
r"$\rm{[
%s
/H]$"
%
(
opt
.
x
)
ylabel
=
r"$\rm{[
%s
/
%s
]}$"
%
(
opt
.
y
,
opt
.
x
)
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
Xgcont
,
Ygcont
-
Xgcont
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'r'
)
datas
.
append
(
data
)
data
=
pt
.
DataPoints
(
Xgj
,
Ygj
-
Xgj
,
tpe
=
'line'
,
label
=
'Xs'
,
color
=
'k'
)
datas
.
append
(
data
)
#data = pt.DataPoints(Xgjcont,Ygjcont-Xgjcont,tpe='line',label='Xs',color='g')
#datas.append(data)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
#xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
xmin
=
-
4
xmax
=
0
ymin
=
-
1
ymax
=
2
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
def
PlotStarMass
(
i
):
xlabel
=
r'$\rm{Time} [\rm{
%s
}]$'
%
(
opt
.
timeunit
)
ylabel
=
r"$\rm{Mass}$"
pt
.
subplot
(
np
,
1
,
i
)
datas
=
[]
data
=
pt
.
DataPoints
(
TIMES
,
Mstas
,
tpe
=
'line'
,
label
=
'Mstas'
,
color
=
'k'
)
datas
.
append
(
data
)
for
d
in
datas
:
if
d
.
tpe
==
'line'
or
d
.
tpe
==
'both'
:
pt
.
plot
(
d
.
x
,
d
.
y
,
color
=
d
.
color
)
log
=
opt
.
log
xmin
,
xmax
,
ymin
,
ymax
=
pt
.
SetLimitsFromDataPoints
(
opt
.
xmin
,
opt
.
xmax
,
opt
.
ymin
,
opt
.
ymax
,
datas
,
log
)
# plot
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
log
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
pdict
=
{
'Y'
:
PlotYields
,
'Yg'
:
PlotGasElts
,
'ESN'
:
PlotEnergy
,
'NSNII'
:
PlotNSNII
,
'CumNSNII'
:
PlotCumNSNII
,
'NSNIa'
:
PlotNSNIa
,
'CumNSNIa'
:
PlotCumNSNIa
,
'NDYIN'
:
PlotNDYIN
,
'CumNDYIN'
:
PlotCumNDYIN
,
'D'
:
PlotDiff
,
'MSTAR'
:
PlotStarMass
,
'EjMass'
:
PlotEjectedMass
,
'CumEjMass'
:
PlotCumEjectedMass
,
'StellarMass'
:
PlotStellarMass
,
'NSNIIDYIN'
:
PlotNSNIIDYIN
,
'Fe'
:
PlotFevsTime
}
ip
=
1
for
obs
in
opt
.
obs
:
pdict
[
obs
](
ip
)
ip
=
ip
+
1
pt
.
EndPlot
(
files
,
opt
)
Event Timeline
Log In to Comment