-
Notifications
You must be signed in to change notification settings - Fork 9
/
AnisCoefficients_pix.py
380 lines (282 loc) · 11 KB
/
AnisCoefficients_pix.py
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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
from __future__ import division
import numpy as np
import healpy as hp
import math
import scipy.linalg as sl, scipy.special as ss
"""
Script to compute the correlation basis-functions for various anisotropic
configurations of the GW background energy-density
-- Rutger van Haasteren (June 2014)
-- Stephen Taylor (modifications, February 2016)
"""
def real_sph_harm(mm, ll, phi, theta):
"""
The real-valued spherical harmonics.
"""
if mm>0:
ans = (1./math.sqrt(2)) * \
(ss.sph_harm(mm, ll, phi, theta) + \
((-1)**mm) * ss.sph_harm(-mm, ll, phi, theta))
elif mm==0:
ans = ss.sph_harm(0, ll, phi, theta)
elif mm<0:
ans = (1./(math.sqrt(2)*complex(0.,1))) * \
(ss.sph_harm(-mm, ll, phi, theta) - \
((-1)**mm) * ss.sph_harm(mm, ll, phi, theta))
return ans.real
def signalResponse_fast(ptheta_a, pphi_a, gwtheta_a, gwphi_a):
"""
Create the signal response matrix FAST
"""
npsrs = len(ptheta_a)
# Create a meshgrid for both phi and theta directions
gwphi, pphi = np.meshgrid(gwphi_a, pphi_a)
gwtheta, ptheta = np.meshgrid(gwtheta_a, ptheta_a)
return createSignalResponse(pphi, ptheta, gwphi, gwtheta)
def createSignalResponse(pphi, ptheta, gwphi, gwtheta):
"""
Create the signal response matrix. All parameters are assumed to be of the
same dimensionality.
@param pphi: Phi of the pulsars
@param ptheta: Theta of the pulsars
@param gwphi: Phi of GW propagation direction
@param gwtheta: Theta of GW propagation direction
@return: Signal response matrix of Earth-term
"""
Fp = createSignalResponse_pol(pphi, ptheta, gwphi, gwtheta, plus=True)
Fc = createSignalResponse_pol(pphi, ptheta, gwphi, gwtheta, plus=False)
# Pixel maps are lumped together, polarization pixels are neighbours
F = np.zeros((Fp.shape[0], 2*Fp.shape[1]))
F[:, 0::2] = Fp
F[:, 1::2] = Fc
return F
def createSignalResponse_pol(pphi, ptheta, gwphi, gwtheta, plus=True, norm=True):
"""
Create the signal response matrix. All parameters are assumed to be of the
same dimensionality.
@param pphi: Phi of the pulsars
@param ptheta: Theta of the pulsars
@param gwphi: Phi of GW propagation direction
@param gwtheta: Theta of GW propagation direction
@param plus: Whether or not this is the plus-polarization
@param norm: Normalise the correlations to equal Jenet et. al (2005)
@return: Signal response matrix of Earth-term
"""
# Create the unit-direction vectors. First dimension will be collapsed later
# Sign convention of Gair et al. (2014)
Omega = np.array([-np.sin(gwtheta)*np.cos(gwphi), \
-np.sin(gwtheta)*np.sin(gwphi), \
-np.cos(gwtheta)])
mhat = np.array([-np.sin(gwphi), np.cos(gwphi), np.zeros(gwphi.shape)])
nhat = np.array([-np.cos(gwphi)*np.cos(gwtheta), \
-np.cos(gwtheta)*np.sin(gwphi), \
np.sin(gwtheta)])
p = np.array([np.cos(pphi)*np.sin(ptheta), \
np.sin(pphi)*np.sin(ptheta), \
np.cos(ptheta)])
# There is a factor of 3/2 difference between the Hellings & Downs
# integral, and the one presented in Jenet et al. (2005; also used by Gair
# et al. 2014). This factor 'normalises' the correlation matrix.
npixels = Omega.shape[2]
if norm:
# Add extra factor of 3/2
c = np.sqrt(1.5) / np.sqrt(npixels)
else:
c = 1.0 / np.sqrt(npixels)
# Calculate the Fplus or Fcross antenna pattern. Definitions as in Gair et
# al. (2014), with right-handed coordinate system
if plus:
# The sum over axis=0 represents an inner-product
Fsig = 0.5 * c * (np.sum(nhat * p, axis=0)**2 - np.sum(mhat * p, axis=0)**2) / \
(1 - np.sum(Omega * p, axis=0))
else:
# The sum over axis=0 represents an inner-product
Fsig = c * np.sum(mhat * p, axis=0) * np.sum(nhat * p, axis=0) / \
(1 - np.sum(Omega * p, axis=0))
return Fsig
def almFromClm(clm):
"""
Given an array of clm values, return an array of complex alm valuex
Note: There is a bug in healpy for the negative m values. This function just
takes the imaginary part of the abs(m) alm index.
"""
maxl = int(np.sqrt(len(clm)))-1
nclm = len(clm)
nalm = hp.Alm.getsize(maxl)
alm = np.zeros((nalm), dtype=np.complex128)
clmindex = 0
for ll in range(0, maxl+1):
for mm in range(-ll, ll+1):
almindex = hp.Alm.getidx(maxl, ll, abs(mm))
if mm == 0:
alm[almindex] += clm[clmindex]
elif mm < 0:
alm[almindex] -= 1j * clm[clmindex] / np.sqrt(2)
elif mm > 0:
alm[almindex] += clm[clmindex] / np.sqrt(2)
clmindex += 1
return alm
def clmFromAlm(alm):
"""
Given an array of clm values, return an array of complex alm valuex
Note: There is a bug in healpy for the negative m values. This function just
takes the imaginary part of the abs(m) alm index.
"""
nalm = len(alm)
maxl = int(np.sqrt(9.0 - 4.0 * (2.0-2.0*nalm))*0.5 - 1.5) # Really?
nclm = (maxl+1)**2
# Check the solution. Went wrong one time..
if nalm != int(0.5 * (maxl+1) * (maxl+2)):
raise ValueError("Check numerical precision. This should not happen")
clm = np.zeros(nclm)
clmindex = 0
for ll in range(0, maxl+1):
for mm in range(-ll, ll+1):
almindex = hp.Alm.getidx(maxl, ll, abs(mm))
if mm == 0:
clm[clmindex] = alm[almindex].real
elif mm < 0:
clm[clmindex] = - alm[almindex].imag * np.sqrt(2)
elif mm > 0:
clm[clmindex] = alm[almindex].real * np.sqrt(2)
clmindex += 1
return clm
def mapFromClm_fast(clm, nside):
"""
Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for
healpix pixelation with nside
@param clm: Array of C_{lm} values (inc. 0,0 element)
@param nside: Nside of the healpix pixelation
return: Healpix pixels
Use Healpix spherical harmonics for computational efficiency
"""
maxl = int(np.sqrt(len(clm)))-1
alm = almFromClm(clm)
h = hp.alm2map(alm, nside, maxl, verbose=False)
return h
def mapFromClm(clm, nside):
"""
Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for
healpix pixelation with nside
@param clm: Array of C_{lm} values (inc. 0,0 element)
@param nside: Nside of the healpix pixelation
return: Healpix pixels
Use real_sph_harm for the map
"""
npixels = hp.nside2npix(nside)
pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
h = np.zeros(npixels)
ind = 0
maxl = int(np.sqrt(len(clm)))-1
for ll in range(maxl+1):
for mm in range(-ll, ll+1):
h += clm[ind] * real_sph_harm(mm, ll, pixels[1], pixels[0])
ind += 1
return h
def clmFromMap_fast(h, lmax):
"""
Given a pixel map, and a maximum l-value, return the corresponding C_{lm}
values.
@param h: Sky power map
@param lmax: Up to which order we'll be expanding
return: clm values
Use Healpix spherical harmonics for computational efficiency
"""
alm = hp.sphtfunc.map2alm(h, lmax=lmax)
alm[0] = np.sum(h) * np.sqrt(4*np.pi) / len(h) # Why doesn't healpy do this?
return clmFromAlm(alm)
def clmFromMap(h, lmax):
"""
Given a pixel map, and a maximum l-value, return the corresponding C_{lm}
values.
@param h: Sky power map
@param lmax: Up to which order we'll be expanding
return: clm values
Use real_sph_harm for the map
"""
npixels = len(h)
nside = hp.npix2nside(npixels)
pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
clm = np.zeros( (lmax+1)**2 )
ind = 0
for ll in range(lmax+1):
for mm in range(-ll, ll+1):
clm[ind] += np.sum(h * real_sph_harm(mm, ll, pixels[1], pixels[0]))
ind += 1
return clm * 4 * np.pi / npixels
def getCov(clm, nside, F_e):
"""
Given a vector of clm values, construct the covariance matrix
@param clm: Array with Clm values
@param nside: Healpix nside resolution
@param F_e: Signal response matrix
@return: Cross-pulsar correlation for this array of clm values
"""
# Create a sky-map (power)
# Use mapFromClm to compare to real_sph_harm. Fast uses Healpix
#sh00 = mapFromClm(clm, nside)
sh00 = mapFromClm_fast(clm, nside)
# Double the power (one for each polarization)
sh = np.array([sh00, sh00]).T.flatten()
# Create the cross-pulsar covariance
hdcov_F = np.dot(F_e * sh, F_e.T)
# The pulsar term is added (only diagonals: uncorrelated)
return hdcov_F + np.diag(np.diag(hdcov_F))
def CorrBasis(psr_locs, lmax, nside=32):
"""
Calculate the correlation basis matrices using the pixel-space
transormations
@param psr_locs: Location of the pulsars [phi, theta]
@param lmax: Maximum l to go up to
@param nside: What nside to use in the pixelation [32]
Note: GW directions are in direction of GW propagation
"""
npsrs = len(psr_locs)
pphi = psr_locs[:,0]
ptheta = psr_locs[:,1]
# Create the pixels
npixels = hp.nside2npix(nside)
pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
gwtheta = pixels[0]
gwphi = pixels[1]
# Create the signal response matrix
F_e = signalResponse_fast(ptheta, pphi, gwtheta, gwphi)
# Loop over all (l,m)
basis = []
nclm = (lmax+1)**2
clmindex = 0
for ll in range(0, lmax+1):
for mm in range(-ll, ll+1):
clm = np.zeros(nclm)
clm[clmindex] = 1.0
basis.append(getCov(clm, nside, F_e))
clmindex += 1
return basis
def orfFromMap_fast(psr_locs, usermap, response=None):
"""
Calculate the correlation basis matrices using the pixel-space
transormations
@param psr_locs: Location of the pulsars [phi, theta]
@param usermap: Provide a healpix map for GW power
Note: GW directions are in direction of GW propagation
"""
if response is None:
npsrs = len(psr_locs)
pphi = psr_locs[:,0]
ptheta = psr_locs[:,1]
# Create the pixels
nside = hp.npix2nside(len(usermap))
npixels = hp.nside2npix(nside)
pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
gwtheta = pixels[0]
gwphi = pixels[1]
# Create the signal response matrix
F_e = signalResponse_fast(ptheta, pphi, gwtheta, gwphi)
elif response is not None:
F_e = response
# Double the power (one for each polarization)
sh = np.array([usermap, usermap]).T.flatten()
# Create the cross-pulsar covariance
hdcov_F = np.dot(F_e * sh, F_e.T)
# The pulsar term is added (only diagonals: uncorrelated)
return hdcov_F + np.diag(np.diag(hdcov_F))