20 """This module provides a Python class for generating antenna response
21 functions. By default the class uses its own implementation of the antenna
22 response functions, but it also provides a wrapper to the equivalent
124 from six
import string_types
128 from lal
import (LALDetectorIndexLHODIFF, LALDetectorIndexLLODIFF,
129 LALDetectorIndexGEO600DIFF, LALDetectorIndexVIRGODIFF,
130 LALDetectorIndexTAMA300DIFF, LALDetectorIndexKAGRADIFF,
131 LALDetectorIndexLIODIFF, LALDetectorIndexE1DIFF,
132 LALDetectorIndexE2DIFF, LALDetectorIndexE3DIFF,
133 CachedDetectors, LIGOTimeGPS, GreenwichMeanSiderealTime,
134 ComputeDetAMResponse, ComputeDetAMResponseExtraModes,
135 DAYSID_SI, Detector, TranslateHMStoRAD, TranslateDMStoRAD)
137 from .
import git_version
139 __author__ =
"Matthew Pitkin <matthew.pitkin@ligo.org>"
140 __version__ = git_version.verbose_msg
141 __date__ = git_version.date
145 DETMAP = {
'H1': LALDetectorIndexLHODIFF,
146 'H2': LALDetectorIndexLHODIFF,
147 'LHO': LALDetectorIndexLHODIFF,
148 'L1': LALDetectorIndexLLODIFF,
149 'LLO': LALDetectorIndexLLODIFF,
150 'G1': LALDetectorIndexGEO600DIFF,
151 'GEO': LALDetectorIndexGEO600DIFF,
152 'GEO600': LALDetectorIndexGEO600DIFF,
153 'V1': LALDetectorIndexVIRGODIFF,
154 'VIRGO': LALDetectorIndexVIRGODIFF,
155 'T1': LALDetectorIndexTAMA300DIFF,
156 'TAMA': LALDetectorIndexTAMA300DIFF,
157 'TAMA300': LALDetectorIndexTAMA300DIFF,
158 'K1': LALDetectorIndexKAGRADIFF,
159 'KAGRA': LALDetectorIndexKAGRADIFF,
160 'LCGT': LALDetectorIndexKAGRADIFF,
161 'I1': LALDetectorIndexLIODIFF,
162 'LIO': LALDetectorIndexLIODIFF,
163 'E1': LALDetectorIndexE1DIFF,
164 'E2': LALDetectorIndexE2DIFF,
165 'E3': LALDetectorIndexE3DIFF}
169 VAR_RANGES = {
'psi': [0., 2.*np.pi],
170 'time': [0., DAYSID_SI],
171 'ra': [0., 2.*np.pi],
176 def __init__(self, detector, ra, dec, psi=0., times=None, tensor=True,
177 vector=False, scalar=False, use_lal=False, lookup=False,
178 lookuppair=None, bins1=100, bins2=100):
180 Calculate the long-wavelength limit antenna response functions for a given
181 ground-based gravitational wave detector. The response can include tensor,
182 vector, and scalar modes.
184 @param detector: (str) a valid detector name (e.g., 'H1') or
186 @param ra: (array_like) the right ascension of the source in radians,
187 or a string of the format 'hh:mm:ss.s'.
188 @param dec: (array_like) the declination of the source in radians, or a
189 string of the format 'dd:mm:ss.s'.
190 @param psi: (array_like) the polarization angle in radians. Defaults to
192 @param times: (array_like) an array of GPS time values at which to
193 calculate the response function.
194 @param tensor: (bool) set to calculate and store the tensor
195 polarization components (plus and cross). Defaults to True.
196 @param vector: (bool) set to calculate and store the vector
197 polarization components ("x" and "y"). Defaults to False.
198 @param scalar: (bool) set to calculate and store the scalar
199 polarization components (longitudinal and breathing). Defaults to
201 @param use_lal: (bool) set to internally use the
202 XLALComputeDetAMResponse() and XLALComputeDetAMResponseExtraModes()
203 functions. Defaults to False.
204 @param lookup: (bool) set to generate and use a look-up table in a pair
205 of parameters for computing the antenna responses. Defaults to
206 False. If using the look-up table, the arrays of values being
207 "looked-up" must be in ascending order.
208 @param lookuppair: (list, tuple) a list of the two parameters that will
209 be used in the look-up table interpolation. Defaults to
210 <tt>['psi', 'time']</tt> (allowed values are <tt>'ra'</tt>,
211 <tt>'dec'</tt>, <tt>'psi'</tt>, or <tt>'time'</tt>)
212 @param bins1: (int) the number of bins in the grid in the first look-up
213 table parameter. Defaults to 100.
214 @param bins2: (int) the number of bins in the grid in the second
215 look-up table parameter. Defaults to 100.
217 Example usage for tensor polarizations:
219 >>> from lal.antenna import AntennaResponse
220 >>> # compute tensor response for a single time
221 >>> resp = AntennaResponse('H1', ra=1.2, dec=-0.3, psi=2.9,
223 >>> print('Fplus: {}'.format(resp.plus))
225 >>> print('Fcross: {}'.format(resp.cross))
226 Fcross: [-0.79809163]
227 >>> # re-use class to get response at multiple new times
228 >>> resp.compute_response([1010101010., 1234567890.])
229 >>> print('Fplus: {}'.format(resp.plus))
230 Fplus: [ 0.09498567 -0.45495654]
231 >>> print('Fcross: {}'.format(resp.cross))
232 Fcross: [0.1706959 0.21690418]
235 Example usage for tensor, vector and scalar polarizations (at a series
238 >>> import numpy as np
239 >>> times = np.linspace(1000000000.0, 1000086340.0, 1440)
240 >>> resp = AntennaResponse('H1', ra=1.2, dec=-0.3, psi=2.9,
241 ...scalar=True, vector=True, times=times)
243 array([0.32427018, 0.32805983, 0.3318344 , ..., 0.32780195, 0.33157755,
246 array([-0.79809163, -0.79607858, -0.79404097, ..., -0.7962166 ,
247 -0.79418066, -0.79212028])
248 >>> resp.x # vector "x" polarization
249 array([-0.46915186, -0.46773594, -0.46627224, ..., -0.46783399,
250 -0.46637354, -0.46486538])
251 >>> resp.y # vector "y" polarization
252 array([-0.17075718, -0.17475991, -0.17875012, ..., -0.17448742,
253 -0.17847849, -0.18245689])
254 >>> resp.b # scalar "breathing" mode
255 array([0.05365678, 0.05573073, 0.05780282, ..., 0.05558939, 0.05766162,
257 >>> resp.l # scalar "longitudinal mode"
258 array([-0.05365678, -0.05573073, -0.05780282, ..., -0.05558939,
259 -0.05766162, -0.05973181])
291 if lookuppair
is None:
294 self.
set_lookup_pairset_lookup_pair(pair=lookuppair, bins1=bins1, bins2=bins2)
307 Set the detector for which to calculate the antenna response.
309 @param det: (str) a valid detector name.
312 if isinstance(det, Detector):
313 self.
_detector_detector = det.frDetector.prefix
315 elif isinstance(det, string_types):
316 if det.upper()
not in DETMAP.keys():
317 raise ValueError(
"Detector is not a valid detector name")
324 raise TypeError(
"Detector must be a string or lal.Detector object")
333 Set the lal.Detector.
335 @param det: (str) a valid detector name.
338 if isinstance(det, Detector):
340 elif isinstance(det, string_types):
342 detector = DETMAP[det.upper()]
344 raise KeyError(
"Key {} is not a valid detector name.".format(det))
346 self.
_laldetector_laldetector = CachedDetectors[detector].response
348 raise TypeError(
"Detector must be a string or lal.Detector object")
357 Set the right ascension.
363 elif isinstance(raval, float)
or isinstance(raval, int):
364 self.
_ra_ra = np.array([raval], dtype=
'float64')
365 elif isinstance(raval, list)
or isinstance(raval, np.ndarray):
366 self.
_ra_ra = np.copy(raval).astype(
'float64')
367 elif isinstance(raval, string_types):
369 rarad = TranslateHMStoRAD(raval)
370 self.
_ra_ra = np.array([rarad], dtype=
'float64')
372 raise ValueError(
"Could not convert '{}' to a right "
373 "ascension".format(raval))
375 raise TypeError(
"Right ascension must be an array")
390 def dec(self, decval):
400 elif isinstance(decval, float)
or isinstance(decval, int):
401 self.
_dec_dec = np.array([decval], dtype=
'float64')
402 elif isinstance(decval, list)
or isinstance(decval, np.ndarray):
403 self.
_dec_dec = np.copy(decval).astype(
'float64')
404 elif isinstance(decval, string_types):
406 decrad = TranslateDMStoRAD(decval)
407 self.
_dec_dec = np.array([decrad], dtype=
'float64')
409 raise ValueError(
"Could not convert '{}' to a "
410 "declination".format(decval))
412 raise TypeError(
"Declination must be an array")
430 def psi(self, psival):
432 Set the value of the gravitational wave polarization angle psi.
434 @param psival: (float) the polarization angle (radians)
442 elif isinstance(psival, float)
or isinstance(psival, int):
443 self.
_psi_psi = np.array([psival], dtype=
'float64')
444 elif isinstance(psival, list)
or isinstance(psival, np.ndarray):
445 self.
_psi_psi = np.copy(psival).astype(
'float64')
447 raise TypeError(
"Polarization must be an array")
449 self.
_psi_psi = np.mod(self.
_psi_psi, 2.*np.pi)
461 except AttributeError:
465 def times(self, timearr):
467 Set array of times and GPS times.
475 elif isinstance(timearr, float)
or isinstance(timearr, int):
476 self.
_times_times = np.array([timearr], dtype=
'float64')
477 elif isinstance(timearr, list)
or isinstance(timearr, np.ndarray):
478 self.
_times_times = np.copy(timearr).astype(
'float64')
480 raise TypeError(
"Times must be an array")
484 for i, time
in enumerate(self.
_times_times):
485 gps = LIGOTimeGPS(time)
486 gmstrad = GreenwichMeanSiderealTime(gps)
514 def tensor(self, tensorval):
516 Set whether to include tensor polarizations.
519 if not isinstance(tensorval, bool):
520 raise TypeError(
"Must be boolean value")
522 self.
_tensor_tensor = tensorval
529 def vector(self, vectorval):
531 Set whether to include vector polarizations.
534 if not isinstance(vectorval, bool):
535 raise TypeError(
"Must be boolean value")
537 self.
_vector_vector = vectorval
544 def scalar(self, scalarval):
546 Set whether to include scalar polarizations.
549 if not isinstance(scalarval, bool):
550 raise TypeError(
"Must be boolean value")
552 self.
_scalar_scalar = scalarval
561 Set whether to use LAL antenna response functions.
564 if not isinstance(val, bool):
565 raise TypeError(
"Must be a boolean value")
575 if not isinstance(pair, (list, tuple)):
576 raise TypeError(
'Pair must be a list or tuple')
579 raise ValueError(
'Pair must only contain two values')
582 if val.lower()
not in self.
parametersparameters:
583 raise ValueError(
'Parameter {} in pair not '
584 'recognized'.format(val))
586 self.
_lookup_pair_lookup_pair = [val.lower()
for val
in pair]
600 Set the pair of parameters to use for the look-up table.
606 for val, nbins
in zip(pair, [bins1, bins2]):
607 if not isinstance(nbins, int):
608 raise TypeError(
"Value must be an integer")
611 raise ValueError(
"There must be at least 2 bins")
615 vrange = VAR_RANGES[val]
621 vararray = np.arcsin(np.linspace(vrange[0], vrange[1], nbins))
623 vararray = np.linspace(vrange[0], vrange[1], nbins)
651 Set the 2d look-up table.
656 if not isinstance(val, bool):
657 raise TypeError(
"Value must be a boolean")
664 from scipy.interpolate
import RectBivariateSpline
666 raise ImportError(
"Cannot import scipy")
673 except AttributeError:
674 raise AttributeError(
"A time must be set for the look-up table")
765 return self.
responseresponse[
'plus']
774 self.
responseresponse[
'plus'] = resp
778 return self.
responseresponse[
'cross']
787 self.
responseresponse[
'cross'] = resp
843 Convert one-dimensional arrays into a mesh. The meshes dimensions are
844 ordered as right ascension, declination, psi, and time.
852 ramesh, decmesh, psimesh, timemesh = np.meshgrid(self.
rararara, self.
decdecdecdec,
868 if isinstance(val, np.ndarray):
870 self.
_ra_mesh_ra_mesh = val.squeeze()
888 if isinstance(val, np.ndarray):
908 if isinstance(val, np.ndarray):
923 if isinstance(val, np.ndarray):
934 Compute the detector response.
936 @param times: (array_like) an array of GPS times at which to compute
937 the response function. If not set the times set at initialization
938 of the class, or using the <tt>times</tt> property.
941 if times
is not None:
955 def _compute_response(self):
957 Compute antenna pattern.
965 for i
in range(len(self.
shapeshape)):
966 einsum_indices += eindices[i]
969 einstr1 =
'i{},j{}->ij{}'.format(*3*[einsum_indices])
970 einstr2 =
'ij{},ij->{}'.format(*2*[einsum_indices])
979 mm = np.einsum(einstr1, M, M)
980 mn = np.einsum(einstr1, M, N)
981 nm = np.einsum(einstr1, N, M)
982 nn = np.einsum(einstr1, N, N)
996 mq = np.einsum(einstr1, M, Q)
997 qm = np.einsum(einstr1, Q, M)
998 nq = np.einsum(einstr1, N, Q)
999 qn = np.einsum(einstr1, Q, N)
1005 qq = np.einsum(einstr1, Q, Q)
1010 def _compute_response_lal(self):
1012 Compute antenna pattern using LAL functions.
1017 slen = np.prod(self.
shapeshape)
if len(self.
shapeshape) != 0
else 1
1021 fp = np.zeros(self.
shapeshape)
1022 fc = np.zeros(self.
shapeshape)
1026 antenna_func = ComputeDetAMResponse
1035 for i
in range(slen):
1036 idxs = np.unravel_index(i, self.
shapeshape)
1047 antenna_func = ComputeDetAMResponseExtraModes
1056 fb = np.zeros(self.
shapeshape)
1057 fl = np.zeros(self.
shapeshape)
1058 fx = np.zeros(self.
shapeshape)
1059 fy = np.zeros(self.
shapeshape)
1061 for i
in range(slen):
1062 idxs = np.unravel_index(i, self.
shapeshape)
1089 def _compute_response_lookup(self):
1094 fp = np.zeros(lushape)
1095 fc = np.zeros(lushape)
1098 fb = np.zeros(lushape)
1099 fl = np.zeros(lushape)
1102 fx = np.zeros(lushape)
1103 fy = np.zeros(lushape)
1107 unsorted_idxs =
None
1110 pairs.append(self.
rararara)
1118 unsorted_idxs = np.argsort(np.argsort(mtimes))
1120 pairs.append(mtimes)
1126 sshape = tuple(sshape)
1129 pos = 4*[slice(
None)]
1138 fp[tuple(pos)] = self.
_lookup_func_lookup_func[
'plus'][i, j](*pairs).reshape(sshape)
1139 fc[tuple(pos)] = self.
_lookup_func_lookup_func[
'cross'][i, j](*pairs).reshape(sshape)
1142 fb[tuple(pos)] = self.
_lookup_func_lookup_func[
'b'][i, j](*pairs).reshape(sshape)
1143 fl[tuple(pos)] = self.
_lookup_func_lookup_func[
'l'][i, j](*pairs).reshape(sshape)
1146 fx[tuple(pos)] = self.
_lookup_func_lookup_func[
'x'][i, j](*pairs).reshape(sshape)
1147 fy[tuple(pos)] = self.
_lookup_func_lookup_func[
'y'][i, j](*pairs).reshape(sshape)
1149 if unsorted_idxs
is not None:
1151 self.
plusplusplusplus = fp[:, :, :, unsorted_idxs].squeeze()
1155 self.
bbbb = fb[:, :, :, unsorted_idxs].squeeze()
1156 self.
llll = fl[:, :, :, unsorted_idxs].squeeze()
1159 self.
xxxx = fx[:, :, :, unsorted_idxs].squeeze()
1160 self.
yyyy = fy[:, :, :, unsorted_idxs].squeeze()
1167 self.
bbbb = fb.squeeze()
1168 self.
llll = fl.squeeze()
1171 self.
xxxx = fx.squeeze()
1172 self.
yyyy = fy.squeeze()
1177 def __call__(self, times, ra=None, dec=None, psi=None, detector=None):
1179 Return the antenna response function as a dictionary.
1182 if detector
is not None:
def vector(self, vectorval)
Set whether to include vector polarizations.
def __init__(self, detector, ra, dec, psi=0., times=None, tensor=True, vector=False, scalar=False, use_lal=False, lookup=False, lookuppair=None, bins1=100, bins2=100)
Calculate the long-wavelength limit antenna response functions for a given ground-based gravitational...
def psi(self, psival)
Set the value of the gravitational wave polarization angle psi.
def ra(self, raval)
Set the right ascension.
def compute_response(self, times=None)
Compute the detector response.
def tensor(self, tensorval)
Set whether to include tensor polarizations.
def detector(self, det)
Set the detector for which to calculate the antenna response.
def lookup(self, val)
Set the 2d look-up table.
def _compute_response_lookup(self)
def __call__(self, times, ra=None, dec=None, psi=None, detector=None)
Return the antenna response function as a dictionary.
def _compute_response(self)
def scalar(self, scalarval)
Set whether to include scalar polarizations.
def times(self, timearr)
Set array of times and GPS times.
def dec(self, decval)
Set the declination.
def use_lal(self, val)
Set whether to use LAL antenna response functions.
def laldetector(self, det)
Set the lal.Detector.
def set_lookup_pair(self, pair=['psi', 'time'], bins1=100, bins2=100)
Set the pair of parameters to use for the look-up table.
def lookup_pair(self, pair)
def _compute_response_lal(self)