-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathstrain.py
237 lines (174 loc) · 8.87 KB
/
strain.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
"""Computes several types of gravitational wave strains"""
import astropy.units as u
import astropy.constants as c
from legwork import utils
import numpy as np
__all__ = ['amplitude_modulation', 'h_0_n', 'h_c_n']
def amplitude_modulation(position, polarisation, inclination):
"""Computes the modulation of the strain due to the orbit averaged response of the detector to the
position, polarisation, and inclination of the source using Cornish+03 Eq.42 and Babak+21.
Note that since the majority of the calculations in LEGWORK are carried out for the full position,
polarisation, and inclination averages, we include a pre-factor of 5/4 on the amplitude modulation
to undo the factor of 4/5 which arises from the averaging of :meth:`legwork.utils.F_plus_squared` and
:meth:`legwork.utils.F_cross_squared`.
Additionally, note that this function does not include the factor of 1/2 from the Cornish+03 paper in
order to remain in the frequency domain. More recent papers (e.g. Babak+21) define the strain as
h(f)~A(f)*e^(2*i*psi(f)) and so the inner product is 1 instead of 1/2 as in Cornish+03.
Parameters
----------
position : `SkyCoord/array`, optional
Sky position of source. Must be specified using Astropy's :class:`astropy.coordinates.SkyCoord` class.
polarisation : `float/array`, optional
GW polarisation angle of the source. Must have astropy angular units.
inclination : `float/array`, optional
Inclination of the source. Must have astropy angular units.
Returns
-------
modulation : `float/array`
modulation to apply to strain from detector response
"""
# the pi/2 subtraction ensures that theta is a co-latitude to match Cornish+03
theta, phi = (np.pi/2 * u.rad) - position.lat, position.lon
# a_plus/a_cross as defined in Robson+19 Eq.15 and Babak+21 Eq. 67
a_plus = (1 + np.cos(inclination)**2) / 2
a_cross = np.cos(inclination)
term1 = a_plus**2 * utils.F_plus_squared(theta, phi, polarisation)
term2 = a_cross**2 * utils.F_cross_squared(theta, phi, polarisation)
# The modulation first undoes the averaging with the 5/4 factor
# This is because the full averaged response is 4/5 (e.g. see Robson+19 Eq.16)
modulation = 5/4 * (term1 + term2)
return modulation
def h_0_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclination=None, interpolated_g=None):
"""Computes strain amplitude
Computes the dimensionless power of a general binary radiating gravitational waves in the quadrupole
approximation at ``n`` harmonics of the orbital frequency
In the docs below, `x` refers to the number of sources, `y` to the number of timesteps and `z`
to the number of harmonics.
Parameters
----------
m_c : `float/array`
Chirp mass of each binary. Shape should be (x,).
f_orb : `float/array`
Orbital frequency of each binary at each timestep. Shape should be (x, y),
or (x,) if only one timestep.
ecc : `float/array`
Eccentricity of each binary at each timestep. Shape should be (x, y), or (x,) if only one timestep.
n : `int/array`
Harmonic(s) at which to calculate the strain. Either a single int or shape should be (z,)
dist : `float/array`
Distance to each binary. Shape should be (x,)
position : `SkyCoord/array`, optional
Sky position of source. Must be specified using Astropy's :class:`astropy.coordinates.SkyCoord` class.
polarisation : `float/array`, optional
GW polarisation angle of the source. Must have astropy angular units.
inclination : `float/array`, optional
Inclination of the source. Must have astropy angular units.
interpolated_g : `function`
A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
from Peters (1964). Default is None and uses exact g(n,e) in this case.
Returns
-------
h_0 : `float/array`
Strain amplitude. Shape is (x, y, z).
"""
# convert to array if necessary
arrayed_args, _ = utils.ensure_array(m_c, f_orb, ecc, n, dist)
m_c, f_orb, ecc, n, dist = arrayed_args
# raise a value error if any n is less than 1
if np.any(n < 1):
raise ValueError("All harmonics must be greater than or equal to 1")
# if one timestep then extend dimensions
if f_orb.ndim != 2:
f_orb = f_orb[:, np.newaxis]
if ecc.ndim != 2:
ecc = ecc[:, np.newaxis]
# extend mass and distance dimensions
m_c = m_c[:, np.newaxis]
dist = dist[:, np.newaxis]
# work out strain for n independent part and broadcast to correct shape
prefac = (2**(28/3) / 5)**(0.5) * c.G**(5/3) / c.c**4
n_independent_part = prefac * m_c**(5/3) * (np.pi * f_orb)**(2/3) / dist
# extend harmonic and eccentricity dimensions to full (x, y, z)
n = n[np.newaxis, np.newaxis, :]
ecc = ecc[..., np.newaxis]
# check whether to interpolate g(n, e)
if interpolated_g is None:
n_dependent_part = utils.peters_g(n, ecc)**(1/2) / n
else:
g_vals = interpolated_g(n, ecc)
# set negative values from cubic fit to 0.0
g_vals[g_vals < 0.0] = 0.0
n_dependent_part = g_vals**(0.5) / n
h_0 = n_independent_part[..., np.newaxis] * n_dependent_part
if position is not None or polarisation is not None or inclination is not None:
amp_mod = amplitude_modulation(position=position, polarisation=polarisation, inclination=inclination)
amp_mod = amp_mod[..., np.newaxis, np.newaxis]
h_0 = amp_mod**0.5 * h_0
return h_0.decompose().value
def h_c_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclination=None, interpolated_g=None):
"""Computes characteristic strain amplitude
Computes the dimensionless characteristic power of a general binary radiating gravitational waves in the
quadrupole approximation at ``n`` harmonics of the orbital frequency
In the docs below, `x` refers to the number of sources, `y` to the number of timesteps and `z`
to the number of harmonics.
Parameters
----------
m_c : `float/array`
Chirp mass of each binary. Shape should be (x,).
f_orb : `float/array`
Orbital frequency of each binary at each timestep. Shape should be (x, y),
or (x,) if only one timestep.
ecc : `float/array`
Eccentricity of each binary at each timestep. Shape should be (x, y), or (x,) if only one timestep.
n : `int/array`
Harmonic(s) at which to calculate the strain. Either a single int or shape should be (z,)
dist : `float/array`
Distance to each binary. Shape should be (x,)
position : `SkyCoord/array`, optional
Sky position of source. Must be specified using Astropy's :class:`astropy.coordinates.SkyCoord` class.
polarisation : `float/array`, optional
GW polarisation angle of the source. Must have astropy angular units.
inclination : `float/array`, optional
Inclination of the source. Must have astropy angular units.
interpolated_g : `function`
A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
from Peters (1964). Default is None and uses exact g(n,e) in this case.
Returns
-------
h_c : `float/array`
Characteristic strain. Shape is (x, y, z).
"""
# convert to array if necessary
arrayed_args, _ = utils.ensure_array(m_c, f_orb, ecc, n, dist)
m_c, f_orb, ecc, n, dist = arrayed_args
# raise a value error if any n is less than 1
if np.any(n < 1):
raise ValueError("All harmonics must be greater than or equal to 1")
# if one timestep then extend dimensions
if f_orb.ndim != 2:
f_orb = f_orb[:, np.newaxis]
if ecc.ndim != 2:
ecc = ecc[:, np.newaxis]
# extend mass and distance dimensions
m_c = m_c[:, np.newaxis]
dist = dist[:, np.newaxis]
# work out strain for n independent part
prefac = (2**(5/3) / (3 * np.pi**(4/3)))**(0.5) * c.G**(5/6) / c.c**(3/2)
n_independent_part = prefac * m_c**(5/6) / dist * f_orb**(-1/6) / utils.peters_f(ecc)**(0.5)
# extend harmonic and eccentricity dimensions to full (x, y, z)
n = n[np.newaxis, np.newaxis, :]
ecc = ecc[..., np.newaxis]
# check whether to interpolate g(n, e)
if interpolated_g is None:
n_dependent_part = (utils.peters_g(n, ecc) / n)**(1/2)
else:
g_vals = interpolated_g(n, ecc)
# set negative values from cubic fit to 0.0
g_vals[g_vals < 0.0] = 0.0
n_dependent_part = (g_vals / n)**(0.5)
h_c = n_independent_part[..., np.newaxis] * n_dependent_part
if position is not None or polarisation is not None or inclination is not None:
amp_mod = amplitude_modulation(position=position, polarisation=polarisation, inclination=inclination)
amp_mod = amp_mod[..., np.newaxis, np.newaxis]
h_c = amp_mod**0.5 * h_c
return h_c.decompose().value