Source code for kuibit.gw_utils

#!/usr/bin/env python3

# Copyright (C) 2020-2024 Gabriele Bozzola
# sYlm: Copyright (C) Christian Reisswig
# ra_dec_to_theta_phi: Copyright (C) 2020-2024 Yuk Tung Liu
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, see <https://www.gnu.org/licenses/>.


"""The module :py:mod:`~.gw_utils` contains convenience functions and
structures to analyze and work with gravitational waves.

First, the ``Detectors`` object is defined. ``Detectors`` is a named tuple
with fields "hanford", "livingston", "virgo". This is used every time we
deal with specific detectors.

The functions provided are:

- :py:func:`~.luminosity_distance_to_redshift`: convert a given luminosity
  distance to a redshift in the ΛCDM cosmology.
- :py:func:`~.sYlm`: return the spin-weighted spherical harmonics at a given
  angle.
- :py:func:`~.ra_dec_to_theta_phi`: convert right ascension and declination to
  spherical coordinates.
- :py:func:`~.antenna_responses`: compute the antenna responses of a given
  angle.
- :py:func:`~.antenna_responses_from_sky_localization`: compute the antenna
  responses for known detectors at a given sky localization.
- :py:func:`~.signal_to_noise_ratio_from_strain`: compute the signal to noise
  for a given signal and a given noise curve.
- :py:func:`~.effective_amplitude_spectral_density`: compute an effective (
  polarization-averaged) amplitude spectral density (ASD) from a strain signal.

"""

import datetime
from collections import namedtuple

import numpy as np
from scipy import integrate, optimize

import kuibit.timeseries as ts
import kuibit.unitconv as uc

# This is just a convenience to avoid having to remember the order of
# the output (and for easy of extension)
# One can access the output as detectos.hanford, or
# detectors['hanford']
Detectors = namedtuple("Detectors", "hanford livingston virgo")


[docs]def luminosity_distance_to_redshift( luminosity_distance, Omega_m=0.309, Omega_L=0.691, initial_guess=0.1, ): r"""Compute redshift from luminosity distance in Mpc assuming the ΛCDM cosmology. This function is useful to correctly reproduce observed signals from cosmological sources (e.g., binary black holes far away). The redshift is computed via root-finding, so an initial guess is needed. :param luminosity_distance: Luminosity distance in megaparsec. :type luminosity_distance: float :param Omega_m: :math:`\Omega_m` (matter) cosmological parameter. :type Omega_m: float :param Omega_L: :math:`\Omega_m` (dark energy) cosmological parameter. :type Omega_L: float :param initial_guess: Initial guess to the redshift for the root-finding routine. :type initial_guess: float :returns z: Redshift :rtype z: float """ distance_in_m = luminosity_distance * uc.MEGAPARSEC_SI # m H0 = uc.H0_SI # 1/s c = uc.C_SI # m/s def function_to_root_find(z): def z_to_DL(z): def DL_integral(z): return 1 / np.sqrt(Omega_m * (1 + z) ** 3 + Omega_L) return c / H0 * (1 + z) * integrate.quad(DL_integral, 0, z)[0] return np.abs(distance_in_m - z_to_DL(z)) redshift = optimize.root(function_to_root_find, initial_guess) if redshift["status"] != 1: raise RuntimeError("Conversion between distance and redshift failed!") return redshift["x"][0]
[docs]def sYlm(ss, ll, mm, theta, phi): """Compute spin-weighted spherical harmonics at the angles ``theta`` and ``phi``. When ``ss = 0``, these are spherical harmonics. :param ss: Spin weight. :type ss: int :param ll: :math:`l` multipolar number. :type ll: int :param mm: :math:`m` multipolar number. :type mm: int :param theta: Meridional angle. :type theta: float :param phi: Azimuthal angle. :type phi: float :returns sYlm: Spin-weighted spherical harmonic evaluated at ``theta`` and ``phi``. :rtype sYlm: float """ # Code by Christian Reisswig # Recursion function for spin-weighted spherical harmonics def s_lambda_lm(local_s, local_l, local_m, x): # Coefficient function for spin-weighted spherical harmonics def sYlm_Cslm(local_s, local_l, local_m): return np.sqrt( local_l * local_l * (4.0 * local_l * local_l - 1.0) / ( (local_l * local_l - local_m * local_m) * (local_l * local_l - local_s * local_s) ) ) Pm = np.power(-0.5, local_m) if local_m != local_s: Pm = Pm * np.power(1.0 + x, (local_m - local_s) * 1.0 / 2) if local_m != -local_s: Pm = Pm * np.power(1.0 - x, (local_m + local_s) * 1.0 / 2) Pm = Pm * np.sqrt( np.math.factorial(2 * local_m + 1) * 1.0 / ( 4.0 * np.pi * np.math.factorial(local_m + local_s) * np.math.factorial(local_m - local_s) ) ) if local_l == local_m: return Pm Pm1 = ( (x + local_s * 1.0 / (local_m + 1)) * sYlm_Cslm(local_s, local_m + 1, local_m) * Pm ) if local_l == local_m + 1: return Pm1 for n in range(local_m + 2, local_l + 1): Pn = (x + local_s * local_m * 1.0 / (n * (n - 1.0))) * sYlm_Cslm( local_s, n, local_m ) * Pm1 - sYlm_Cslm(local_s, n, local_m) * 1.0 / sYlm_Cslm( local_s, n - 1, local_m ) * Pm Pm = Pm1 Pm1 = Pn return Pn Pm = 1.0 mult_l = ll mult_m = mm mult_s = ss if mult_l < 0: return 0 if abs(mult_m) > mult_l or mult_l < abs(mult_s): return 0 if abs(mm) < abs(ss): mult_s = mm mult_m = ss if (mult_m + mult_s) % 2: Pm = -Pm if mult_m < 0: mult_s = -mult_s mult_m = -mult_m if (mult_m + mult_s) % 2: Pm = -Pm result = Pm * s_lambda_lm(mult_s, mult_l, mult_m, np.cos(theta)) return complex(result * np.cos(mm * phi), result * np.sin(mm * phi))
[docs]def ra_dec_to_theta_phi(right_ascension, declination, time_utc): """Compute the spherical angles ``theta`` and ``phi`` for Hanford, Livingston and Virgo for a given source localization. ``utc_time`` has to have the following formatting: ``%Y-%m-%d %H:%M``, (eg ``2015-09-14 09:50:45``) :param right_ascension: Right ascension of the source in degrees. :type right_ascension: float :param declination: Declination of the source in degrees. :type declination: float :param time_utc: UTC time of the event. :type declination: str :returns spherical coordinates: ``Theta``, ``phi`` for the different detectors. :rtype: namedtuple with fields hanford, livingston, and virgo """ # NOTE: This function can make use of some cleanup and double-checking # This code is based on the one written by Yuk Tung Liu. # RA and DEC of the GW source: alpha_hours = right_ascension delta_degrees = declination # degrees to radian deg_to_rad = np.pi / 180.0 alpha = alpha_hours * 15.0 * deg_to_rad delta = delta_degrees * deg_to_rad # LIGO detector geometry from # http://www.ligo.org/scientists/GW100916/GW100916-geometry.html # and also # http://www.ligo.org/scientists/GRB051103/GRB051103-geometry.php # The two pages are consistent. # LIGO Hanford (H) detector geometry # Latitude, Longitude and azimuth of the x-arm # Note: Longitude is measured positively westward from Greenwich # Azimuth is measured from the South point, turning positive to the West. lat_H = (46 + 27.0 / 60 + 19.0 / 3600) * deg_to_rad long_H = (119 + 24.0 / 60 + 28.0 / 3600) * deg_to_rad xazi_H = (180 - 36) * deg_to_rad # LIGO Livingston (L) detector geometry # Latitude, Longitude and azimuth of the x-arm lat_L = (30 + 33.0 / 60 + 46.0 / 3600) * deg_to_rad long_L = (90 + 46.0 / 60 + 27.0 / 3600) * deg_to_rad xazi_L = (90 - 18) * deg_to_rad # Virgo detector geometry # https://arxiv.org/pdf/1706.09505.pdf (Table 3.1) lat_V = (43 + 37.0 / 60 + 53.0 / 3600) * deg_to_rad long_V = -(10 + 30.0 / 60 + 16.0 / 3600) * deg_to_rad xazi_V = -(180 - 19) * deg_to_rad # Time of detection of GW150914: 2015-09-14 09:50:45 UTC # The Greenwich sidereal time theta_G is calculated by # the formula at https://en.wikipedia.org/wiki/Sidereal_time base_date = datetime.datetime.strptime( "2000-01-01 12:00:00", "%Y-%m-%d %H:%M:%S" ) date = datetime.datetime.strptime(time_utc, "%Y-%m-%d %H:%M:%S") # Days between DATE and 2000-01-01 12:00 D = (date - base_date).total_seconds() / 86400 theta_G = 18.697374558 + 24.06570982441908 * D # theta_G mod 24h theta_G = theta_G - np.floor(theta_G / 24.0) * 24 theta_G = theta_G * 15 * deg_to_rad # Equatorial -> horizontal using the formulae at # https://en.wikipedia.org/wiki/Celestial_coordinate_system # Hour angles h_H = theta_G - long_H - alpha h_L = theta_G - long_L - alpha h_V = theta_G - long_V - alpha # Zenith distance theta and azimuth A theta_H = np.arccos( np.sin(lat_H) * np.sin(delta) + np.cos(lat_H) * np.cos(delta) * np.cos(h_H) ) theta_L = np.arccos( np.sin(lat_L) * np.sin(delta) + np.cos(lat_L) * np.cos(delta) * np.cos(h_L) ) theta_V = np.arccos( np.sin(lat_V) * np.sin(delta) + np.cos(lat_V) * np.cos(delta) * np.cos(h_V) ) A_H = np.arctan2( np.cos(delta) * np.sin(h_H), np.cos(delta) * np.cos(h_H) * np.sin(lat_H) - np.sin(delta) * np.cos(lat_H), ) phi_H = xazi_H - A_H A_L = np.arctan2( np.cos(delta) * np.sin(h_L), np.cos(delta) * np.cos(h_L) * np.sin(lat_L) - np.sin(delta) * np.cos(lat_L), ) phi_L = xazi_L - A_L A_V = np.arctan2( np.cos(delta) * np.sin(h_V), np.cos(delta) * np.cos(h_V) * np.sin(lat_V) - np.sin(delta) * np.cos(lat_V), ) phi_V = xazi_V - A_V coords = Detectors( hanford=(theta_H, phi_H), livingston=(theta_L, phi_L), virgo=(theta_V, phi_V), ) return coords
[docs]def antenna_responses(theta, phi, polarization=0): """Return the antenna response pattern of a detector on the z = 0 plane with the arms on the x and y directions for a given localization defined by the spherical angles ``theta`` and ``phi``. :param theta: Meridional angle. :type theta: float :param phi: Azimuthal angle. :type phi: float :param polarization: Polarization angle of the wave. :type polarization: float :returns: Antenna response for cross and plus polarizations (in this order). :rtype: tuple of floats. """ # http://research.physics.illinois.edu/cta/movies/bhbh_sim/wavestrain.html Fp = 0.5 * (1 + np.cos(theta) * np.cos(theta)) * np.cos(2 * phi) * np.cos( 2 * polarization ) - np.cos(theta) * np.sin(2 * phi) * np.sin(2 * polarization) Fc = 0.5 * (1 + np.cos(theta) * np.cos(theta)) * np.cos(2 * phi) * np.sin( 2 * polarization ) + np.cos(theta) * np.sin(2 * phi) * np.cos(2 * polarization) return (Fc, Fp)
[docs]def antenna_responses_from_sky_localization( right_ascension, declination, time_utc, polarization=0 ): """Return the antenna responses for Hanford, Livingston and Virgo for a given source. See, http://research.physics.illinois.edu/cta/movies/bhbh_sim/wavestrain.html. ``utc_time`` has to have the following formatting: ``%Y-%m-%d %H:%M``, (eg ``2015-09-14 09:50:45``) :param right_ascension: Right ascension of the source in degrees :type right_ascension: float :param declination: Declination of the source in degrees :type declination: float :param time_utc: UTC time of the event :type declination: str :param polarization: Polarization of the wave :type polarization: float :returns antenna_pattern: Cross and plus antenna pattern for the different interferometers. :rtype: namedtuple with fields hanford, livingston, and virgo """ coords = ra_dec_to_theta_phi(right_ascension, declination, time_utc) Fc_H, Fp_H = antenna_responses(*coords.hanford, polarization) Fc_L, Fp_L = antenna_responses(*coords.livingston, polarization) Fc_V, Fp_V = antenna_responses(*coords.virgo, polarization) antenna = Detectors( hanford=(Fc_H, Fp_H), livingston=(Fc_L, Fp_L), virgo=(Fc_V, Fp_V), ) return antenna
# Extrapolation of GWs at infinity. # Follows what done in the NRAR collaboration (1307.5307) # This is what we are going to do: # 1. Assume that the spacetime is almost Schwarzschild with mass M. # 2. Choose a set of retarded times u_i where GWs have to be evaluated # 3. Compute the coordinate times t_i that correspond to the retarded # times u_i at the radius r. This uses tortoise coordinates. # 4. Interpolate the waveforms at the coordinate times corresponding to # the retarded times u_i for different extraction radii. # Now our waves are so that they are evaluate at different t_i but at # the same u_i. # 5. We do this process for a bunch of extraction radii (rex1, rex2, rex3, ...) # So, we should have wave1, wave2, wave3, with all the same number of # points that correspond to the same retarded time. # 6. For each element in u_i, fit all the waves (wave1, wave2, ...) with a # polynomial of the form a_n/r^n.
[docs]def Schwarzschild_radius_to_tortoise(radii, mass): """Transform radial coordinates ``radii`` to tortoise coordinates assuming mass ``mass``. Equation (26) in 1307.5307. :param radii: Radius in Schwarzschild coordinates. :type radii: float or 1D NumPy array :param mass: ADM mass. :type mass: float :returns: Tortoise radii. :rtype: float or 1D NumPy array """ return radii + 2 * mass * np.log(radii / (2 * mass) - 1)
[docs]def retarded_times_to_coordinate_times(retarded_times, radii, mass): """Compute the coordinate times corresponding to the retarded times at the coordinate radii. First, the tortoise radius is computed from ``radii``, then the coordinate times are computed with :math:`t = u + r_{tortoise}(radii)`. This function is used to extrapolate gravitational waves to infinity. :param retarded_times: Retarded times. :type retarded_times: float or 1D NumPy array :param radii: Radii of evaluation. :type radii: float or 1D NumPy array :param mass: ADM mass (needed to compute tortoise radius). :type mass: float :returns: Coordinate times corresponding to the given retarded times and evaluation radii. :rtype: float or 1D NumPy array """ return retarded_times + Schwarzschild_radius_to_tortoise(radii, mass)
def _coordinate_times_to_retarded_times(coordinate_times, radii, mass): """Compute the coordinate times corresponding to the retarded times at the given coordinate radii. This function is used to extrapolate gravitational waves to infinity. :param retarded_times: Coordinate times. :type retarded_times: float or 1D NumPy array :param radii: Radii (it can be just one) :type radii: float or 1D NumPy array :param mass: ADM mass :type mass: float :returns: Retarded times corresponding to the given coordinate times and evaluation radii. :rtype: float or 1D NumPy array """ return coordinate_times - Schwarzschild_radius_to_tortoise(radii, mass)
[docs]def signal_to_noise_ratio_from_strain( h, *args, noise=None, fmin=0, fmax=np.inf, window_function=None, **kwargs ): r"""Return the signal to noise ratio given a strain and a power spectal density distribution for a detector. If ``window_function`` is not None, the window will be applied to the signal. All the unknown arguments are passed to the window function. The SNR is computed as :math:`sqrt of 4 \int_fmin^fmax |\tilde{h} f|^2 / Sn(f) d f` (equation from 1408.0740) :param h: Gravitational-wave strain. :type h: :py:class:`~.TimeSeries` :param noise: Power spectral density of the noise of the detector. :type noise: :py:class:`~.FrequencySeries` :param fmin: Minimum frequency over which to compute the SNR. :type fmin: float :param fmax: Maximum frequency over which to compute the SNR. :type fmax: float :param window_function: If not None, apply ``window_function`` to the series before computing the strain. :type window_function: callable, str, or None :param args, kwargs: All the additional parameters are passed to the window function. :returns: Signal-to-noise ratio. :rtype: float """ if not isinstance(h, ts.TimeSeries): raise TypeError("Strain has to be a TimeSeries") # First, we window to avoid spectral leakage h_win = h.windowed(window_function, *args, **kwargs) # Then, we take the Fourier transform h_fft = h_win.to_FrequencySeries() # The SNR is just the inner product of h_fft with itself return np.sqrt( h_fft.inner_product(h_fft, noises=noise, fmin=fmin, fmax=fmax) )
[docs]def effective_amplitude_spectral_density( strain, *args, window_function=None, **kwargs ): r"""Return the effective amplitude spectral density for a given strain. If ``window_function`` is not None, the window will be applied to the signal. All the unknown arguments are passed to the window function. The effective amplitude spectral density is computed for a strain (h) as :math:`h_{\rm eff}(f) = f * \sqrt{(|h_{+}|^2 + |h_{\times}|^2) / 2.0}` (see for instance 1604.00246) :param strain: Gravitational-wave strain. :type strain: :py:class:`~.TimeSeries` :param window_function: If not None, apply ``window_function`` to the series before computing the strain. :type window_function: callable, str, or None :param args, kwargs: All the additional parameters are passed to the window function. :returns: Effective amplitude spectral density :rtype: :py:class:`~.FrequencySeries` """ if not isinstance(strain, ts.TimeSeries): raise TypeError("Strain has to be a TimeSeries") # First, ensure we operate on a regularly sampled copy strain_regular = strain.regular_resampled() # Next, we window to avoid spectral leakage if window_function is not None: strain_regular.window(window_function, *args, **kwargs) # We extract the plus and cross components of the strain # Note: for an unaltered kuibit strain, this will in reality be # r * hp, r * hc respectively h_plus = strain_regular.real() h_cross = -strain_regular.imag() # Then, we take the Fourier transform. By extracting hp and hc seperately we # obtain the real signal only in the positive frequency space hp_fft = h_plus.to_FrequencySeries() hc_fft = h_cross.to_FrequencySeries() # Finally, we can compute the effective strain amplitude power spectral # density akin to [Eq (8-9) in 1604.00246]. To save some space, we write # the return value is `hp_fft`. h_eff = hp_fft h_eff.fft = h_eff.f * np.sqrt((hp_fft.amp**2 + hc_fft.amp**2) / 2.0) return h_eff