diff --git a/bilby/core/utils.py b/bilby/core/utils.py
index d631610915bfb24010ab3d0637aaad3f65466283..4ab90eae13547a477d756d024ec3b0904fe5e1f4 100644
--- a/bilby/core/utils.py
+++ b/bilby/core/utils.py
@@ -311,6 +311,12 @@ def ra_dec_to_theta_phi(ra, dec, gmst):
return theta, phi
+def theta_phi_to_ra_dec(theta, phi, gmst):
+ ra = phi + gmst
+ dec = np.pi / 2 - theta
+ return ra, dec
+
+
def gps_time_to_gmst(gps_time):
"""
Convert gps time to Greenwich mean sidereal time in radians
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index 469eef80e594f6a06013be2945d1262366855fcd..d154577659184ba3fab01842bd6c0b0ee90976fa 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -719,6 +719,20 @@ def _generate_all_cbc_parameters(sample, defaults, base_conversion,
['luminosity_distance', 'phase', 'geocent_time']):
if getattr(likelihood, '{}_marginalization'.format(par), False):
priors[name] = likelihood.priors[name]
+
+ if (
+ not getattr(likelihood, "reference_frame", "sky") == "sky"
+ or not getattr(likelihood, "time_reference", "geocenter") == "geocenter"
+ ):
+ try:
+ generate_sky_frame_parameters(
+ samples=output_sample, likelihood=likelihood
+ )
+ except TypeError:
+ logger.info(
+ "Failed to generate sky frame parameters for type {}"
+ .format(type(output_sample))
+ )
for key, func in zip(["mass", "spin", "source frame"], [
generate_mass_parameters, generate_spin_parameters,
generate_source_frame_parameters]):
@@ -1115,3 +1129,21 @@ def generate_posterior_samples_from_marginalized_likelihood(
samples['luminosity_distance'] = new_distance_samples
samples['phase'] = new_phase_samples
return samples
+
+
+def generate_sky_frame_parameters(samples, likelihood):
+ if isinstance(samples, dict):
+ likelihood.parameters.update(samples)
+ samples.update(likelihood.get_sky_frame_parameters())
+ elif not isinstance(samples, DataFrame):
+ raise ValueError
+
+ logger.info('Generating sky frame parameters.')
+ new_samples = list()
+ for ii in tqdm(range(len(samples)), file=sys.stdout):
+ sample = dict(samples.iloc[ii]).copy()
+ likelihood.parameters.update(sample)
+ new_samples.append(likelihood.get_sky_frame_parameters())
+ new_samples = DataFrame(new_samples)
+ for key in new_samples:
+ samples[key] = new_samples[key]
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index bdfb8de545e9d2050a3e0059ad3520e76e5bf903..e2ca541ae2f9b03553c7571089ac7f67bd0491ba 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -21,10 +21,12 @@ from ..core.utils import (
logger, UnsortedInterp2d, create_frequency_series, create_time_series,
speed_of_light, radius_of_earth)
from ..core.prior import Interped, Prior, Uniform
-from .detector import InterferometerList
+from .detector import InterferometerList, get_empty_interferometer
from .prior import BBHPriorDict, CBCPriorDict, Cosmological
from .source import lal_binary_black_hole
-from .utils import noise_weighted_inner_product, build_roq_weights, blockwise_dot_product
+from .utils import (
+ noise_weighted_inner_product, build_roq_weights, blockwise_dot_product,
+ kappa_eta_to_ra_dec)
from .waveform_generator import WaveformGenerator
from collections import namedtuple
@@ -94,11 +96,12 @@ class GravitationalWaveTransient(Likelihood):
'complex_matched_filter_snr',
'd_inner_h_squared_tc_array'])
- def __init__(self, interferometers, waveform_generator,
- time_marginalization=False, distance_marginalization=False,
- phase_marginalization=False, priors=None,
- distance_marginalization_lookup_table=None,
- jitter_time=True):
+ def __init__(
+ self, interferometers, waveform_generator, time_marginalization=False,
+ distance_marginalization=False, phase_marginalization=False, priors=None,
+ distance_marginalization_lookup_table=None, jitter_time=True,
+ reference_frame="sky", time_reference="geocenter"
+ ):
self.waveform_generator = waveform_generator
super(GravitationalWaveTransient, self).__init__(dict())
@@ -109,6 +112,17 @@ class GravitationalWaveTransient(Likelihood):
self.priors = priors
self._check_set_duration_and_sampling_frequency_of_waveform_generator()
self.jitter_time = jitter_time
+ self.reference_frame = reference_frame
+ if "geocent" not in time_reference:
+ self.time_reference = time_reference
+ self.reference_ifo = get_empty_interferometer(self.time_reference)
+ if self.time_marginalization:
+ logger.info("Cannot marginalise over non-geocenter time.")
+ self.time_marginalization = False
+ self.jitter_time = False
+ else:
+ self.time_reference = "geocent"
+ self.reference_ifo = None
if self.time_marginalization:
self._check_prior_is_set(key='geocent_time')
@@ -261,6 +275,8 @@ class GravitationalWaveTransient(Likelihood):
waveform_polarizations =\
self.waveform_generator.frequency_domain_strain(self.parameters)
+ self.parameters.update(self.get_sky_frame_parameters())
+
if waveform_polarizations is None:
return np.nan_to_num(-np.inf)
@@ -702,6 +718,40 @@ class GravitationalWaveTransient(Likelihood):
for mode in signal:
signal[mode] *= self._ref_dist / new_distance
+ @property
+ def reference_frame(self):
+ return self._reference_frame
+
+ @reference_frame.setter
+ def reference_frame(self, frame):
+ if frame == "sky":
+ self._reference_frame = frame
+ elif isinstance(frame, InterferometerList):
+ self._reference_frame = frame[:2]
+ elif isinstance(frame, list):
+ self._reference_frame = InterferometerList(frame[:2])
+ elif isinstance(frame, str):
+ self._reference_frame = InterferometerList([frame[:2], frame[2:4]])
+
+ def get_sky_frame_parameters(self):
+ time = self.parameters['{}_time'.format(self.time_reference)]
+ if not self.reference_frame == "sky":
+ ra, dec = kappa_eta_to_ra_dec(
+ self.parameters['kappa'], self.parameters['eta'],
+ time, self.reference_frame)
+ else:
+ ra = self.parameters["ra"]
+ dec = self.parameters["dec"]
+ if "geocent" not in self.time_reference:
+ geocent_time = (
+ time - self.reference_ifo.time_delay_from_geocenter(
+ ra=ra, dec=dec, time=time
+ )
+ )
+ else:
+ geocent_time = self.parameters["geocent_time"]
+ return dict(ra=ra, dec=dec, geocent_time=geocent_time)
+
@property
def lal_version(self):
try:
@@ -859,11 +909,15 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
'.distance_marginalization_lookup_dmin{}_dmax{}_n{}.npz'
"""
- def __init__(self, interferometers, waveform_generator, priors,
- weights=None, linear_matrix=None, quadratic_matrix=None,
- roq_params=None, roq_params_check=True, roq_scale_factor=1,
- distance_marginalization=False, phase_marginalization=False,
- distance_marginalization_lookup_table=None):
+ def __init__(
+ self, interferometers, waveform_generator, priors,
+ weights=None, linear_matrix=None, quadratic_matrix=None,
+ roq_params=None, roq_params_check=True, roq_scale_factor=1,
+ distance_marginalization=False, phase_marginalization=False,
+ distance_marginalization_lookup_table=None,
+ reference_frame="sky", time_reference="geocenter"
+
+ ):
super(ROQGravitationalWaveTransient, self).__init__(
interferometers=interferometers,
waveform_generator=waveform_generator, priors=priors,
@@ -871,7 +925,10 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
phase_marginalization=phase_marginalization,
time_marginalization=False,
distance_marginalization_lookup_table=distance_marginalization_lookup_table,
- jitter_time=False)
+ jitter_time=False,
+ reference_frame=reference_frame,
+ time_reference=time_reference
+ )
self.roq_params_check = roq_params_check
self.roq_scale_factor = roq_scale_factor
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index dc5ea7964555a01d97f4bfdfc1815c733a308277..a9c70a2e0daf0f0a24f0d37dca866e60c035cc4b 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -10,7 +10,7 @@ import matplotlib.pyplot as plt
from ..core.utils import (ra_dec_to_theta_phi,
speed_of_light, logger, run_commandline,
check_directory_exists_and_if_not_mkdir,
- SamplesSummary)
+ SamplesSummary, theta_phi_to_ra_dec)
try:
from gwpy.timeseries import TimeSeries
@@ -279,6 +279,98 @@ def optimal_snr_squared(signal, power_spectral_density, duration):
return noise_weighted_inner_product(signal, signal, power_spectral_density, duration)
+def euler_rotation(delta_x):
+ """
+ Calculate the rotation matrix mapping the vector (0, 0, 1) to delta_x
+ while preserving the origin of the azimuthal angle.
+
+ This is decomposed into three Euler angle, alpha, beta, gamma, which rotate
+ about the z-, y-, and z- axes respectively.
+
+ Parameters
+ ----------
+ delta_x: array-like (3,)
+ Vector onto which (0, 0, 1) should be mapped.
+
+ Returns
+ -------
+ total_rotation: array-like (3,3)
+ Rotation matrix which maps vectors from the frame in which delta_x is
+ aligned with the z-axis to the target frame.
+ """
+ delta_x = delta_x / np.sum(delta_x**2)**0.5
+ alpha = np.arctan(- delta_x[1] * delta_x[2] / delta_x[0])
+ beta = np.arccos(delta_x[2])
+ gamma = np.arctan(delta_x[1] / delta_x[0])
+ rotation_1 = np.array([
+ [np.cos(alpha), -np.sin(alpha), 0], [np.sin(alpha), np.cos(alpha), 0],
+ [0, 0, 1]])
+ rotation_2 = np.array([
+ [np.cos(beta), 0, - np.sin(beta)], [0, 1, 0],
+ [np.sin(beta), 0, np.cos(beta)]])
+ rotation_3 = np.array([
+ [np.cos(gamma), -np.sin(gamma), 0], [np.sin(gamma), np.cos(gamma), 0],
+ [0, 0, 1]])
+ total_rotation = np.einsum(
+ 'ij,jk,kl->il', rotation_3, rotation_2, rotation_1)
+ return total_rotation
+
+
+def kappa_eta_to_theta_phi(kappa, eta, ifos):
+ """
+ Convert from the 'detector frame' to the Earth frame.
+
+ Parameters
+ kappa: float
+ The zenith angle in the detector frame
+ eta: float
+ The azimuthal angle in the detector frame
+ ifos: list
+ List of Interferometer objects defining the detector frame
+
+ Returns
+ -------
+ theta, phi: float
+ The zenith and azimuthal angles in the earth frame.
+ """
+ delta_x = (ifos[0].vertex_position_geocentric() -
+ ifos[1].vertex_position_geocentric())
+ omega_prime = np.array([
+ np.sin(kappa) * np.cos(eta), np.sin(kappa) * np.sin(eta),
+ np.cos(kappa)])
+ rotation_matrix = euler_rotation(delta_x)
+ omega = np.einsum('ij,j->i', rotation_matrix, omega_prime)
+ theta = np.arccos(omega[2])
+ phi = np.arctan2(omega[1], omega[0]) % (2 * np.pi)
+ return theta, phi
+
+
+def kappa_eta_to_ra_dec(kappa, eta, geocent_time, ifos):
+ """
+ Convert from the 'detector frame' to the Earth frame.
+
+ Parameters
+ kappa: float
+ The zenith angle in the detector frame
+ eta: float
+ The azimuthal angle in the detector frame
+ geocent_time: float
+ GPS time at geocenter
+ ifos: list
+ List of Interferometer objects defining the detector frame
+
+ Returns
+ -------
+ ra, dec: float
+ The zenith and azimuthal angles in the sky frame.
+ """
+ theta, phi = kappa_eta_to_theta_phi(kappa, eta, ifos)
+ gmst = lal.GreenwichMeanSiderealTime(geocent_time)
+ ra, dec = theta_phi_to_ra_dec(theta, phi, gmst)
+ ra = ra % (2 * np.pi)
+ return ra, dec
+
+
def get_event_time(event):
"""
Get the merger time for known GW events.