Coverage for bilby/gw/source.py: 81%
274 statements
« prev ^ index » next coverage.py v7.6.1, created at 2025-05-06 04:57 +0000
1import numpy as np
3from ..core import utils
4from ..core.utils import logger
5from .conversion import bilby_to_lalsimulation_spins
6from .utils import (lalsim_GetApproximantFromString,
7 lalsim_SimInspiralFD,
8 lalsim_SimInspiralChooseFDWaveform,
9 lalsim_SimInspiralChooseFDWaveformSequence)
11UNUSED_KWARGS_MESSAGE = """There are unused waveform kwargs. This is deprecated behavior and will
12result in an error in future releases. Make sure all of the waveform kwargs are correctly
13spelled.
15Unused waveform_kwargs: {waveform_kwargs}
16"""
19def gwsignal_binary_black_hole(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
20 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
21 """
22 A binary black hole waveform model using GWsignal
24 Parameters
25 ==========
26 frequency_array: array_like
27 The frequencies at which we want to calculate the strain
28 mass_1: float
29 The mass of the heavier object in solar masses
30 mass_2: float
31 The mass of the lighter object in solar masses
32 luminosity_distance: float
33 The luminosity distance in megaparsec
34 a_1: float
35 Dimensionless primary spin magnitude
36 tilt_1: float
37 Primary tilt angle
38 phi_12: float
39 Azimuthal angle between the two component spins
40 a_2: float
41 Dimensionless secondary spin magnitude
42 tilt_2: float
43 Secondary tilt angle
44 phi_jl: float
45 Azimuthal angle between the total binary angular momentum and the
46 orbital angular momentum
47 theta_jn: float
48 Angle between the total binary angular momentum and the line of sight
49 phase: float
50 The phase at coalescence
51 kwargs: dict
52 Optional keyword arguments
53 Supported arguments:
55 - waveform_approximant
56 - reference_frequency
57 - minimum_frequency
58 - maximum_frequency
59 - catch_waveform_errors
60 - pn_amplitude_order
61 - mode_array:
62 Activate a specific mode array and evaluate the model using those
63 modes only. e.g. waveform_arguments =
64 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2]])
65 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
66 specify modes that are included in that particular model. e.g.
67 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
68 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
69 55 modes are not included in this model. Be aware that some models
70 only take positive modes and return the positive and the negative
71 mode together, while others need to call both. e.g.
72 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
73 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
74 However, waveform_arguments =
75 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
76 returns the 22 and 4-4 of IMRPhenomXHM.
78 Returns
79 =======
80 dict: A dictionary with the plus and cross polarisation strain modes
82 Notes
83 =====
84 This function is a temporary wrapper to the interface that will
85 likely be significantly changed or removed in a future release.
86 This version is only intended to be used with `SEOBNRv5HM` and `SEOBNRv5PHM` and
87 does not have full functionality for other waveform models.
88 """
90 from lalsimulation.gwsignal import GenerateFDWaveform
91 from lalsimulation.gwsignal.models import gwsignal_get_waveform_generator
92 import astropy.units as u
94 waveform_kwargs = dict(
95 waveform_approximant="SEOBNRv5PHM",
96 reference_frequency=50.0,
97 minimum_frequency=20.0,
98 maximum_frequency=frequency_array[-1],
99 catch_waveform_errors=False,
100 mode_array=None,
101 pn_amplitude_order=0,
102 )
103 waveform_kwargs.update(kwargs)
105 waveform_approximant = waveform_kwargs['waveform_approximant']
106 if waveform_approximant not in ["SEOBNRv5HM", "SEOBNRv5PHM"]:
107 if waveform_approximant == "IMRPhenomXPHM":
108 logger.warning("The new waveform interface is unreviewed for this model" +
109 "and it is only intended for testing.")
110 else:
111 raise ValueError("The new waveform interface is unreviewed for this model.")
112 reference_frequency = waveform_kwargs['reference_frequency']
113 minimum_frequency = waveform_kwargs['minimum_frequency']
114 maximum_frequency = waveform_kwargs['maximum_frequency']
115 catch_waveform_errors = waveform_kwargs['catch_waveform_errors']
116 mode_array = waveform_kwargs['mode_array']
117 pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
119 if pn_amplitude_order != 0:
120 # This is to mimic the behaviour in
121 # https://git.ligo.org/lscsoft/lalsuite/-/blob/master/lalsimulation/lib/LALSimInspiral.c#L5542
122 if pn_amplitude_order == -1:
123 if waveform_approximant in ["SpinTaylorT4", "SpinTaylorT5"]:
124 pn_amplitude_order = 3 # Equivalent to MAX_PRECESSING_AMP_PN_ORDER in LALSimulation
125 else:
126 pn_amplitude_order = 6 # Equivalent to MAX_NONPRECESSING_AMP_PN_ORDER in LALSimulation
127 start_frequency = minimum_frequency * 2. / (pn_amplitude_order + 2)
128 else:
129 start_frequency = minimum_frequency
131 # Call GWsignal generator
132 wf_gen = gwsignal_get_waveform_generator(waveform_approximant)
134 delta_frequency = frequency_array[1] - frequency_array[0]
136 frequency_bounds = ((frequency_array >= minimum_frequency) *
137 (frequency_array <= maximum_frequency))
139 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
140 theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
141 phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1 * utils.solar_mass, mass_2=mass_2 * utils.solar_mass,
142 reference_frequency=reference_frequency, phase=phase)
144 eccentricity = 0.0
145 longitude_ascending_nodes = 0.0
146 mean_per_ano = 0.0
148 # Check if conditioning is needed
149 condition = 0
150 if wf_gen.metadata["implemented_domain"] == 'time':
151 condition = 1
153 # Create dict for gwsignal generator
154 gwsignal_dict = {'mass1' : mass_1 * u.solMass,
155 'mass2' : mass_2 * u.solMass,
156 'spin1x' : spin_1x * u.dimensionless_unscaled,
157 'spin1y' : spin_1y * u.dimensionless_unscaled,
158 'spin1z' : spin_1z * u.dimensionless_unscaled,
159 'spin2x' : spin_2x * u.dimensionless_unscaled,
160 'spin2y' : spin_2y * u.dimensionless_unscaled,
161 'spin2z' : spin_2z * u.dimensionless_unscaled,
162 'deltaF' : delta_frequency * u.Hz,
163 'f22_start' : start_frequency * u.Hz,
164 'f_max': maximum_frequency * u.Hz,
165 'f22_ref': reference_frequency * u.Hz,
166 'phi_ref' : phase * u.rad,
167 'distance' : luminosity_distance * u.Mpc,
168 'inclination' : iota * u.rad,
169 'eccentricity' : eccentricity * u.dimensionless_unscaled,
170 'longAscNodes' : longitude_ascending_nodes * u.rad,
171 'meanPerAno' : mean_per_ano * u.rad,
172 # 'ModeArray': mode_array,
173 'condition': condition
174 }
176 if mode_array is not None:
177 gwsignal_dict.update(ModeArray=mode_array)
179 # Pass extra waveform arguments to gwsignal
180 extra_args = waveform_kwargs.copy()
182 for key in [
183 "waveform_approximant",
184 "reference_frequency",
185 "minimum_frequency",
186 "maximum_frequency",
187 "catch_waveform_errors",
188 "mode_array",
189 "pn_spin_order",
190 "pn_amplitude_order",
191 "pn_tidal_order",
192 "pn_phase_order",
193 "numerical_relativity_file",
194 ]:
195 if key in extra_args.keys():
196 del extra_args[key]
198 gwsignal_dict.update(extra_args)
200 try:
201 hpc = GenerateFDWaveform(gwsignal_dict, wf_gen)
202 except Exception as e:
203 if not catch_waveform_errors:
204 raise
205 else:
206 EDOM = (
207 "Internal function call failed: Input domain error" in e.args[0]
208 ) or "Input domain error" in e.args[
209 0
210 ]
211 if EDOM:
212 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
213 spin_1=(spin_1x, spin_1y, spin_1z),
214 spin_2=(spin_2x, spin_2y, spin_2z),
215 luminosity_distance=luminosity_distance,
216 iota=iota, phase=phase,
217 eccentricity=eccentricity,
218 start_frequency=minimum_frequency)
219 logger.warning("Evaluating the waveform failed with error: {}\n".format(e) +
220 "The parameters were {}\n".format(failed_parameters) +
221 "Likelihood will be set to -inf.")
222 return None
223 else:
224 raise
226 hplus = hpc.hp
227 hcross = hpc.hc
229 h_plus = np.zeros_like(frequency_array, dtype=complex)
230 h_cross = np.zeros_like(frequency_array, dtype=complex)
232 if len(hplus) > len(frequency_array):
233 logger.debug("GWsignal waveform longer than bilby's `frequency_array`" +
234 "({} vs {}), ".format(len(hplus), len(frequency_array)) +
235 "probably because padded with zeros up to the next power of two length." +
236 " Truncating GWsignal array.")
237 h_plus = hplus[:len(h_plus)]
238 h_cross = hcross[:len(h_cross)]
239 else:
240 h_plus[:len(hplus)] = hplus
241 h_cross[:len(hcross)] = hcross
243 h_plus *= frequency_bounds
244 h_cross *= frequency_bounds
246 if condition:
247 dt = 1 / hplus.df.value + hplus.epoch.value
248 time_shift = np.exp(-1j * 2 * np.pi * dt * frequency_array[frequency_bounds])
249 h_plus[frequency_bounds] *= time_shift
250 h_cross[frequency_bounds] *= time_shift
252 return dict(plus=h_plus, cross=h_cross)
255def lal_binary_black_hole(
256 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
257 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
258 """ A Binary Black Hole waveform model using lalsimulation
260 Parameters
261 ==========
262 frequency_array: array_like
263 The frequencies at which we want to calculate the strain
264 mass_1: float
265 The mass of the heavier object in solar masses
266 mass_2: float
267 The mass of the lighter object in solar masses
268 luminosity_distance: float
269 The luminosity distance in megaparsec
270 a_1: float
271 Dimensionless primary spin magnitude
272 tilt_1: float
273 Primary tilt angle
274 phi_12: float
275 Azimuthal angle between the two component spins
276 a_2: float
277 Dimensionless secondary spin magnitude
278 tilt_2: float
279 Secondary tilt angle
280 phi_jl: float
281 Azimuthal angle between the total binary angular momentum and the
282 orbital angular momentum
283 theta_jn: float
284 Angle between the total binary angular momentum and the line of sight
285 phase: float
286 The phase at reference frequency or peak amplitude (depends on waveform)
287 kwargs: dict
288 Optional keyword arguments
289 Supported arguments:
291 - waveform_approximant
292 - reference_frequency
293 - minimum_frequency
294 - maximum_frequency
295 - catch_waveform_errors
296 - pn_spin_order
297 - pn_tidal_order
298 - pn_phase_order
299 - pn_amplitude_order
300 - mode_array:
301 Activate a specific mode array and evaluate the model using those
302 modes only. e.g. waveform_arguments =
303 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
304 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
305 specify modes that are included in that particular model. e.g.
306 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
307 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
308 55 modes are not included in this model. Be aware that some models
309 only take positive modes and return the positive and the negative
310 mode together, while others need to call both. e.g.
311 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
312 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
313 However, waveform_arguments =
314 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
315 returns the 22 and 4-4 of IMRPhenomXHM.
316 - lal_waveform_dictionary:
317 A dictionary (lal.Dict) of arguments passed to the lalsimulation
318 waveform generator. The arguments are specific to the waveform used.
320 Returns
321 =======
322 dict: A dictionary with the plus and cross polarisation strain modes
323 """
324 waveform_kwargs = dict(
325 waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
326 minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
327 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
328 pn_phase_order=-1, pn_amplitude_order=0)
329 waveform_kwargs.update(kwargs)
330 return _base_lal_cbc_fd_waveform(
331 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
332 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
333 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
334 phi_jl=phi_jl, **waveform_kwargs)
337def lal_binary_neutron_star(
338 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
339 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, lambda_1, lambda_2,
340 **kwargs):
341 """ A Binary Neutron Star waveform model using lalsimulation
343 Parameters
344 ==========
345 frequency_array: array_like
346 The frequencies at which we want to calculate the strain
347 mass_1: float
348 The mass of the heavier object in solar masses
349 mass_2: float
350 The mass of the lighter object in solar masses
351 luminosity_distance: float
352 The luminosity distance in megaparsec
353 a_1: float
354 Dimensionless primary spin magnitude
355 tilt_1: float
356 Primary tilt angle
357 phi_12: float
358 Azimuthal angle between the two component spins
359 a_2: float
360 Dimensionless secondary spin magnitude
361 tilt_2: float
362 Secondary tilt angle
363 phi_jl: float
364 Azimuthal angle between the total binary angular momentum and the
365 orbital angular momentum
366 theta_jn: float
367 Orbital inclination
368 phase: float
369 The phase at reference frequency or peak amplitude (depends on waveform)
370 lambda_1: float
371 Dimensionless tidal deformability of mass_1
372 lambda_2: float
373 Dimensionless tidal deformability of mass_2
374 kwargs: dict
375 Optional keyword arguments
376 Supported arguments:
378 - waveform_approximant
379 - reference_frequency
380 - minimum_frequency
381 - maximum_frequency
382 - catch_waveform_errors
383 - pn_spin_order
384 - pn_tidal_order
385 - pn_phase_order
386 - pn_amplitude_order
387 - mode_array:
388 Activate a specific mode array and evaluate the model using those
389 modes only. e.g. waveform_arguments =
390 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
391 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
392 specify modes that are included in that particular model. e.g.
393 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
394 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
395 55 modes are not included in this model. Be aware that some models
396 only take positive modes and return the positive and the negative
397 mode together, while others need to call both. e.g.
398 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
399 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
400 However, waveform_arguments =
401 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
402 returns the 22 and 4-4 of IMRPhenomXHM.
404 Returns
405 =======
406 dict: A dictionary with the plus and cross polarisation strain modes
407 """
408 waveform_kwargs = dict(
409 waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0,
410 minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
411 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
412 pn_phase_order=-1, pn_amplitude_order=0)
413 waveform_kwargs.update(kwargs)
414 return _base_lal_cbc_fd_waveform(
415 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
416 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
417 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
418 phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
421def lal_eccentric_binary_black_hole_no_spins(
422 frequency_array, mass_1, mass_2, eccentricity, luminosity_distance,
423 theta_jn, phase, **kwargs):
424 """ Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
426 Parameters
427 ==========
428 frequency_array: array_like
429 The frequencies at which we want to calculate the strain
430 mass_1: float
431 The mass of the heavier object in solar masses
432 mass_2: float
433 The mass of the lighter object in solar masses
434 eccentricity: float
435 The orbital eccentricity of the system
436 luminosity_distance: float
437 The luminosity distance in megaparsec
438 theta_jn: float
439 Orbital inclination
440 phase: float
441 The phase at reference frequency or peak amplitude (depends on waveform)
442 kwargs: dict
443 Optional keyword arguments
444 Supported arguments:
446 - waveform_approximant
447 - reference_frequency
448 - minimum_frequency
449 - maximum_frequency
450 - catch_waveform_errors
451 - pn_spin_order
452 - pn_tidal_order
453 - pn_phase_order
454 - pn_amplitude_order
455 - mode_array:
456 Activate a specific mode array and evaluate the model using those
457 modes only. e.g. waveform_arguments =
458 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
459 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
460 specify modes that are included in that particular model. e.g.
461 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
462 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
463 55 modes are not included in this model. Be aware that some models
464 only take positive modes and return the positive and the negative
465 mode together, while others need to call both. e.g.
466 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
467 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
468 However, waveform_arguments =
469 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
470 returns the 22 and 4-4 of IMRPhenomXHM.
472 Returns
473 =======
474 dict: A dictionary with the plus and cross polarisation strain modes
475 """
476 waveform_kwargs = dict(
477 waveform_approximant='EccentricFD', reference_frequency=10.0,
478 minimum_frequency=10.0, maximum_frequency=frequency_array[-1],
479 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
480 pn_phase_order=-1, pn_amplitude_order=0)
481 waveform_kwargs.update(kwargs)
482 return _base_lal_cbc_fd_waveform(
483 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
484 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
485 eccentricity=eccentricity, **waveform_kwargs)
488def set_waveform_dictionary(waveform_kwargs, lambda_1=0, lambda_2=0):
489 """
490 Add keyword arguments to the :code:`LALDict` object.
492 Parameters
493 ==========
494 waveform_kwargs: dict
495 A dictionary of waveform kwargs. This is modified in place to remove used arguments.
496 lambda_1: float
497 Dimensionless tidal deformability of the primary object.
498 lambda_2: float
499 Dimensionless tidal deformability of the primary object.
501 Returns
502 =======
503 waveform_dictionary: lal.LALDict
504 The lal waveform dictionary. This is either taken from the waveform_kwargs or created
505 internally.
506 """
507 import lalsimulation as lalsim
508 from lal import CreateDict
509 waveform_dictionary = waveform_kwargs.pop('lal_waveform_dictionary', CreateDict())
510 waveform_kwargs["TidalLambda1"] = lambda_1
511 waveform_kwargs["TidalLambda2"] = lambda_2
512 waveform_kwargs["NumRelData"] = waveform_kwargs.pop("numerical_relativity_data", None)
514 for key in [
515 "pn_spin_order", "pn_tidal_order", "pn_phase_order", "pn_amplitude_order"
516 ]:
517 waveform_kwargs[key[:2].upper() + key[3:].title().replace('_', '')] = waveform_kwargs.pop(key)
519 for key in list(waveform_kwargs.keys()).copy():
520 func = getattr(lalsim, f"SimInspiralWaveformParamsInsert{key}", None)
521 if func is None:
522 continue
523 value = waveform_kwargs.pop(key)
524 if func is not None and value is not None:
525 func(waveform_dictionary, value)
527 mode_array = waveform_kwargs.pop("mode_array", None)
528 if mode_array is not None:
529 mode_array_lal = lalsim.SimInspiralCreateModeArray()
530 for mode in mode_array:
531 lalsim.SimInspiralModeArrayActivateMode(mode_array_lal, mode[0], mode[1])
532 lalsim.SimInspiralWaveformParamsInsertModeArray(waveform_dictionary, mode_array_lal)
533 return waveform_dictionary
536def _base_lal_cbc_fd_waveform(
537 frequency_array, mass_1, mass_2, luminosity_distance, theta_jn, phase,
538 a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0,
539 lambda_1=0.0, lambda_2=0.0, eccentricity=0.0, **waveform_kwargs):
540 """ Generate a cbc waveform model using lalsimulation
542 Parameters
543 ==========
544 frequency_array: array_like
545 The frequencies at which we want to calculate the strain
546 mass_1: float
547 The mass of the heavier object in solar masses
548 mass_2: float
549 The mass of the lighter object in solar masses
550 luminosity_distance: float
551 The luminosity distance in megaparsec
552 a_1: float
553 Dimensionless primary spin magnitude
554 tilt_1: float
555 Primary tilt angle
556 phi_12: float
557 Azimuthal angle between the component spins
558 a_2: float
559 Dimensionless secondary spin magnitude
560 tilt_2: float
561 Secondary tilt angle
562 phi_jl: float
563 Azimuthal angle between the total and orbital angular momenta
564 theta_jn: float
565 Orbital inclination
566 phase: float
567 The phase at reference frequency or peak amplitude (depends on waveform)
568 eccentricity: float
569 Binary eccentricity
570 lambda_1: float
571 Tidal deformability of the more massive object
572 lambda_2: float
573 Tidal deformability of the less massive object
574 kwargs: dict
575 Optional keyword arguments
577 Returns
578 =======
579 dict: A dictionary with the plus and cross polarisation strain modes
580 """
581 import lalsimulation as lalsim
583 waveform_approximant = waveform_kwargs.pop('waveform_approximant')
584 reference_frequency = waveform_kwargs.pop('reference_frequency')
585 minimum_frequency = waveform_kwargs.pop('minimum_frequency')
586 maximum_frequency = waveform_kwargs.pop('maximum_frequency')
587 catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors')
588 pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
590 waveform_dictionary = set_waveform_dictionary(waveform_kwargs, lambda_1, lambda_2)
591 approximant = lalsim_GetApproximantFromString(waveform_approximant)
593 if pn_amplitude_order != 0:
594 start_frequency = lalsim.SimInspiralfLow2fStart(
595 float(minimum_frequency), int(pn_amplitude_order), approximant
596 )
597 else:
598 start_frequency = minimum_frequency
600 delta_frequency = frequency_array[1] - frequency_array[0]
602 frequency_bounds = ((frequency_array >= minimum_frequency) *
603 (frequency_array <= maximum_frequency))
605 luminosity_distance = luminosity_distance * 1e6 * utils.parsec
606 mass_1 = mass_1 * utils.solar_mass
607 mass_2 = mass_2 * utils.solar_mass
609 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
610 theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
611 phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2,
612 reference_frequency=reference_frequency, phase=phase)
614 longitude_ascending_nodes = 0.0
615 mean_per_ano = 0.0
617 if lalsim.SimInspiralImplementedFDApproximants(approximant):
618 wf_func = lalsim_SimInspiralChooseFDWaveform
619 else:
620 wf_func = lalsim_SimInspiralFD
621 try:
622 hplus, hcross = wf_func(
623 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
624 spin_2z, luminosity_distance, iota, phase,
625 longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
626 start_frequency, maximum_frequency, reference_frequency,
627 waveform_dictionary, approximant)
628 except Exception as e:
629 if not catch_waveform_errors:
630 raise
631 else:
632 EDOM = (e.args[0] == 'Internal function call failed: Input domain error')
633 if EDOM:
634 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
635 spin_1=(spin_1x, spin_1y, spin_1z),
636 spin_2=(spin_2x, spin_2y, spin_2z),
637 luminosity_distance=luminosity_distance,
638 iota=iota, phase=phase,
639 eccentricity=eccentricity,
640 start_frequency=start_frequency)
641 logger.warning("Evaluating the waveform failed with error: {}\n".format(e) +
642 "The parameters were {}\n".format(failed_parameters) +
643 "Likelihood will be set to -inf.")
644 return None
645 else:
646 raise
648 h_plus = np.zeros_like(frequency_array, dtype=complex)
649 h_cross = np.zeros_like(frequency_array, dtype=complex)
651 if len(hplus.data.data) > len(frequency_array):
652 logger.debug("LALsim waveform longer than bilby's `frequency_array`" +
653 "({} vs {}), ".format(len(hplus.data.data), len(frequency_array)) +
654 "probably because padded with zeros up to the next power of two length." +
655 " Truncating lalsim array.")
656 h_plus = hplus.data.data[:len(h_plus)]
657 h_cross = hcross.data.data[:len(h_cross)]
658 else:
659 h_plus[:len(hplus.data.data)] = hplus.data.data
660 h_cross[:len(hcross.data.data)] = hcross.data.data
662 h_plus *= frequency_bounds
663 h_cross *= frequency_bounds
665 if wf_func == lalsim_SimInspiralFD:
666 dt = 1 / hplus.deltaF + (hplus.epoch.gpsSeconds + hplus.epoch.gpsNanoSeconds * 1e-9)
667 time_shift = np.exp(-1j * 2 * np.pi * dt * frequency_array[frequency_bounds])
668 h_plus[frequency_bounds] *= time_shift
669 h_cross[frequency_bounds] *= time_shift
671 if len(waveform_kwargs) > 0:
672 logger.warning(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs))
674 return dict(plus=h_plus, cross=h_cross)
677def binary_black_hole_roq(
678 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
679 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
680 waveform_kwargs = dict(
681 waveform_approximant='IMRPhenomPv2', reference_frequency=20.0,
682 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
683 pn_phase_order=-1, pn_amplitude_order=0)
684 waveform_kwargs.update(waveform_arguments)
685 return _base_roq_waveform(
686 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
687 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
688 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
689 phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
692def binary_neutron_star_roq(
693 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
694 phi_12, a_2, tilt_2, phi_jl, lambda_1, lambda_2, theta_jn, phase,
695 **waveform_arguments):
696 waveform_kwargs = dict(
697 waveform_approximant='IMRPhenomD_NRTidal', reference_frequency=20.0,
698 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
699 pn_phase_order=-1, pn_amplitude_order=0)
700 waveform_kwargs.update(waveform_arguments)
701 return _base_roq_waveform(
702 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
703 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
704 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
705 phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
708def lal_binary_black_hole_relative_binning(
709 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
710 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, fiducial, **kwargs):
711 """ Source model to go with RelativeBinningGravitationalWaveTransient likelihood.
713 All parameters are the same as in the usual source models, except `fiducial`
715 fiducial: float
716 If fiducial=1, waveform evaluated on the full frequency grid is returned.
717 If fiducial=0, waveform evaluated at waveform_kwargs["frequency_bin_edges"]
718 is returned.
719 """
721 waveform_kwargs = dict(
722 waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
723 minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
724 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
725 pn_phase_order=-1, pn_amplitude_order=0)
726 waveform_kwargs.update(kwargs)
728 if fiducial == 1:
729 _ = waveform_kwargs.pop("frequency_bin_edges", None)
730 return _base_lal_cbc_fd_waveform(
731 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
732 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
733 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
734 phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
736 else:
737 _ = waveform_kwargs.pop("minimum_frequency", None)
738 _ = waveform_kwargs.pop("maximum_frequency", None)
739 waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
740 return _base_waveform_frequency_sequence(
741 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
742 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
743 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
744 phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
747def lal_binary_neutron_star_relative_binning(
748 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
749 phi_12, a_2, tilt_2, phi_jl, lambda_1, lambda_2, theta_jn, phase,
750 fiducial, **kwargs):
751 """ Source model to go with RelativeBinningGravitationalWaveTransient likelihood.
753 All parameters are the same as in the usual source models, except `fiducial`
755 fiducial: float
756 If fiducial=1, waveform evaluated on the full frequency grid is returned.
757 If fiducial=0, waveform evaluated at waveform_kwargs["frequency_bin_edges"]
758 is returned.
759 """
761 waveform_kwargs = dict(
762 waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0,
763 minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
764 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
765 pn_phase_order=-1, pn_amplitude_order=0)
766 waveform_kwargs.update(kwargs)
768 if fiducial == 1:
769 return _base_lal_cbc_fd_waveform(
770 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
771 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
772 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
773 phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
774 else:
775 _ = waveform_kwargs.pop("minimum_frequency", None)
776 _ = waveform_kwargs.pop("maximum_frequency", None)
777 waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
778 return _base_waveform_frequency_sequence(
779 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
780 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
781 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
782 phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
785def _base_roq_waveform(
786 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
787 phi_12, a_2, tilt_2, lambda_1, lambda_2, phi_jl, theta_jn, phase,
788 **waveform_arguments):
789 """ Base source model for ROQGravitationalWaveTransient, which evaluates
790 waveform values at frequency nodes contained in waveform_arguments. This
791 requires that waveform_arguments contain all of 'frequency_nodes',
792 'linear_indices', and 'quadratic_indices', or both 'frequency_nodes_linear' and
793 'frequency_nodes_quadratic'.
795 Parameters
796 ==========
797 frequency_array: np.array
798 This input is ignored for the roq source model
799 mass_1: float
800 The mass of the heavier object in solar masses
801 mass_2: float
802 The mass of the lighter object in solar masses
803 luminosity_distance: float
804 The luminosity distance in megaparsec
805 a_1: float
806 Dimensionless primary spin magnitude
807 tilt_1: float
808 Primary tilt angle
809 phi_12: float
811 a_2: float
812 Dimensionless secondary spin magnitude
813 tilt_2: float
814 Secondary tilt angle
815 phi_jl: float
817 theta_jn: float
818 Orbital inclination
819 phase: float
820 The phase at reference frequency or peak amplitude (depends on waveform)
822 Waveform arguments
823 ===================
824 Non-sampled extra data used in the source model calculation
825 frequency_nodes_linear: np.array
826 frequency nodes for linear likelihood part
827 frequency_nodes_quadratic: np.array
828 frequency nodes for quadratic likelihood part
829 frequency_nodes: np.array
830 unique frequency nodes for linear and quadratic likelihood parts
831 linear_indices: np.array
832 indices to recover frequency nodes for linear part from unique frequency nodes
833 quadratic_indices: np.array
834 indices to recover frequency nodes for quadratic part from unique frequency nodes
835 reference_frequency: float
836 approximant: str
838 Note: for the frequency_nodes_linear and frequency_nodes_quadratic arguments,
839 if using data from https://git.ligo.org/lscsoft/ROQ_data, this should be
840 loaded as `np.load(filename).T`.
842 Returns
843 =======
844 waveform_polarizations: dict
845 Dict containing plus and cross modes evaluated at the linear and
846 quadratic frequency nodes.
847 """
848 if 'frequency_nodes' not in waveform_arguments:
849 size_linear = len(waveform_arguments['frequency_nodes_linear'])
850 frequency_nodes_combined = np.hstack(
851 (waveform_arguments.pop('frequency_nodes_linear'),
852 waveform_arguments.pop('frequency_nodes_quadratic'))
853 )
854 frequency_nodes_unique, original_indices = np.unique(
855 frequency_nodes_combined, return_inverse=True
856 )
857 linear_indices = original_indices[:size_linear]
858 quadratic_indices = original_indices[size_linear:]
859 waveform_arguments['frequencies'] = frequency_nodes_unique
860 else:
861 linear_indices = waveform_arguments.pop("linear_indices")
862 quadratic_indices = waveform_arguments.pop("quadratic_indices")
863 for key in ["frequency_nodes_linear", "frequency_nodes_quadratic"]:
864 _ = waveform_arguments.pop(key, None)
865 waveform_arguments['frequencies'] = waveform_arguments.pop('frequency_nodes')
866 waveform_polarizations = _base_waveform_frequency_sequence(
867 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
868 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
869 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
870 phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_arguments)
872 return {
873 'linear': {
874 'plus': waveform_polarizations['plus'][linear_indices],
875 'cross': waveform_polarizations['cross'][linear_indices]
876 },
877 'quadratic': {
878 'plus': waveform_polarizations['plus'][quadratic_indices],
879 'cross': waveform_polarizations['cross'][quadratic_indices]
880 }
881 }
884def binary_black_hole_frequency_sequence(
885 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
886 phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
887 """ A Binary Black Hole waveform model using lalsimulation. This generates
888 a waveform only on specified frequency points. This is useful for
889 likelihood requiring waveform values at a subset of all the frequency
890 samples. For example, this is used for MBGravitationalWaveTransient.
892 Parameters
893 ==========
894 frequency_array: array_like
895 The input is ignored.
896 mass_1: float
897 The mass of the heavier object in solar masses
898 mass_2: float
899 The mass of the lighter object in solar masses
900 luminosity_distance: float
901 The luminosity distance in megaparsec
902 a_1: float
903 Dimensionless primary spin magnitude
904 tilt_1: float
905 Primary tilt angle
906 phi_12: float
907 Azimuthal angle between the two component spins
908 a_2: float
909 Dimensionless secondary spin magnitude
910 tilt_2: float
911 Secondary tilt angle
912 phi_jl: float
913 Azimuthal angle between the total binary angular momentum and the
914 orbital angular momentum
915 theta_jn: float
916 Angle between the total binary angular momentum and the line of sight
917 phase: float
918 The phase at reference frequency or peak amplitude (depends on waveform)
919 kwargs: dict
920 Required keyword arguments
921 - frequencies:
922 ndarray of frequencies at which waveforms are evaluated
924 Optional keyword arguments
925 - waveform_approximant
926 - reference_frequency
927 - catch_waveform_errors
928 - pn_spin_order
929 - pn_tidal_order
930 - pn_phase_order
931 - pn_amplitude_order
932 - mode_array:
933 Activate a specific mode array and evaluate the model using those
934 modes only. e.g. waveform_arguments =
935 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
936 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
937 specify modes that are included in that particular model. e.g.
938 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
939 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
940 55 modes are not included in this model. Be aware that some models
941 only take positive modes and return the positive and the negative
942 mode together, while others need to call both. e.g.
943 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
944 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
945 However, waveform_arguments =
946 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
947 returns the 22 and 4-4 of IMRPhenomXHM.
949 Returns
950 =======
951 dict: A dictionary with the plus and cross polarisation strain modes
952 """
953 waveform_kwargs = dict(
954 waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
955 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
956 pn_phase_order=-1, pn_amplitude_order=0)
957 waveform_kwargs.update(kwargs)
958 return _base_waveform_frequency_sequence(
959 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
960 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
961 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
962 phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
965def binary_neutron_star_frequency_sequence(
966 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
967 phi_12, a_2, tilt_2, phi_jl, lambda_1, lambda_2, theta_jn, phase,
968 **kwargs):
969 """ A Binary Neutron Star waveform model using lalsimulation. This generates
970 a waveform only on specified frequency points. This is useful for
971 likelihood requiring waveform values at a subset of all the frequency
972 samples. For example, this is used for MBGravitationalWaveTransient.
974 Parameters
975 ==========
976 frequency_array: array_like
977 The input is ignored.
978 mass_1: float
979 The mass of the heavier object in solar masses
980 mass_2: float
981 The mass of the lighter object in solar masses
982 luminosity_distance: float
983 The luminosity distance in megaparsec
984 a_1: float
985 Dimensionless primary spin magnitude
986 tilt_1: float
987 Primary tilt angle
988 phi_12: float
989 Azimuthal angle between the two component spins
990 a_2: float
991 Dimensionless secondary spin magnitude
992 tilt_2: float
993 Secondary tilt angle
994 phi_jl: float
995 Azimuthal angle between the total binary angular momentum and the
996 orbital angular momentum
997 lambda_1: float
998 Dimensionless tidal deformability of mass_1
999 lambda_2: float
1000 Dimensionless tidal deformability of mass_2
1001 theta_jn: float
1002 Angle between the total binary angular momentum and the line of sight
1003 phase: float
1004 The phase at reference frequency or peak amplitude (depends on waveform)
1005 kwargs: dict
1006 Required keyword arguments
1007 - frequencies:
1008 ndarray of frequencies at which waveforms are evaluated
1010 Optional keyword arguments
1011 - waveform_approximant
1012 - reference_frequency
1013 - catch_waveform_errors
1014 - pn_spin_order
1015 - pn_tidal_order
1016 - pn_phase_order
1017 - pn_amplitude_order
1018 - mode_array:
1019 Activate a specific mode array and evaluate the model using those
1020 modes only. e.g. waveform_arguments =
1021 dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
1022 returns the 22 and 2-2 modes only of IMRPhenomHM. You can only
1023 specify modes that are included in that particular model. e.g.
1024 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
1025 mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
1026 55 modes are not included in this model. Be aware that some models
1027 only take positive modes and return the positive and the negative
1028 mode together, while others need to call both. e.g.
1029 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
1030 mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
1031 However, waveform_arguments =
1032 dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
1033 returns the 22 and 4-4 of IMRPhenomXHM.
1035 Returns
1036 =======
1037 dict: A dictionary with the plus and cross polarisation strain modes
1038 """
1039 waveform_kwargs = dict(
1040 waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0,
1041 catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
1042 pn_phase_order=-1, pn_amplitude_order=0)
1043 waveform_kwargs.update(kwargs)
1044 return _base_waveform_frequency_sequence(
1045 frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
1046 luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
1047 a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
1048 phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
1051def _base_waveform_frequency_sequence(
1052 frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
1053 phi_12, a_2, tilt_2, lambda_1, lambda_2, phi_jl, theta_jn, phase,
1054 **waveform_kwargs):
1055 """ Generate a cbc waveform model on specified frequency samples
1057 Parameters
1058 ----------
1059 frequency_array: np.array
1060 This input is ignored
1061 mass_1: float
1062 The mass of the heavier object in solar masses
1063 mass_2: float
1064 The mass of the lighter object in solar masses
1065 luminosity_distance: float
1066 The luminosity distance in megaparsec
1067 a_1: float
1068 Dimensionless primary spin magnitude
1069 tilt_1: float
1070 Primary tilt angle
1071 phi_12: float
1072 a_2: float
1073 Dimensionless secondary spin magnitude
1074 tilt_2: float
1075 Secondary tilt angle
1076 phi_jl: float
1077 theta_jn: float
1078 Orbital inclination
1079 phase: float
1080 The phase at reference frequency or peak amplitude (depends on waveform)
1081 waveform_kwargs: dict
1082 Optional keyword arguments
1084 Returns
1085 -------
1086 waveform_polarizations: dict
1087 Dict containing plus and cross modes evaluated at the linear and
1088 quadratic frequency nodes.
1089 """
1090 frequencies = waveform_kwargs.pop('frequencies')
1091 reference_frequency = waveform_kwargs.pop('reference_frequency')
1092 approximant = waveform_kwargs.pop('waveform_approximant')
1093 catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors')
1095 waveform_dictionary = set_waveform_dictionary(waveform_kwargs, lambda_1, lambda_2)
1096 approximant = lalsim_GetApproximantFromString(approximant)
1098 luminosity_distance = luminosity_distance * 1e6 * utils.parsec
1099 mass_1 = mass_1 * utils.solar_mass
1100 mass_2 = mass_2 * utils.solar_mass
1102 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
1103 theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
1104 phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2,
1105 reference_frequency=reference_frequency, phase=phase)
1107 try:
1108 h_plus, h_cross = lalsim_SimInspiralChooseFDWaveformSequence(
1109 phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
1110 spin_2z, reference_frequency, luminosity_distance, iota,
1111 waveform_dictionary, approximant, frequencies)
1112 except Exception as e:
1113 if not catch_waveform_errors:
1114 raise
1115 else:
1116 EDOM = (e.args[0] == 'Internal function call failed: Input domain error')
1117 if EDOM:
1118 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
1119 spin_1=(spin_1x, spin_1y, spin_1z),
1120 spin_2=(spin_2x, spin_2y, spin_2z),
1121 luminosity_distance=luminosity_distance,
1122 iota=iota, phase=phase)
1123 logger.warning("Evaluating the waveform failed with error: {}\n".format(e) +
1124 "The parameters were {}\n".format(failed_parameters) +
1125 "Likelihood will be set to -inf.")
1126 return None
1127 else:
1128 raise
1130 if len(waveform_kwargs) > 0:
1131 logger.warning(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs))
1133 return dict(plus=h_plus.data.data, cross=h_cross.data.data)
1136def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
1137 r"""
1138 A frequency-domain sine-Gaussian burst source model.
1140 .. math::
1142 \tau &= \frac{Q}{\sqrt{2}\pi f_{0}} \\
1143 \alpha &= \frac{Q}{4\sqrt{\pi} f_{0}} \\
1144 h_{+} &=
1145 \frac{h_{\rm rss}\sqrt{\pi}\tau}{2\sqrt{\alpha (1 + e^{-Q^2})}}
1146 \left(
1147 e^{-\pi^2 \tau^2 (f + f_{0})^2}
1148 + e^{-\pi^2 \tau^2 (f - f_{0})^2}
1149 \right) \\
1150 h_{\times} &=
1151 \frac{i h_{\rm rss}\sqrt{\pi}\tau}{2\sqrt{\alpha (1 - e^{-Q^2})}}
1152 \left(
1153 e^{-\pi^2 \tau^2 (f + f_{0})^2}
1154 - e^{-\pi^2 \tau^2 (f - f_{0})^2}
1155 \right)
1157 Parameters
1158 ----------
1159 frequency_array: array-like
1160 The frequencies at which to compute the model.
1161 hrss: float
1162 The amplitude of the signal.
1163 Q: float
1164 The quality factor of the burst, determines the decay time.
1165 frequency: float
1166 The peak frequency of the burst.
1167 kwargs: dict
1168 UNUSED
1170 Returns
1171 -------
1172 dict:
1173 Dictionary containing the plus and cross components of the strain.
1174 """
1175 tau = Q / (np.sqrt(2.0) * np.pi * frequency)
1176 temp = Q / (4.0 * np.sqrt(np.pi) * frequency)
1177 fm = frequency_array - frequency
1178 fp = frequency_array + frequency
1180 h_plus = ((hrss / np.sqrt(temp * (1 + np.exp(-Q**2)))) *
1181 ((np.sqrt(np.pi) * tau) / 2.0) *
1182 (np.exp(-fm**2 * np.pi**2 * tau**2) +
1183 np.exp(-fp**2 * np.pi**2 * tau**2)))
1185 h_cross = (-1j * (hrss / np.sqrt(temp * (1 - np.exp(-Q**2)))) *
1186 ((np.sqrt(np.pi) * tau) / 2.0) *
1187 (np.exp(-fm**2 * np.pi**2 * tau**2) -
1188 np.exp(-fp**2 * np.pi**2 * tau**2)))
1190 return {'plus': h_plus, 'cross': h_cross}
1193def supernova(frequency_array, luminosity_distance, **kwargs):
1194 """
1195 A source model that reads a simulation from a text file.
1197 This was originally intended for use with supernova simulations, but can
1198 be applied to any source class.
1200 Parameters
1201 ----------
1202 frequency_array: array-like
1203 Unused
1204 file_path: str
1205 Path to the file containing the NR simulation. The format of this file
1206 should be readable by :code:`numpy.loadtxt` and have four columns
1207 containing the real and imaginary components of the plus and cross
1208 polarizations.
1209 luminosity_distance: float
1210 The distance to the source in kpc, this scales the amplitude of the
1211 signal. The simulation is assumed to be at 10kpc.
1212 kwargs:
1213 extra keyword arguments, this should include the :code:`file_path`
1215 Returns
1216 -------
1217 dict:
1218 A dictionary containing the plus and cross components of the signal.
1219 """
1221 file_path = kwargs["file_path"]
1222 data = np.genfromtxt(file_path)
1224 # waveform in file at 10kpc
1225 scaling = 1e-3 * (10.0 / luminosity_distance)
1227 h_plus = scaling * (data[:, 0] + 1j * data[:, 1])
1228 h_cross = scaling * (data[:, 2] + 1j * data[:, 3])
1229 return {'plus': h_plus, 'cross': h_cross}
1232def supernova_pca_model(
1233 frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5, luminosity_distance, **kwargs
1234):
1235 r"""
1236 Signal model based on a five-component principal component decomposition
1237 of a model.
1239 While this was initially intended for modelling supernova signal, it is
1240 applicable to any situation using such a principal component decomposition.
1242 .. math::
1244 h_{A} = \frac{10^{-22}}{d_{L}} \sum_{i=1}^{5} c_{i} h_{i}
1246 Parameters
1247 ----------
1248 frequency_array: UNUSED
1249 pc_coeff1: float
1250 The first principal component coefficient.
1251 pc_coeff2: float
1252 The second principal component coefficient.
1253 pc_coeff3: float
1254 The third principal component coefficient.
1255 pc_coeff4: float
1256 The fourth principal component coefficient.
1257 pc_coeff5: float
1258 The fifth principal component coefficient.
1259 luminosity_distance: float
1260 The distance to the source, the amplitude is scaled such that the
1261 amplitude at 10 kpc is 1e-23.
1262 kwargs: dict
1263 Dictionary containing numpy arrays with the real and imaginary
1264 components of the principal component decomposition.
1266 Returns
1267 -------
1268 dict:
1269 The plus and cross polarizations of the signal
1270 """
1272 principal_components = kwargs["realPCs"] + 1j * kwargs["imagPCs"]
1273 coefficients = [pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5]
1275 strain = np.sum(
1276 [coeff * principal_components[:, ii] for ii, coeff in enumerate(coefficients)],
1277 axis=0
1278 )
1280 # file at 10kpc
1281 scaling = 1e-23 * (10.0 / luminosity_distance)
1283 h_plus = scaling * strain
1284 h_cross = scaling * strain
1286 return {'plus': h_plus, 'cross': h_cross}
1289precession_only = {
1290 "tilt_1", "tilt_2", "phi_12", "phi_jl", "chi_1_in_plane", "chi_2_in_plane",
1291}
1293spin = {
1294 "a_1", "a_2", "tilt_1", "tilt_2", "phi_12", "phi_jl", "chi_1", "chi_2",
1295 "chi_1_in_plane", "chi_2_in_plane",
1296}
1297mass = {
1298 "chirp_mass", "mass_ratio", "total_mass", "mass_1", "mass_2",
1299 "symmetric_mass_ratio",
1300}
1301primary_spin_and_q = {
1302 "a_1", "chi_1", "mass_ratio"
1303}
1304tidal = {
1305 "lambda_1", "lambda_2", "lambda_tilde", "delta_lambda_tilde"
1306}
1307phase = {
1308 "phase", "delta_phase",
1309}
1310extrinsic = {
1311 "azimuth", "zenith", "luminosity_distance", "psi", "theta_jn",
1312 "cos_theta_jn", "geocent_time", "time_jitter", "ra", "dec",
1313 "H1_time", "L1_time", "V1_time",
1314}
1315sky = {
1316 "azimuth", "zenith", "ra", "dec",
1317}
1318distance_inclination = {
1319 "luminosity_distance", "redshift", "theta_jn", "cos_theta_jn",
1320}
1321measured_spin = {
1322 "chi_1", "chi_2", "a_1", "a_2", "chi_1_in_plane"
1323}
1325PARAMETER_SETS = dict(
1326 spin=spin, mass=mass, phase=phase, extrinsic=extrinsic,
1327 tidal=tidal, primary_spin_and_q=primary_spin_and_q,
1328 intrinsic=spin.union(mass).union(phase).union(tidal),
1329 precession_only=precession_only,
1330 sky=sky, distance_inclination=distance_inclination,
1331 measured_spin=measured_spin,
1332)