Align Einstein Telescope geometry with LAL

bilby/gw/detector/detectors/ET.interferometer

@@ -1,6 +1,7 @@
 # The proposed Einstein Telescope located at the site of Virgo.
 # LIGO-T980044-10
 # https://dcc.ligo.org/LIGO-P1600143/public
+# These are the parameters for E1, the other detectors are derived according to T1400308.
 name = 'ET'
 power_spectral_density = PowerSpectralDensity(psd_file='ET_D_psd.txt')
 length = 10

bilby/gw/detector/geometry.py

@@ -2,6 +2,7 @@
 from ..geometry import calculate_arm, detector_tensor
 
 from .. import utils as gwutils
+from ...compat.utils import array_module
 
 
 class InterferometerGeometry(object):
@@ -207,9 +208,12 @@ def vertex(self):
         =======
         array_like: A 3D array representation of the vertex
         """
+
+        xp = self.array_backend
+
         if not self._vertex_updated:
-            self._vertex = gwutils.get_vertex_position_geocentric(self._latitude, self._longitude,
-                                                                  self.elevation)
+            self._vertex = gwutils.get_vertex_position_geocentric(xp.asarray([self._latitude, self._longitude,
+                                                                  self.elevation]))
             self._vertex_updated = True
         return self._vertex
 
@@ -305,7 +309,12 @@ def unit_vector_along_arm(self, arm):
         else:
             raise ValueError("Arm must either be 'x' or 'y'.")
 
-    def set_array_backend(self, xp):
+    @property
+    def array_backend(self):
+        return array_module(self.length)
+
+    @array_backend.setter
+    def array_backend(self, xp):
         self.length = xp.asarray(self.length)
         self.latitude = xp.asarray(self.latitude)
         self.longitude = xp.asarray(self.longitude)

bilby/gw/detector/interferometer.py

@@ -586,9 +586,11 @@ def vertex_position_geocentric(self):
         =======
         array_like: A 3D array representation of the vertex
         """
-        return gwutils.get_vertex_position_geocentric(self.geometry.latitude_radians,
-                                                      self.geometry.longitude_radians,
-                                                      self.geometry.elevation)
+        xp = self.array_backend
+
+        return gwutils.get_vertex_position_geocentric(xp.asarray([self.geometry.latitude_radians,
+                                                                  self.geometry.longitude_radians,
+                                                                  self.geometry.elevation]))
 
     def optimal_snr_squared(self, signal):
         """
@@ -939,10 +941,11 @@ def from_pickle(cls, filename=None):
             raise TypeError('The loaded object is not an Interferometer')
         return res
 
-    def set_array_backend(self, xp):
-        self.geometry.set_array_backend(xp=xp)
-        self.power_spectral_density.set_array_backend(xp=xp)
-
     @property
     def array_backend(self):
-        return array_module(self.geometry.length)
+        return self.geometry.array_backend
+
+    @array_backend.setter
+    def array_backend(self, xp):
+        self.geometry.array_backend = xp
+        self.power_spectral_density.array_backend = xp

bilby/gw/detector/networks.py

@@ -1,13 +1,15 @@
 import os
 
-import numpy as np
 import math
+from scipy.spatial.transform import Rotation
 
+from ...compat.utils import xp_wrap
 from ...core import utils
 from ...core.utils import logger, safe_file_dump
 from ..geometry import zenith_azimuth_to_theta_phi
 from .interferometer import Interferometer
 from .psd import PowerSpectralDensity
+from ..utils import get_vertex_position_geocentric, get_vertex_position_ellipsoid
 
 
 class InterferometerList(list):
@@ -343,14 +345,15 @@ def from_pickle(cls, filename=None):
     )
     from_pickle.__doc__ = _load_docstring.format(format="pickle")
 
-    def set_array_backend(self, xp):
-        for ifo in self:
-            ifo.set_array_backend(xp)
-
     @property
     def array_backend(self):
         return self[0].array_backend
 
+    @array_backend.setter
+    def array_backend(self, xp):
+        for ifo in self:
+            ifo.array_backend = xp
+
 
 class TriangularInterferometer(InterferometerList):
     def __init__(
@@ -378,9 +381,8 @@ def __init__(
         if isinstance(maximum_frequency, float) or isinstance(maximum_frequency, int):
             maximum_frequency = [maximum_frequency] * 3
 
-        brng = 90 - xarm_azimuth
-
         for ii in range(3):
+
             self.append(
                 Interferometer(
                     "{}{}".format(name, ii + 1),
@@ -398,23 +400,84 @@ def __init__(
                 )
             )
 
-            phi1 = np.radians(latitude)
-            phi2 = np.arcsin(
-                np.sin(phi1) * np.cos(length * 1e3 / utils.radius_of_earth) +
-                np.cos(phi1) * np.sin(length * 1e3 / utils.radius_of_earth) * np.cos(np.radians(brng))
-            )
-            latitude = np.degrees(phi2)
+            latitude, longitude, elevation, xarm_azimuth, \
+                yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
 
-            lam1 = np.radians(longitude)
-            lam2 = lam1 + np.arctan2(
-                np.sin(np.radians(brng)) * np.sin(length * 1e3 / utils.radius_of_earth) * np.cos(phi1),
-                np.cos(length * 1e3 / utils.radius_of_earth) - np.sin(phi1) * np.sin(phi2)
-            )
-            longitude = np.degrees(lam2)
+    @xp_wrap
+    def _get_next_vertex_parameters(self, current_index, *, xp):
+        """
+        Get the parameters for the next vertex (counterclockwise) in the triangular interferometer.
+        The location of the next vertex is calculated by moving along the x-arm of the current
+        interferometer in geocentric coordinates. The new latitude, longitude, and elevation are
+        calculated through a coordinate transformation from geocentric to ellipsoidal coordinates.
+        The new unit vectors for the x and y arms are calculated by rotating the current unit vectors
+        by 120 degrees around the normal vector. The new azimuths and tilts are calculated by projecting
+        the new unit vectors onto the local normal, north, and east vectors.
+
+        Parameters
+        ==========
+        current_index: int
+            The index of the current interferometer in the list.
 
-            brng += 240
-            xarm_azimuth += 240
-            yarm_azimuth += 240
+        Returns
+        =======
+        latitude: float
+            The latitude of the next vertex in degrees.
+        longitude: float
+            The longitude of the next vertex in degrees.
+        elevation: float
+            The elevation of the next vertex in meters.
+        xarm_azimuth: float
+            The azimuth of the x-arm of the next interferometer in degrees.
+        yarm_azimuth: float
+            The azimuth of the y-arm of the next interferometer in degrees.
+        xarm_tilt: float
+            The tilt of the x-arm of the next interferometer in radians.
+        yarm_tilt: float
+            The tilt of the y-arm of the next interferometer in radians.
+
+        """
+        current_ifo = self[current_index]
+        unit_vector_x = current_ifo.geometry.unit_vector_along_arm("x")
+        unit_vector_y = current_ifo.geometry.unit_vector_along_arm("y")
+
+        vertex_geocentric = get_vertex_position_geocentric(xp.asarray([current_ifo.latitude_radians,
+                                                                       current_ifo.longitude_radians,
+                                                                       current_ifo.elevation]))
+        next_vertex_geocentric = vertex_geocentric + current_ifo.length * 1000 * unit_vector_x
+        next_vertex_ellipsoid = get_vertex_position_ellipsoid(next_vertex_geocentric)
+        next_latitude_rad, next_longitude_rad, next_elevation = next_vertex_ellipsoid
+
+        rotation_vector = xp.cross(unit_vector_x, unit_vector_y)
+        rotation_vector /= xp.linalg.norm(rotation_vector)
+        rotation_angle = 2 / 3 * xp.pi
+        rotation = Rotation.from_rotvec(rotation_angle * rotation_vector)
+        next_unit_vector_x = rotation.apply(unit_vector_x)
+        next_unit_vector_y = rotation.apply(unit_vector_y)
+
+        next_local_normal_vector = xp.asarray([xp.cos(next_latitude_rad) * xp.cos(next_longitude_rad),
+                                               xp.cos(next_latitude_rad) * xp.sin(next_longitude_rad),
+                                               xp.sin(next_latitude_rad)])
+        next_local_north_vector = xp.asarray([-xp.sin(next_latitude_rad) * xp.cos(next_longitude_rad),
+                                              -xp.sin(next_latitude_rad) * xp.sin(next_longitude_rad),
+                                              xp.cos(next_latitude_rad)])
+        next_local_east_vector = xp.asarray([-xp.sin(next_longitude_rad),
+                                             xp.cos(next_longitude_rad), 0])
+
+        next_xarm_tilt_rad = xp.arcsin(xp.dot(next_unit_vector_x, next_local_normal_vector))
+        next_yarm_tilt_rad = xp.arcsin(xp.dot(next_unit_vector_y, next_local_normal_vector))
+        next_xarm_azimuth_rad = xp.arctan2(xp.dot(next_unit_vector_x, next_local_north_vector),
+                                           xp.dot(next_unit_vector_x, next_local_east_vector))
+        next_yarm_azimuth_rad = xp.arctan2(xp.dot(next_unit_vector_y, next_local_north_vector),
+                                           xp.dot(next_unit_vector_y, next_local_east_vector))
+
+        return (xp.rad2deg(next_latitude_rad),
+                xp.rad2deg(next_longitude_rad),
+                next_elevation,
+                xp.rad2deg(next_xarm_azimuth_rad),
+                xp.rad2deg(next_yarm_azimuth_rad),
+                next_xarm_tilt_rad,
+                next_yarm_tilt_rad)
 
 
 def get_empty_interferometer(name):

bilby/gw/detector/psd.py

@@ -4,6 +4,7 @@
 import array_api_extra as xpx
 import numpy as np
 
+from ...compat.utils import array_module
 from ...core.utils.calculus import interp1d
 from ...core import utils
 from ...core.utils import logger
@@ -376,7 +377,12 @@ def get_noise_realisation(self, number_of_samples, duration, *, random_state=Non
         frequency_domain_strain = xpx.at(frequency_domain_strain, out_of_bounds).set(0j)
         return xp.nan_to_num(frequency_domain_strain), xp.asarray(frequencies)
 
-    def set_array_backend(self, xp):
+    @property
+    def array_backend(self):
+        return array_module(self.frequency_array)
+
+    @array_backend.setter
+    def array_backend(self, xp):
         """ Set the array backend for the cached arrays
 
         Parameters

bilby/gw/likelihood/base.py

@@ -158,7 +158,7 @@ def __init__(
         self.waveform_generator = waveform_generator
         super(GravitationalWaveTransient, self).__init__()
         self.interferometers = InterferometerList(interferometers)
-        self.interferometers.set_array_backend(interferometers.array_backend)
+        self.interferometers.array_backend = interferometers.array_backend
         self.time_marginalization = time_marginalization
         self.distance_marginalization = distance_marginalization
         self.phase_marginalization = phase_marginalization
@@ -171,7 +171,7 @@ def __init__(
         if "geocent" not in time_reference:
             self.time_reference = time_reference
             self.reference_ifo = get_empty_interferometer(self.time_reference)
-            self.reference_ifo.set_array_backend(self.interferometers.array_backend)
+            self.reference_ifo.array_backend = self.interferometers.array_backend
             if self.time_marginalization:
                 logger.info("Cannot marginalise over non-geocenter time.")
                 self.time_marginalization = False
@@ -1086,7 +1086,7 @@ def reference_frame(self, frame):
         else:
             raise ValueError("Unable to parse reference frame {}".format(frame))
         if isinstance(self._reference_frame, InterferometerList):
-            self._reference_frame.set_array_backend(self.interferometers.array_backend)
+            self._reference_frame.array_backend = self.interferometers.array_backend
 
     def get_sky_frame_parameters(self, parameters):
         """

bilby/gw/utils.py

@@ -8,7 +8,7 @@
 
 from .geometry import zenith_azimuth_to_theta_phi
 from .time import greenwich_mean_sidereal_time
-from ..compat.utils import array_module, xp_wrap
+from ..compat.utils import xp_wrap
 from ..core.utils import (logger, run_commandline,
                           check_directory_exists_and_if_not_mkdir,
                           SamplesSummary, theta_phi_to_ra_dec)
@@ -56,7 +56,8 @@ def psd_from_freq_series(freq_data, df, *, xp=None):
     return asd_from_freq_series(freq_data, df, xp=xp) ** 2
 
 
-def get_vertex_position_geocentric(latitude, longitude, elevation):
+@xp_wrap
+def get_vertex_position_geocentric(position, _longitude=None, _elevation=None, *, xp=None):
     """
     Calculate the position of the IFO vertex in geocentric coordinates in meters.
 
@@ -65,19 +66,26 @@ def get_vertex_position_geocentric(latitude, longitude, elevation):
 
     Parameters
     ==========
-    latitude: float
-        Latitude in radians
-    longitude:
-        Longitude in radians
-    elevation:
-        Elevation in meters
+    position: array_like
+        3D representation of the ellipsoidal vertex position (latitude [rad], longitude [rad], elevation [m])
 
     Returns
     =======
     array_like: A 3D representation of the geocentric vertex position
 
     """
-    xp = array_module(latitude)
+
+    if _longitude is not None and _elevation is not None:
+        logger.warning("Syntax of get_vertex_position_geocentric changed in Bilby 3.0.0."
+                       "Please pass latitude, longitude and elevation as the single array-like argument.")
+        latitude = position
+        longitude = _longitude
+        elevation = _elevation
+    elif len(position) == 3:
+        latitude, longitude, elevation = position
+    else:
+        raise ValueError("Position must be an array-like of length 3 containing latitude, longitude, and elevation.")
+
     semi_major_axis = 6378137  # for ellipsoid model of Earth, in m
     semi_minor_axis = 6356752.314  # in m
     radius = semi_major_axis**2 * (semi_major_axis**2 * xp.cos(latitude)**2 +
@@ -88,6 +96,46 @@ def get_vertex_position_geocentric(latitude, longitude, elevation):
     return xp.asarray([x_comp, y_comp, z_comp])
 
 
+@xp_wrap
+def get_vertex_position_ellipsoid(position, *, xp=None):
+    """
+    Calculate the position of the IFO vertex in ellipsoidal coordinates.
+
+    Implementation of Borkowski's algorithm for the exact inverse of get_vertex_position_geocentric().
+    See http://www.astro.uni.torun.pl/~kb/Papers/ASS/Geod-ASS.htm for details.
+    For the correct version of formula 11a, see: http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-GK.htm#GEOD
+    See https://ieeexplore.ieee.org/document/303772 for an overview of inversion algorithms.
+
+    Parameters
+    ==========
+    position: array_like
+        3D representation of the geocentric vertex position in meters (x, y, z)
+
+    Returns
+    =======
+    array_like: A 3D representation of the ellipsoidal vertex position (latitude [rad], longitude [rad], elevation [m])
+    """
+    semi_major_axis = 6378137  # for ellipsoid model of Earth, in m
+    semi_minor_axis = 6356752.314  # in m
+
+    r = xp.sqrt(position[0] ** 2 + position[1] ** 2)
+    E = (semi_minor_axis * position[2] - (semi_major_axis ** 2 - semi_minor_axis ** 2)) / (r * semi_major_axis)
+    F = (semi_minor_axis * position[2] + (semi_major_axis ** 2 - semi_minor_axis ** 2)) / (r * semi_major_axis)
+    P = 4 / 3 * (E * F + 1)
+    Q = 2 * (E ** 2 - F ** 2)
+    D = P ** 3 + Q ** 2
+    v = (xp.sqrt(D) - Q) ** (1 / 3) - (xp.sqrt(D) + Q) ** (1 / 3)
+
+    # Calculate solution in first quadrant and then adjust based on the sign of the original z_comp
+    G = 1 / 2 * (xp.sqrt(E ** 2 + v) + E)
+    t = xp.sqrt(G ** 2 + (F - v * G) / (2 * G - E)) - G
+    latitude = xp.sign(position[2]) * xp.arctan((semi_major_axis * (1 - t ** 2)) / (2 * semi_minor_axis * t))
+    longitude = xp.arctan2(position[1], position[0])
+    elevation = (r - semi_major_axis * t) * xp.cos(latitude) + (position[2] - semi_minor_axis) * xp.sin(latitude)
+
+    return xp.asarray([latitude, longitude, elevation])
+
+
 def inner_product(aa, bb, frequency, PSD):
     """
     Calculate the inner product defined in the matched filter statistic

docs/array_api.rst

@@ -162,7 +162,7 @@ Gravitational-Wave Likelihoods
 The :code:`Bilby` implementation of gravitational-wave likelihood is compatible with the Array API,
 however this requires access to waveform models that support the provided array backend.
 The desired array backend must be explicitly specified for the data,
-using :code:`bilby.gw.detector.networks.InterferometerList.set_array_backend`.
+using :code:`bilby.gw.detector.networks.InterferometerList.array_backend = <backend>`.
 Below is an example using the :code:`ripplegw` package for waveform generation.
 Here, an injection is performed using the standard :code:`LALSimulation` waveform generator,
 and the analysis is then performed using the JIT-compiled likelihood.
@@ -193,7 +193,7 @@ and the analysis is then performed using the JIT-compiled likelihood.
     ifos.inject_signal(parameters=injection_parameters, waveform_generator=injection_wfg)
 
     # set the array backend after the injection
-    ifos.set_array_backend(jnp)
+    ifos.array_backend = jnp
 
     ripple_wfg = bilby.gw.waveform_generator.WaveformGenerator(
         duration=4,