|
1
|
|
|
"Some utility methods" |
|
2
|
|
|
from __future__ import division |
|
3
|
|
|
|
|
4
|
|
|
import numpy as np |
|
5
|
|
|
from numpy import cos |
|
6
|
|
|
from numpy import pi |
|
7
|
|
|
from numpy import sin |
|
8
|
|
|
|
|
9
|
|
|
|
|
10
|
|
|
def parse_units(val, units, name, dtype=float): |
|
11
|
|
|
""" |
|
12
|
|
|
Convert, as necessary, `val` to `units` as `dtype` |
|
13
|
|
|
""" |
|
14
|
|
|
|
|
15
|
|
|
try: |
|
16
|
|
|
return val.to(units).value |
|
17
|
|
|
except AttributeError: |
|
18
|
|
|
pass |
|
19
|
|
|
except Exception as inst: |
|
20
|
|
|
raise IOError("Bad format for input {}. ({})". |
|
21
|
|
|
format(name, inst)) |
|
22
|
|
|
try: |
|
23
|
|
|
parsed_val = np.array(val).astype(dtype) |
|
24
|
|
|
if parsed_val.size == 1: |
|
25
|
|
|
return parsed_val.flatten()[0] |
|
26
|
|
|
return parsed_val |
|
27
|
|
|
except Exception as inst: |
|
28
|
|
|
raise IOError("Bad format for input {}. ({})". |
|
29
|
|
|
format(name, inst)) |
|
30
|
|
|
|
|
31
|
|
|
|
|
32
|
|
|
def parse_DM(DM): |
|
33
|
|
|
""" |
|
34
|
|
|
Convert, as necessary, DM into float |
|
35
|
|
|
""" |
|
36
|
|
|
return parse_units(DM, 'pc/cm**3', 'DM') |
|
37
|
|
|
|
|
38
|
|
|
|
|
39
|
|
|
def parse_lbd(gal_l, gal_b, distance): |
|
40
|
|
|
""" |
|
41
|
|
|
Convert, as necessary, l,b,d into floats |
|
42
|
|
|
|
|
43
|
|
|
Parameters |
|
44
|
|
|
---------- |
|
45
|
|
|
in_l : float or Angle |
|
46
|
|
|
Galactic longitude; assumed deg if unitless |
|
47
|
|
|
in_b : float |
|
48
|
|
|
Galactic latitude; assumed deg if unitless |
|
49
|
|
|
in_d : float or Quantity |
|
50
|
|
|
Distance to source; assumed kpc if unitless |
|
51
|
|
|
|
|
52
|
|
|
Returns |
|
53
|
|
|
------- |
|
54
|
|
|
gal_l : float |
|
55
|
|
|
Galactic longitude in deg |
|
56
|
|
|
gal_b : float |
|
57
|
|
|
Galactic latitude in deg |
|
58
|
|
|
distance : float |
|
59
|
|
|
Distance in kpc |
|
60
|
|
|
|
|
61
|
|
|
""" |
|
62
|
|
|
l = parse_units(gal_l, 'deg', 'Galactic longitude') |
|
63
|
|
|
b = parse_units(gal_b, 'deg', 'Galactic latitude') |
|
64
|
|
|
d = parse_units(distance, 'kpc', 'distance') |
|
65
|
|
|
return l, b, d |
|
66
|
|
|
|
|
67
|
|
|
|
|
68
|
|
|
def galactic_to_galactocentric(l, b, distance, xyz_sun): |
|
69
|
|
|
""" Convert galactic coordiantes to galactocentric |
|
70
|
|
|
Parameters |
|
71
|
|
|
---------- |
|
72
|
|
|
l : float |
|
73
|
|
|
latitude |
|
74
|
|
|
b : float |
|
75
|
|
|
longitude |
|
76
|
|
|
distance : float |
|
77
|
|
|
kpc |
|
78
|
|
|
xyz_sun : ndarray |
|
79
|
|
|
positions of the Sun in kpc |
|
80
|
|
|
|
|
81
|
|
|
Returns |
|
82
|
|
|
------- |
|
83
|
|
|
xyz_c : ndarray |
|
84
|
|
|
x,y,z positions along the sightline |
|
85
|
|
|
|
|
86
|
|
|
""" |
|
87
|
|
|
slc = sin(l/180*pi) |
|
88
|
|
|
clc = cos(l/180*pi) |
|
89
|
|
|
sbc = sin(b/180*pi) |
|
90
|
|
|
cbc = cos(b/180*pi) |
|
91
|
|
|
rgalc = distance*cbc |
|
92
|
|
|
xc = xyz_sun[0] + rgalc*slc |
|
93
|
|
|
yc = xyz_sun[1] - rgalc*clc |
|
94
|
|
|
zc = xyz_sun[-1] + distance*sbc |
|
95
|
|
|
return np.array([xc, yc, zc]) |
|
96
|
|
|
|
|
97
|
|
|
|
|
98
|
|
|
def lzproperty(attribute): |
|
99
|
|
|
""" |
|
100
|
|
|
Lazy property: evaluate property only once |
|
101
|
|
|
""" |
|
102
|
|
|
save_att = '_' + attribute.__name__ |
|
103
|
|
|
|
|
104
|
|
|
@property |
|
105
|
|
|
def _get(self): |
|
106
|
|
|
try: |
|
107
|
|
|
return getattr(self, save_att) |
|
108
|
|
|
except AttributeError: |
|
109
|
|
|
setattr(self, save_att, attribute(self)) |
|
110
|
|
|
return getattr(self, save_att) |
|
111
|
|
|
return _get |
|
112
|
|
|
|
|
113
|
|
|
|
|
114
|
|
|
def rotation(theta, axis=-1): |
|
115
|
|
|
""" |
|
116
|
|
|
Return a rotation matrix around axis |
|
117
|
|
|
0:x, 1:y, 2:z |
|
118
|
|
|
""" |
|
119
|
|
|
ct = cos(theta) |
|
120
|
|
|
st = sin(theta) |
|
121
|
|
|
|
|
122
|
|
|
if axis in (0, -3): |
|
123
|
|
|
return np.array([[1, 0, 0], |
|
124
|
|
|
[0, ct, st], |
|
125
|
|
|
[0, -st, ct]]) |
|
126
|
|
|
|
|
127
|
|
|
if axis in (1, -2): |
|
128
|
|
|
return np.array([[ct, 0, st], |
|
129
|
|
|
[0, 1, 0], |
|
130
|
|
|
[-st, 0, ct]]) |
|
131
|
|
|
|
|
132
|
|
|
if axis in (2, -1): |
|
133
|
|
|
return np.array([[ct, st, 0], |
|
134
|
|
|
[-st, ct, 0], |
|
135
|
|
|
[0, 0, 1]]) |
|
136
|
|
|
|
|
137
|
|
|
|
|
138
|
|
|
def rad3d2(xyz): |
|
139
|
|
|
return xyz[0]**2 + xyz[1]**2 + xyz[-1]**2 |
|
140
|
|
|
|
|
141
|
|
|
|
|
142
|
|
|
def rad2d2(xyz): |
|
143
|
|
|
return xyz[0]**2 + xyz[1]**2 |
|
144
|
|
|
|
|
145
|
|
|
|
|
146
|
|
|
def matmul(a, b): |
|
147
|
|
|
try: |
|
148
|
|
|
return a.__matmul__(b) |
|
149
|
|
|
except AttributeError: |
|
150
|
|
|
return np.matmul(a, b) |
|
151
|
|
|
|
benbaror /
ne2001
