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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
|
#!/usr/bin/env python3
#
# Deproject the 2D surface brightness profile (SBP) into the 3D emission
# measure (EM) profile, using the non-parametric approach.
# And the 3D gas density profile can be further derived through the 3D
# EM profile by taking into account the variation of the cooling function
# Lambda(T, Z) with radius.
#
#
# References:
# [1] Croston et al. 2006, A&A, 459, 1007-1019
# [2] McLaughlin, 1999, ApJ, 117, 2398-2427
#
#
# Weitian LI
# Created: 2016-06-10
# Updated: 2016-06-10
#
import argparse
import numpy as np
class Projection:
"""
Class that deals with projection from 3D volume density to 2D
surface density and vice versa.
The inner-most shell/cylinder is assumed to at the center with inner
radius of ZERO.
"""
# number of shells/cylinders
N = 0
# inner and outer radii of each spherical shell or cylinders
rin = None
rout = None
# projection matrix from 3D volume density to 2D surface density
proj_mat = None
def __init__(self, rout):
self.N = len(rout)
self.rout = np.array(rout, dtype=float)
self.rin = np.concatenate([[0.0], self.rout[:-1]])
self.calc_projection_matrix()
def __str__(self):
return "%s: #%d shells: Rout(%s)" % (self.__class__.__name__,
self.N, self.rout)
def calc_projection_matrix(self):
"""
Calculate the projection matrix according to the given outer radii.
Arguments:
* rout: (vector) outer radius of each SB annulus or spherical shell
Return:
* proj_mat: (matrix) an upper triangular matrix with element
[i, j] indicate the fraction of the emission from
shell j that is observed in annulus i.
N(R_{i-1}, R_i) * \pi * (R^2_i - R^2_{i-1}) =
\sum_{j=i}^{m} (n(R_{j-1}, R_j) *
V_int(R_{j-1}, R_j; R_{i-1}, R_i))
References:
* ref.[1], eq.(1)
* ref.[2], eq.(A2)
"""
proj_mat = np.zeros((self.N, self.N))
for i in range(self.N):
# loop over each annulus
rin = self.rin[i]
rout = self.rout[i]
area = np.pi * (rout**2 - rin**2)
for j in range(i, self.N):
# calculate the contribution from each shell to annulus i
rin2 = self.rin[j]
rout2 = self.rout[j]
v_int = self.intersection_volume(rin2, rout2, rin, rout)
proj_mat[i, j] = v_int / area
self.proj_mat = proj_mat
def project(self, densities):
"""
Project the given 3D (volume) densities to 2D (surface) densities,
using the calculated projection matrix: 'proj_mat'.
"""
densities = np.array(densities)
if self.rout.shape != densities.shape:
raise ValueError("different shapes of rout and given densities")
return self.proj_mat.dot(densities.T)
def deproject(self, densities):
"""
Revert the projection procedure, i.e., deproject the given 2D
(surface) densities to derive the 3D (volume) densities.
\curl{N}(R_{i-1}, R_i) = N(R_{i-1}, R_i) * \pi * (R^2_i - R^2_{i-1})
n(R_{i-1}, R_i) =
(N(R_{i-1}, R_i) * \pi * (R^2_i - R^2_{i-1}) /
V_int(R_{i-1}, R_i; R_{i-1}, R_i)) -
\sum_{j=i+1}^{m} (n(R_{j-1}, R_j) *
V_int(R_{j-1}, R_j; R_{i-1}, R_i) /
V_int(R_{i-1}, R_i; R_{i-1}, R_i))
Reference: ref.[2], eq.(A2)
"""
densities = np.array(densities)
if self.rout.shape != densities.shape:
raise ValueError("different shapes of rout and given densities")
n_3d = np.zeros(densities.shape)
# peel the onion: from outside inward
for i in reversed(range(self.N)):
rin = self.rin[i]
rout = self.rout[i]
area = np.pi * (rout**2 - rin**2)
v_int = self.intersection_volume(rin, rout, rin, rout)
n_3d[i] = densities[i] * area / v_int
# subtract the projections from the outer shells
for j in range(i+1, self.N):
rin2 = self.rin[j]
rout2 = self.rout[j]
v_int2 = self.intersection_volume(rin2, rout2, rin, rout)
n_3d[i] -= n_3d[j] * v_int2 / v_int
return n_3d
@staticmethod
def intersection_volume(r1, r2, R1, R2):
"""
Calculate the volume of intersection between the spherical shell of
r1 <= r <= r2 and the cylinder of R1 <= R <= R2.
Reference: ref.[2], eq.(A1)
"""
def trunc_pow(x, p):
if x <= 0.0:
return 0
else:
return x ** p
#
v_int = (4.0*np.pi/3.0) * (trunc_pow((r2**2 - R1**2), 1.5) -
trunc_pow((r2**2 - R2**2), 1.5) +
trunc_pow((r1**2 - R2**2), 1.5) -
trunc_pow((r1**2 - R1**2), 1.5))
return v_int
def testProjection():
rout = np.array([1, 2, 3, 4, 5], dtype=float)
proj = Projection(rout)
n1 = np.array([1, 1, 1, 1, 1], dtype=float)
np.testing.assert_array_almost_equal(proj.deproject(proj.project(n1)), n1)
s2 = np.array([1, 1, 1, 1, 1], dtype=float)
np.testing.assert_array_almost_equal(proj.project(proj.deproject(s2)), s2)
print("All tests PASSED!")
def main():
parser = argparse.ArgumentParser(
description="Deproject the surface brightness profile (SBP)")
args = parser.parse_args()
if __name__ == "__main__":
main()
|
