diff --git a/bilby/gw/detector/geometry.py b/bilby/gw/detector/geometry.py
index 5ed8ec2e789c86396e32f4a12ae1f1f19412e41c..d7e1433decbccaef03990c8dc9f1fa99a2350c2e 100644
--- a/bilby/gw/detector/geometry.py
+++ b/bilby/gw/detector/geometry.py
@@ -1,4 +1,5 @@
import numpy as np
+from bilby_cython.geometry import calculate_arm, detector_tensor
from .. import utils as gwutils
@@ -263,7 +264,7 @@ class InterferometerGeometry(object):
if not self._x_updated or not self._y_updated:
_, _ = self.x, self.y # noqa
if not self._detector_tensor_updated:
- self._detector_tensor = 0.5 * (np.einsum('i,j->ij', self.x, self.x) - np.einsum('i,j->ij', self.y, self.y))
+ self._detector_tensor = detector_tensor(x=self.x, y=self.y)
self._detector_tensor_updated = True
return self._detector_tensor
@@ -288,19 +289,18 @@ class InterferometerGeometry(object):
"""
if arm == 'x':
- return self._calculate_arm(self._xarm_tilt, self._xarm_azimuth)
+ return calculate_arm(
+ arm_tilt=self._xarm_tilt,
+ arm_azimuth=self._xarm_azimuth,
+ longitude=self._longitude,
+ latitude=self._latitude
+ )
elif arm == 'y':
- return self._calculate_arm(self._yarm_tilt, self._yarm_azimuth)
+ return calculate_arm(
+ arm_tilt=self._yarm_tilt,
+ arm_azimuth=self._yarm_azimuth,
+ longitude=self._longitude,
+ latitude=self._latitude
+ )
else:
raise ValueError("Arm must either be 'x' or 'y'.")
-
- def _calculate_arm(self, arm_tilt, arm_azimuth):
- e_long = np.array([-np.sin(self._longitude), np.cos(self._longitude), 0])
- e_lat = np.array([-np.sin(self._latitude) * np.cos(self._longitude),
- -np.sin(self._latitude) * np.sin(self._longitude), np.cos(self._latitude)])
- e_h = np.array([np.cos(self._latitude) * np.cos(self._longitude),
- np.cos(self._latitude) * np.sin(self._longitude), np.sin(self._latitude)])
-
- return (np.cos(arm_tilt) * np.cos(arm_azimuth) * e_long +
- np.cos(arm_tilt) * np.sin(arm_azimuth) * e_lat +
- np.sin(arm_tilt) * e_h)
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index 06b298386e3d2d8a8887f89af3d65f3e58919efd..08290f2e578a97e7afcfbd49b32120425f52b48c 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -1,6 +1,11 @@
import os
import numpy as np
+from bilby_cython.geometry import (
+ get_polarization_tensor,
+ three_by_three_matrix_contraction,
+ time_delay_from_geocenter,
+)
from ...core import utils
from ...core.utils import docstring, logger, PropertyAccessor
@@ -268,11 +273,11 @@ class Interferometer(object):
Returns
=======
- array_like: A 3x3 array representation of the antenna response for the specified mode
+ float: The antenna response for the specified mode and time/location
"""
- polarization_tensor = gwutils.get_polarization_tensor(ra, dec, time, psi, mode)
- return np.einsum('ij,ij->', self.geometry.detector_tensor, polarization_tensor)
+ polarization_tensor = get_polarization_tensor(ra, dec, time, psi, mode)
+ return three_by_three_matrix_contraction(self.geometry.detector_tensor, polarization_tensor)
def get_detector_response(self, waveform_polarizations, parameters):
""" Get the detector response for a particular waveform
@@ -527,7 +532,7 @@ class Interferometer(object):
=======
float: The time delay from geocenter in seconds
"""
- return gwutils.time_delay_geocentric(self.geometry.vertex, np.array([0, 0, 0]), ra, dec, time)
+ return time_delay_from_geocenter(self.geometry.vertex, ra, dec, time)
def vertex_position_geocentric(self):
"""
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 34211e82b9f31e0a1c9bcb8e2d1ff28c39b269bb..da9956d2c9b9e1d8e24163e1d183077bd61c8d48 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -1,14 +1,16 @@
import json
import os
from functools import lru_cache
-from math import fmod
import numpy as np
from scipy.interpolate import interp1d
from scipy.special import i0e
+from bilby_cython.geometry import (
+ zenith_azimuth_to_theta_phi as _zenith_azimuth_to_theta_phi,
+)
+from bilby_cython.time import greenwich_mean_sidereal_time
-from ..core.utils import (ra_dec_to_theta_phi,
- speed_of_light, logger, run_commandline,
+from ..core.utils import (logger, run_commandline,
check_directory_exists_and_if_not_mkdir,
SamplesSummary, theta_phi_to_ra_dec)
from ..core.utils.constants import solar_mass
@@ -53,90 +55,6 @@ def psd_from_freq_series(freq_data, df):
return np.power(asd_from_freq_series(freq_data, df), 2)
-def time_delay_geocentric(detector1, detector2, ra, dec, time):
- """
- Calculate time delay between two detectors in geocentric coordinates based on XLALArrivaTimeDiff in TimeDelay.c
-
- Parameters
- ==========
- detector1: array_like
- Cartesian coordinate vector for the first detector in the geocentric frame
- generated by the Interferometer class as self.vertex.
- detector2: array_like
- Cartesian coordinate vector for the second detector in the geocentric frame.
- To get time delay from Earth center, use detector2 = np.array([0,0,0])
- ra: float
- Right ascension of the source in radians
- dec: float
- Declination of the source in radians
- time: float
- GPS time in the geocentric frame
-
- Returns
- =======
- float: Time delay between the two detectors in the geocentric frame
-
- """
- gmst = fmod(greenwich_mean_sidereal_time(time), 2 * np.pi)
- theta, phi = ra_dec_to_theta_phi(ra, dec, gmst)
- omega = np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])
- delta_d = detector2 - detector1
- return np.dot(omega, delta_d) / speed_of_light
-
-
-def get_polarization_tensor(ra, dec, time, psi, mode):
- """
- Calculate the polarization tensor for a given sky location and time
-
- See Nishizawa et al. (2009) arXiv:0903.0528 for definitions of the polarisation tensors.
- [u, v, w] represent the Earth-frame
- [m, n, omega] represent the wave-frame
- Note: there is a typo in the definition of the wave-frame in Nishizawa et al.
-
- Parameters
- ==========
- ra: float
- right ascension in radians
- dec: float
- declination in radians
- time: float
- geocentric GPS time
- psi: float
- binary polarisation angle counter-clockwise about the direction of propagation
- mode: str
- polarisation mode
-
- Returns
- =======
- array_like: A 3x3 representation of the polarization_tensor for the specified mode.
-
- """
- gmst = fmod(greenwich_mean_sidereal_time(time), 2 * np.pi)
- theta, phi = ra_dec_to_theta_phi(ra, dec, gmst)
- u = np.array([np.cos(phi) * np.cos(theta), np.cos(theta) * np.sin(phi), -np.sin(theta)])
- v = np.array([-np.sin(phi), np.cos(phi), 0])
- m = -u * np.sin(psi) - v * np.cos(psi)
- n = -u * np.cos(psi) + v * np.sin(psi)
-
- if mode.lower() == 'plus':
- return np.einsum('i,j->ij', m, m) - np.einsum('i,j->ij', n, n)
- elif mode.lower() == 'cross':
- return np.einsum('i,j->ij', m, n) + np.einsum('i,j->ij', n, m)
- elif mode.lower() == 'breathing':
- return np.einsum('i,j->ij', m, m) + np.einsum('i,j->ij', n, n)
-
- # Calculating omega here to avoid calculation when model in [plus, cross, breathing]
- omega = np.cross(m, n)
- if mode.lower() == 'longitudinal':
- return np.einsum('i,j->ij', omega, omega)
- elif mode.lower() == 'x':
- return np.einsum('i,j->ij', m, omega) + np.einsum('i,j->ij', omega, m)
- elif mode.lower() == 'y':
- return np.einsum('i,j->ij', n, omega) + np.einsum('i,j->ij', omega, n)
- else:
- raise ValueError("{} not a polarization mode!".format(mode))
-
-
def get_vertex_position_geocentric(latitude, longitude, elevation):
"""
Calculate the position of the IFO vertex in geocentric coordinates in meters.
@@ -310,56 +228,6 @@ def overlap(signal_a, signal_b, power_spectral_density=None, delta_frequency=Non
return sum(integral).real
-__cached_euler_matrix = None
-__cached_delta_x = None
-
-
-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.
- """
- global __cached_delta_x
- global __cached_euler_matrix
-
- delta_x = delta_x / np.sum(delta_x**2)**0.5
- if np.array_equal(delta_x, __cached_delta_x):
- return __cached_euler_matrix
- else:
- __cached_delta_x = delta_x
- 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)
- __cached_delta_x = delta_x
- __cached_euler_matrix = total_rotation
- return total_rotation
-
-
def zenith_azimuth_to_theta_phi(zenith, azimuth, ifos):
"""
Convert from the 'detector frame' to the Earth frame.
@@ -379,15 +247,7 @@ def zenith_azimuth_to_theta_phi(zenith, azimuth, ifos):
The zenith and azimuthal angles in the earth frame.
"""
delta_x = ifos[0].geometry.vertex - ifos[1].geometry.vertex
- omega_prime = np.array([
- np.sin(zenith) * np.cos(azimuth),
- np.sin(zenith) * np.sin(azimuth),
- np.cos(zenith)])
- rotation_matrix = euler_rotation(delta_x)
- omega = np.dot(rotation_matrix, omega_prime)
- theta = np.arccos(omega[2])
- phi = np.arctan2(omega[1], omega[0]) % (2 * np.pi)
- return theta, phi
+ return _zenith_azimuth_to_theta_phi(zenith, azimuth, delta_x)
def zenith_azimuth_to_ra_dec(zenith, azimuth, geocent_time, ifos):
@@ -1005,27 +865,6 @@ def plot_spline_pos(log_freqs, samples, nfreqs=100, level=0.9, color='k', label=
plt.xlim(freq_points.min() - .5, freq_points.max() + 50)
-def greenwich_mean_sidereal_time(time):
- """
- Compute the greenwich mean sidereal time from a GPS time.
-
- This is just a wrapper around :code:`lal.GreenwichMeanSiderealTime` .
-
- Parameters
- ----------
- time: float
- The GPS time to convert.
-
- Returns
- -------
- float
- The sidereal time.
- """
- from lal import GreenwichMeanSiderealTime
- time = float(time)
- return GreenwichMeanSiderealTime(time)
-
-
def ln_i0(value):
"""
A numerically stable method to evaluate ln(I_0) a modified Bessel function
diff --git a/requirements.txt b/requirements.txt
index 70e457fcaca3e0d459bcf524f47300af8f1d9e93..b69a3c7ce373c9df519b9d75a35b3e7260bf43cc 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
+bilby.cython>=0.3.0
dynesty<1.1
emcee
corner
diff --git a/test/gw/detector/geometry_test.py b/test/gw/detector/geometry_test.py
index 3960c196471dfe4fd6856f0ecebd2aafc401f187..358825b237cda9729cecf635c3179c02d38d7b8d 100644
--- a/test/gw/detector/geometry_test.py
+++ b/test/gw/detector/geometry_test.py
@@ -138,81 +138,56 @@ class TestInterferometerGeometry(unittest.TestCase):
self.geometry.latitude = 0
self.assertTrue(np.array_equal(self.geometry.y, np.array([1])))
- def test_detector_tensor_without_update(self):
- _ = self.geometry.detector_tensor
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- expected = np.array(
- [
- [-9.24529394e-06, 1.02425803e-04, 3.24550668e-04],
- [1.02425803e-04, 1.37390844e-03, -8.61137566e-03],
- [3.24550668e-04, -8.61137566e-03, -1.36466315e-03],
- ]
- )
- self.assertTrue(np.allclose(expected, self.geometry.detector_tensor))
-
def test_detector_tensor_with_x_azimuth_update(self):
- _ = self.geometry.detector_tensor
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- self.geometry.xarm_azimuth = 1
- self.assertEqual(0, self.geometry.detector_tensor)
+ original = self.geometry.detector_tensor
+ self.geometry.xarm_azimuth += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_detector_tensor_with_y_azimuth_update(self):
- _ = self.geometry.detector_tensor
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- self.geometry.yarm_azimuth = 1
- self.assertEqual(0, self.geometry.detector_tensor)
+ original = self.geometry.detector_tensor
+ self.geometry.yarm_azimuth += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_detector_tensor_with_x_tilt_update(self):
- _ = self.geometry.detector_tensor
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- self.geometry.xarm_tilt = 1
- self.assertEqual(0, self.geometry.detector_tensor)
+ original = self.geometry.detector_tensor
+ self.geometry.xarm_tilt += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_detector_tensor_with_y_tilt_update(self):
- _ = self.geometry.detector_tensor
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- self.geometry.yarm_tilt = 1
- self.assertEqual(0, self.geometry.detector_tensor)
+ original = self.geometry.detector_tensor
+ self.geometry.yarm_tilt += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_detector_tensor_with_longitude_update(self):
- with mock.patch("numpy.einsum") as m:
- m.return_value = 1
- self.geometry.longitude = 1
- self.assertEqual(0, self.geometry.detector_tensor)
+ original = self.geometry.detector_tensor
+ self.geometry.longitude += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_detector_tensor_with_latitude_update(self):
- with mock.patch("numpy.einsum") as m:
- _ = self.geometry.detector_tensor
- m.return_value = 1
- self.geometry.latitude = 1
- self.assertEqual(self.geometry.detector_tensor, 0)
+ original = self.geometry.detector_tensor
+ self.geometry.latitude += 1
+ self.assertNotEqual(np.max(abs(self.geometry.detector_tensor - original)), 0)
def test_unit_vector_along_arm_default(self):
with self.assertRaises(ValueError):
self.geometry.unit_vector_along_arm("z")
def test_unit_vector_along_arm_x(self):
- with mock.patch("numpy.array") as m:
- m.return_value = 1
- self.geometry.xarm_tilt = 0
- self.geometry.xarm_azimuth = 0
- self.geometry.yarm_tilt = 0
- self.geometry.yarm_azimuth = 90
- self.assertAlmostEqual(self.geometry.unit_vector_along_arm("x"), 1)
+ self.geometry.longitude = 0
+ self.geometry.latitude = 0
+ self.geometry.xarm_tilt = 0
+ self.geometry.xarm_azimuth = 0
+ arm = self.geometry.unit_vector_along_arm("x")
+ self.assertTrue(np.allclose(arm, np.array([0, 1, 0])))
def test_unit_vector_along_arm_y(self):
- with mock.patch("numpy.array") as m:
- m.return_value = 1
- self.geometry.xarm_tilt = 0
- self.geometry.xarm_azimuth = 90
- self.geometry.yarm_tilt = 0
- self.geometry.yarm_azimuth = 180
- self.assertAlmostEqual(self.geometry.unit_vector_along_arm("y"), -1)
+ self.geometry.longitude = 0
+ self.geometry.latitude = 0
+ self.geometry.yarm_tilt = 0
+ self.geometry.yarm_azimuth = 90
+ arm = self.geometry.unit_vector_along_arm("y")
+ print(arm)
+ self.assertTrue(np.allclose(arm, np.array([0, 0, 1])))
def test_repr(self):
expected = (
diff --git a/test/gw/detector/interferometer_test.py b/test/gw/detector/interferometer_test.py
index ee575317a0c0f075bff6cbfe04e83711c7c2a038..ad324e00726980b555fec731eb61eeaa2418fffe 100644
--- a/test/gw/detector/interferometer_test.py
+++ b/test/gw/detector/interferometer_test.py
@@ -97,21 +97,6 @@ class TestInterferometer(unittest.TestCase):
def test_max_freq_setting(self):
self.assertEqual(self.ifo.strain_data.maximum_frequency, self.maximum_frequency)
- def test_antenna_response_default(self):
- with mock.patch("bilby.gw.utils.get_polarization_tensor") as m:
- with mock.patch("numpy.einsum") as n:
- m.return_value = 0
- n.return_value = 1
- self.assertEqual(self.ifo.antenna_response(234, 52, 54, 76, "plus"), 1)
-
- def test_antenna_response_einsum(self):
- with mock.patch("bilby.gw.utils.get_polarization_tensor") as m:
- m.return_value = np.ones((3, 3))
- self.assertAlmostEqual(
- self.ifo.antenna_response(234, 52, 54, 76, "plus"),
- self.ifo.detector_tensor.sum(),
- )
-
def test_get_detector_response_default_behaviour(self):
self.ifo.antenna_response = mock.MagicMock(return_value=1)
self.ifo.time_delay_from_geocenter = mock.MagicMock(return_value=0)
@@ -315,16 +300,6 @@ class TestInterferometer(unittest.TestCase):
with self.assertRaises(ValueError):
self.ifo.inject_signal(injection_polarizations=None, parameters=None)
- def test_time_delay_from_geocenter(self):
- with mock.patch("bilby.gw.utils.time_delay_geocentric") as m:
- m.return_value = 1
- self.assertEqual(self.ifo.time_delay_from_geocenter(1, 2, 3), 1)
-
- def test_vertex_position_geocentric(self):
- with mock.patch("bilby.gw.utils.get_vertex_position_geocentric") as m:
- m.return_value = 1
- self.assertEqual(self.ifo.vertex_position_geocentric(), 1)
-
def test_optimal_snr_squared(self):
"""
Merely checks parameters are given in the right order and the frequency
diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index c5f3af46bb84558ab5be67f832de5657f027720d..7946a99c83b99958160636f482f4593103f6c7fe 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -36,34 +36,6 @@ class TestGWUtils(unittest.TestCase):
psd = gwutils.psd_from_freq_series(freq_data, df)
self.assertTrue(np.all(psd == (freq_data * 2 * df ** 0.5) ** 2))
- def test_time_delay_from_geocenter(self):
- """
- The difference in the two detector case is due to rounding error.
- Different hardware gives different numbers in the last decimal place.
- """
- det1 = np.array([0.1, 0.2, 0.3])
- det2 = np.array([0.1, 0.2, 0.5])
- ra = 0.5
- dec = 0.2
- time = 10
- self.assertEqual(gwutils.time_delay_geocentric(det1, det1, ra, dec, time), 0)
- self.assertAlmostEqual(
- gwutils.time_delay_geocentric(det1, det2, ra, dec, time),
- 1.3253791114055397e-10,
- 14,
- )
-
- def test_get_polarization_tensor(self):
- ra = 1
- dec = 2.0
- time = 10
- psi = 0.1
- for mode in ["plus", "cross", "breathing", "longitudinal", "x", "y"]:
- p = gwutils.get_polarization_tensor(ra, dec, time, psi, mode)
- self.assertEqual(p.shape, (3, 3))
- with self.assertRaises(ValueError):
- gwutils.get_polarization_tensor(ra, dec, time, psi, "not-a-mode")
-
def test_inner_product(self):
aa = np.array([1, 2, 3])
bb = np.array([5, 6, 7])