1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
#!/usr/bin/env python3
#
# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT license
#
"""
Cosmology calculator with support of Chandra ACIS-specific quantities.
"""
import sys
import argparse
from collections import OrderedDict
from _context import acispy
from acispy.cosmo import Calculator
# Supported quantities
QUANTITIES = OrderedDict([
("luminosity_distance", {
"unit": "Mpc",
"label": "Luminosity distance",
"kwargs": ["z", "unit"],
}),
("angular_diameter_distance", {
"unit": "Mpc",
"label": "Angular diameter distance",
"kwargs": ["z", "unit"],
}),
("kpc_per_arcsec", {
"unit": None,
"label": "kpc/arcsec (DA)",
"kwargs": ["z"],
}),
("kpc_per_pix", {
"unit": None,
"label": "kpc/pix (DA)",
"kwargs": ["z"],
}),
("cm_per_pix", {
"unit": None,
"label": "cm/pix (DA)",
"kwargs": ["z"],
}),
("norm_apec", {
"unit": "cm^-5",
"label": "norm (APEC)",
"kwargs": ["z"],
}),
("evolution_factor", {
"unit": None,
"label": "Evolution factor",
"kwargs": ["z"],
}),
])
def get_quantities(args):
# convert ``argparse.Namespace`` to a dictionary
args = vars(args)
q_all = list(QUANTITIES.keys())
q_active = [q for q in q_all if args[q]]
if len(q_active) == 0:
q_active = q_all
return q_active
def calc_quantity(q, calculator, args):
args = vars(args)
kwargs = {arg: args[arg] for arg in QUANTITIES[q]["kwargs"]
if args[arg] is not None}
value = getattr(calculator, q)(**kwargs)
label = QUANTITIES[q]["label"]
unit = args["unit"] if args["unit"] is not None else QUANTITIES[q]["unit"]
if args["brief"]:
print(value)
else:
print("%s: %s # [%s]" % (label, value, unit))
def main():
# Hubble parameter at present day
H0 = 71.0 # [km/s/Mpc]
# Present-day matter density
Om0 = 0.27
# Chandra ACIS pixel size
pixelsize = 0.492 # [arcsec]
parser = argparse.ArgumentParser(
description="Cosmology calculator with Chandra-specific quantities")
parser.add_argument("-H", "--hubble", dest="H0",
type=float, default=H0,
help="present-day Hubble parameter " +
"(default: %s [km/s/Mpc])" % H0)
parser.add_argument("-M", "--omega-m", dest="Om0",
type=float, default=Om0,
help="present-day matter density (default: %s)" % Om0)
parser.add_argument("-b", "--brief", dest="brief", action="store_true",
help="be brief")
parser.add_argument("-u", "--unit", dest="unit",
help="set the unit for output quantity if supported")
parser.add_argument("-E", "--evolution-factor",
dest="evolution_factor",
action="store_true",
help="calculate the redshift evolution factor E(z)")
parser.add_argument("-L", "--luminosity-distance",
dest="luminosity_distance",
action="store_true",
help="calculate the luminosity distance (DL)")
parser.add_argument("-A", "--angular-diameter-distance",
dest="angular_diameter_distance",
action="store_true",
help="calculate the angular diameter distance (DA)")
parser.add_argument("--kpc-per-arcsec", dest="kpc_per_arcsec",
action="store_true",
help="calculate the transversal length [kpc] " +
"w.r.t. 1 arcsec at DA(z)")
parser.add_argument("--kpc-per-pix", dest="kpc_per_pix",
action="store_true",
help="calculate the transversal length [kpc] " +
"w.r.t. 1 pixel (%s [arcsec]) at DA(z)" % pixelsize)
parser.add_argument("--cm-per-pix", dest="cm_per_pix",
action="store_true",
help="calculate the transversal length [cm] " +
"w.r.t. 1 pixel (%s [arcsec]) at DA(z)" % pixelsize)
parser.add_argument("--norm-apec", dest="norm_apec",
action="store_true",
help="calculate the normalization factor " +
"of the XSPEC APEC model assuming EM=1")
parser.add_argument("z", type=float, help="redshift")
args = parser.parse_args()
cosmocalc = Calculator(H0=args.H0, Om0=args.Om0)
q_active = get_quantities(args)
if len(q_active) > 1:
# Multiple quantities to be calculated, thus unit specification
# and brief output are disabled.
if args.unit is not None:
args.unit = None
print("WARNING: ignored argument --unit", file=sys.stderr)
if args.brief:
args.brief = False
print("WARNING: ignored argument --brief", file=sys.stderr)
for q in q_active:
calc_quantity(q=q, calculator=cosmocalc, args=args)
if __name__ == "__main__":
main()
|
