From a5f498e208e3b50c645873b635d55e61ec6611b1 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Mon, 13 Jan 2025 13:40:36 +0100
Subject: [PATCH 01/15] naming ET detectors anti-clockwise

---
 bilby/gw/detector/detectors/ET.interferometer | 3 +++
 bilby/gw/detector/networks.py                 | 6 ++++--
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/bilby/gw/detector/detectors/ET.interferometer b/bilby/gw/detector/detectors/ET.interferometer
index e89384095..09952a891 100644
--- a/bilby/gw/detector/detectors/ET.interferometer
+++ b/bilby/gw/detector/detectors/ET.interferometer
@@ -1,6 +1,8 @@
 # The proposed Einstein Telescope located at the site of Virgo.
 # LIGO-T980044-10
 # https://dcc.ligo.org/LIGO-P1600143/public
+# fixme: based on the internal document T1400308
+# These are the parameters for E1, the other sub-detectors are derived from this.
 name = 'ET'
 power_spectral_density = PowerSpectralDensity(psd_file='ET_D_psd.txt')
 length = 10
@@ -12,3 +14,4 @@ elevation = 51.884
 xarm_azimuth = 70.5674
 yarm_azimuth = 130.5674
 shape = 'Triangle'
+clockwise = False
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index 39a6cee33..a0ffdcf06 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -310,6 +310,7 @@ def __init__(
         yarm_azimuth,
         xarm_tilt=0.0,
         yarm_tilt=0.0,
+        clockwise=True,
     ):
         super(TriangularInterferometer, self).__init__([])
         self.name = name
@@ -322,6 +323,7 @@ def __init__(
             maximum_frequency = [maximum_frequency] * 3
 
         for ii in range(3):
+
             self.append(
                 Interferometer(
                     "{}{}".format(name, ii + 1),
@@ -339,8 +341,8 @@ def __init__(
                 )
             )
 
-            xarm_azimuth += 240
-            yarm_azimuth += 240
+            xarm_azimuth += 240 if clockwise else 120
+            yarm_azimuth += 240 if clockwise else 120
 
             latitude += (
                 np.arctan(

From ca46c1e2214d35b195805fed9cfd6d7b3b9b7367 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Thu, 16 Jan 2025 11:32:50 +0100
Subject: [PATCH 02/15] Implemented transform ECEF to WGS84

---
 bilby/gw/utils.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index bbbc6036a..9aabc75b5 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -86,6 +86,49 @@ def get_vertex_position_geocentric(latitude, longitude, elevation):
     return np.array([x_comp, y_comp, z_comp])
 
 
+def get_vertex_position_ellipsoid(x_comp, y_comp, z_comp):
+    """
+    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
+    ==========
+    x_comp: float
+        Geocentric x-coordinate in meters
+    y_comp: float
+        Geocentric y-coordinate in meters
+    z_comp: float
+        Geocentric z-coordinate in meters
+
+    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 = np.sqrt(x_comp ** 2 + y_comp ** 2)
+    E = (semi_minor_axis * z_comp - (semi_major_axis ** 2 - semi_minor_axis ** 2)) / (r * semi_major_axis)
+    F = (semi_minor_axis * z_comp + (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 = (np.sqrt(D) - Q) ** (1 / 3) - (np.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 * (np.sqrt(E ** 2 + v) + E)
+    t = np.sqrt(G ** 2 + (F - v * G) / (2 * G - E)) - G
+    latitude = np.sign(z_comp) * np.arctan((semi_major_axis * (1 - t ** 2)) / (2 * semi_minor_axis * t))
+    longitude = np.arctan2(y_comp, x_comp)
+    elevation = (r - semi_major_axis * t) * np.cos(latitude) + (z_comp - semi_minor_axis) * np.sin(latitude)
+
+    return np.array([latitude, longitude, elevation])
+
+
 def inner_product(aa, bb, frequency, PSD):
     """
     Calculate the inner product defined in the matched filter statistic

From 3fa49f0c58ffd61a9e75db746077f5d012f577e3 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Thu, 16 Jan 2025 15:22:33 +0100
Subject: [PATCH 03/15] added test for ECEF to WGS84 conversion

---
 test/gw/utils_test.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index bb842cc80..384988734 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -243,6 +243,16 @@ def test_lalsim_SimInspiralChooseFDWaveform(self):
                 1.5,
             )
 
+    def test_get_vertex_position_geocentric_ellipsoid_conversion(self):
+        latitude_true = np.deg2rad(46 + 27. / 60 + 18.528 / 3600)
+        longitude_true = np.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
+        elevation_true = 142.554
+        vertex = gwutils.get_vertex_position_geocentric(latitude_true, longitude_true, elevation_true)
+        latitude, longitude, elevation = gwutils.get_vertex_position_ellipsoid(vertex[0], vertex[1], vertex[2])
+        self.assertAlmostEqual(latitude, latitude_true, 5)
+        self.assertAlmostEqual(longitude, longitude_true, 5)
+        self.assertAlmostEqual(elevation, elevation_true, 5)
+
 
 class TestSkyFrameConversion(unittest.TestCase):
 

From 058d588ad97d03f8562eea2edc10514794741109 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Fri, 17 Jan 2025 17:36:41 +0100
Subject: [PATCH 04/15] Changed triangular config in line with LAL definition
 of ET

---
 bilby/gw/detector/detectors/ET.interferometer |  3 +-
 bilby/gw/detector/networks.py                 | 67 ++++++++++++-------
 2 files changed, 45 insertions(+), 25 deletions(-)

diff --git a/bilby/gw/detector/detectors/ET.interferometer b/bilby/gw/detector/detectors/ET.interferometer
index 09952a891..79b81d1d0 100644
--- a/bilby/gw/detector/detectors/ET.interferometer
+++ b/bilby/gw/detector/detectors/ET.interferometer
@@ -1,8 +1,7 @@
 # The proposed Einstein Telescope located at the site of Virgo.
 # LIGO-T980044-10
 # https://dcc.ligo.org/LIGO-P1600143/public
-# fixme: based on the internal document T1400308
-# These are the parameters for E1, the other sub-detectors are derived from this.
+# 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
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index a0ffdcf06..3edb87bd0 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -2,11 +2,13 @@
 
 import numpy as np
 import math
+from scipy.spatial.transform import Rotation
 
 from ...core import utils
 from ...core.utils import logger, safe_file_dump
 from .interferometer import Interferometer
 from .psd import PowerSpectralDensity
+from ..utils import get_vertex_position_geocentric, get_vertex_position_ellipsoid
 
 
 class InterferometerList(list):
@@ -341,29 +343,48 @@ def __init__(
                 )
             )
 
-            xarm_azimuth += 240 if clockwise else 120
-            yarm_azimuth += 240 if clockwise else 120
-
-            latitude += (
-                np.arctan(
-                    length
-                    * np.sin(xarm_azimuth * np.pi / 180)
-                    * 1e3
-                    / utils.radius_of_earth
-                )
-                * 180
-                / np.pi
-            )
-            longitude += (
-                np.arctan(
-                    length
-                    * np.cos(xarm_azimuth * np.pi / 180)
-                    * 1e3
-                    / utils.radius_of_earth
-                )
-                * 180
-                / np.pi
-            )
+            unit_vector_x = self[ii].geometry.unit_vector_along_arm("x")
+            unit_vector_y = self[ii].geometry.unit_vector_along_arm("y")
+
+            vertex_geocentric = get_vertex_position_geocentric(self[ii].latitude_radians,
+                                                               self[ii].longitude_radians,
+                                                               self[ii].elevation)
+            next_vertex_geocentric = vertex_geocentric + length * 1000 * (unit_vector_y if clockwise else unit_vector_x)
+            next_vertex_ellipsoid = get_vertex_position_ellipsoid(next_vertex_geocentric[0],
+                                                                  next_vertex_geocentric[1],
+                                                                  next_vertex_geocentric[2])
+            next_latitude_rad, next_longitude_rad, next_elevation = next_vertex_ellipsoid
+
+            rotation_vector = np.cross(unit_vector_x, unit_vector_y)
+            rotation_vector /= np.linalg.norm(rotation_vector)
+            rotation_angle = - 2 / 3 * np.pi if clockwise else 2 / 3 * np.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 = np.array([np.cos(next_latitude_rad) * np.cos(next_longitude_rad),
+                                                 np.cos(next_latitude_rad) * np.sin(next_longitude_rad),
+                                                 np.sin(next_latitude_rad)])
+            next_local_north_vector = np.array([-np.sin(next_latitude_rad) * np.cos(next_longitude_rad),
+                                                -np.sin(next_latitude_rad) * np.sin(next_longitude_rad),
+                                                np.cos(next_latitude_rad)])
+            next_local_east_vector = np.array([-np.sin(next_longitude_rad),
+                                               np.cos(next_longitude_rad), 0])
+
+            next_xarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_x, next_local_normal_vector))
+            next_yarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_y, next_local_normal_vector))
+            next_xarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_x, next_local_north_vector),
+                                               np.dot(next_unit_vector_x, next_local_east_vector))
+            next_yarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_y, next_local_north_vector),
+                                               np.dot(next_unit_vector_y, next_local_east_vector))
+
+            latitude = np.rad2deg(next_latitude_rad)
+            longitude = np.rad2deg(next_longitude_rad)
+            elevation = next_elevation
+            xarm_azimuth = np.rad2deg(next_xarm_azimuth_rad)
+            yarm_azimuth = np.rad2deg(next_yarm_azimuth_rad)
+            xarm_tilt = next_xarm_tilt_rad
+            yarm_tilt = next_yarm_tilt_rad
 
 
 def get_empty_interferometer(name):

From c12e5769aa9c608052f267bc3c49824cb2779442 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Mon, 14 Sep 2026 15:17:07 +0200
Subject: [PATCH 05/15] detector: removed the clockwise option in
 TriangularInterferometer

---
 bilby/gw/detector/detectors/ET.interferometer | 1 -
 bilby/gw/detector/networks.py                 | 5 ++---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/bilby/gw/detector/detectors/ET.interferometer b/bilby/gw/detector/detectors/ET.interferometer
index 79b81d1d0..be55be9b8 100644
--- a/bilby/gw/detector/detectors/ET.interferometer
+++ b/bilby/gw/detector/detectors/ET.interferometer
@@ -13,4 +13,3 @@ elevation = 51.884
 xarm_azimuth = 70.5674
 yarm_azimuth = 130.5674
 shape = 'Triangle'
-clockwise = False
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index 3edb87bd0..9cb1ca938 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -312,7 +312,6 @@ def __init__(
         yarm_azimuth,
         xarm_tilt=0.0,
         yarm_tilt=0.0,
-        clockwise=True,
     ):
         super(TriangularInterferometer, self).__init__([])
         self.name = name
@@ -349,7 +348,7 @@ def __init__(
             vertex_geocentric = get_vertex_position_geocentric(self[ii].latitude_radians,
                                                                self[ii].longitude_radians,
                                                                self[ii].elevation)
-            next_vertex_geocentric = vertex_geocentric + length * 1000 * (unit_vector_y if clockwise else unit_vector_x)
+            next_vertex_geocentric = vertex_geocentric + length * 1000 * unit_vector_x
             next_vertex_ellipsoid = get_vertex_position_ellipsoid(next_vertex_geocentric[0],
                                                                   next_vertex_geocentric[1],
                                                                   next_vertex_geocentric[2])
@@ -357,7 +356,7 @@ def __init__(
 
             rotation_vector = np.cross(unit_vector_x, unit_vector_y)
             rotation_vector /= np.linalg.norm(rotation_vector)
-            rotation_angle = - 2 / 3 * np.pi if clockwise else 2 / 3 * np.pi
+            rotation_angle = 2 / 3 * np.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)

From 97cc3cd4ec3ba60f465ff12d0dc8e4bb184989c8 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Wed, 16 Sep 2026 13:09:23 +0200
Subject: [PATCH 06/15] Changed function signature of get_vertex_* functions

---
 bilby/gw/detector/geometry.py       |  4 +--
 bilby/gw/detector/interferometer.py |  6 ++--
 bilby/gw/detector/networks.py       | 10 +++----
 bilby/gw/utils.py                   | 43 +++++++++++++++--------------
 test/gw/utils_test.py               |  4 +--
 5 files changed, 34 insertions(+), 33 deletions(-)

diff --git a/bilby/gw/detector/geometry.py b/bilby/gw/detector/geometry.py
index d7e1433de..3ae3ddc22 100644
--- a/bilby/gw/detector/geometry.py
+++ b/bilby/gw/detector/geometry.py
@@ -208,8 +208,8 @@ def vertex(self):
         array_like: A 3D array representation of the vertex
         """
         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(np.array([self._latitude, self._longitude,
+                                                                  self.elevation]))
             self._vertex_updated = True
         return self._vertex
 
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index ac0c14886..87d20f91f 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -558,9 +558,9 @@ 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)
+        return gwutils.get_vertex_position_geocentric(np.array([self.geometry.latitude_radians,
+                                                                self.geometry.longitude_radians,
+                                                                self.geometry.elevation]))
 
     def optimal_snr_squared(self, signal):
         """
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index 9cb1ca938..2612b851c 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -345,13 +345,11 @@ def __init__(
             unit_vector_x = self[ii].geometry.unit_vector_along_arm("x")
             unit_vector_y = self[ii].geometry.unit_vector_along_arm("y")
 
-            vertex_geocentric = get_vertex_position_geocentric(self[ii].latitude_radians,
-                                                               self[ii].longitude_radians,
-                                                               self[ii].elevation)
+            vertex_geocentric = get_vertex_position_geocentric(np.array([self[ii].latitude_radians,
+                                                                         self[ii].longitude_radians,
+                                                                         self[ii].elevation]))
             next_vertex_geocentric = vertex_geocentric + length * 1000 * unit_vector_x
-            next_vertex_ellipsoid = get_vertex_position_ellipsoid(next_vertex_geocentric[0],
-                                                                  next_vertex_geocentric[1],
-                                                                  next_vertex_geocentric[2])
+            next_vertex_ellipsoid = get_vertex_position_ellipsoid(next_vertex_geocentric)
             next_latitude_rad, next_longitude_rad, next_elevation = next_vertex_ellipsoid
 
             rotation_vector = np.cross(unit_vector_x, unit_vector_y)
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 9aabc75b5..83944a7ec 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -55,7 +55,7 @@ def psd_from_freq_series(freq_data, df):
     return np.power(asd_from_freq_series(freq_data, df), 2)
 
 
-def get_vertex_position_geocentric(latitude, longitude, elevation):
+def get_vertex_position_geocentric(position, _longitude=None, _elevation=None):
     """
     Calculate the position of the IFO vertex in geocentric coordinates in meters.
 
@@ -64,18 +64,25 @@ 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
 
     """
+
+    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 * np.cos(latitude)**2 +
@@ -86,7 +93,7 @@ def get_vertex_position_geocentric(latitude, longitude, elevation):
     return np.array([x_comp, y_comp, z_comp])
 
 
-def get_vertex_position_ellipsoid(x_comp, y_comp, z_comp):
+def get_vertex_position_ellipsoid(position):
     """
     Calculate the position of the IFO vertex in ellipsoidal coordinates.
 
@@ -97,12 +104,8 @@ def get_vertex_position_ellipsoid(x_comp, y_comp, z_comp):
 
     Parameters
     ==========
-    x_comp: float
-        Geocentric x-coordinate in meters
-    y_comp: float
-        Geocentric y-coordinate in meters
-    z_comp: float
-        Geocentric z-coordinate in meters
+    position: array_like
+        3D representation of the geocentric vertex position in meters (x, y, z)
 
     Returns
     =======
@@ -111,9 +114,9 @@ def get_vertex_position_ellipsoid(x_comp, y_comp, z_comp):
     semi_major_axis = 6378137  # for ellipsoid model of Earth, in m
     semi_minor_axis = 6356752.314  # in m
 
-    r = np.sqrt(x_comp ** 2 + y_comp ** 2)
-    E = (semi_minor_axis * z_comp - (semi_major_axis ** 2 - semi_minor_axis ** 2)) / (r * semi_major_axis)
-    F = (semi_minor_axis * z_comp + (semi_major_axis ** 2 - semi_minor_axis ** 2)) / (r * semi_major_axis)
+    r = np.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
@@ -122,9 +125,9 @@ def get_vertex_position_ellipsoid(x_comp, y_comp, z_comp):
     # Calculate solution in first quadrant and then adjust based on the sign of the original z_comp
     G = 1 / 2 * (np.sqrt(E ** 2 + v) + E)
     t = np.sqrt(G ** 2 + (F - v * G) / (2 * G - E)) - G
-    latitude = np.sign(z_comp) * np.arctan((semi_major_axis * (1 - t ** 2)) / (2 * semi_minor_axis * t))
-    longitude = np.arctan2(y_comp, x_comp)
-    elevation = (r - semi_major_axis * t) * np.cos(latitude) + (z_comp - semi_minor_axis) * np.sin(latitude)
+    latitude = np.sign(position[2]) * np.arctan((semi_major_axis * (1 - t ** 2)) / (2 * semi_minor_axis * t))
+    longitude = np.arctan2(position[1], position[0])
+    elevation = (r - semi_major_axis * t) * np.cos(latitude) + (position[2] - semi_minor_axis) * np.sin(latitude)
 
     return np.array([latitude, longitude, elevation])
 
diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index 384988734..2eb2a0644 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -247,8 +247,8 @@ def test_get_vertex_position_geocentric_ellipsoid_conversion(self):
         latitude_true = np.deg2rad(46 + 27. / 60 + 18.528 / 3600)
         longitude_true = np.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
         elevation_true = 142.554
-        vertex = gwutils.get_vertex_position_geocentric(latitude_true, longitude_true, elevation_true)
-        latitude, longitude, elevation = gwutils.get_vertex_position_ellipsoid(vertex[0], vertex[1], vertex[2])
+        vertex = gwutils.get_vertex_position_geocentric(np.array([latitude_true, longitude_true, elevation_true]))
+        latitude, longitude, elevation = gwutils.get_vertex_position_ellipsoid(vertex)
         self.assertAlmostEqual(latitude, latitude_true, 5)
         self.assertAlmostEqual(longitude, longitude_true, 5)
         self.assertAlmostEqual(elevation, elevation_true, 5)

From 7d654817136963cfb06ada6e340d1f5eb4bf27e8 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Thu, 17 Sep 2026 15:24:19 +0200
Subject: [PATCH 07/15] moved the update logic to a separate function

---
 bilby/gw/detector/networks.py | 115 ++++++++++++++++++++++------------
 1 file changed, 75 insertions(+), 40 deletions(-)

diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index ebd98293b..b1dda598a 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -399,46 +399,81 @@ def __init__(
                 )
             )
 
-            unit_vector_x = self[ii].geometry.unit_vector_along_arm("x")
-            unit_vector_y = self[ii].geometry.unit_vector_along_arm("y")
-
-            vertex_geocentric = get_vertex_position_geocentric(np.array([self[ii].latitude_radians,
-                                                                         self[ii].longitude_radians,
-                                                                         self[ii].elevation]))
-            next_vertex_geocentric = vertex_geocentric + 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 = np.cross(unit_vector_x, unit_vector_y)
-            rotation_vector /= np.linalg.norm(rotation_vector)
-            rotation_angle = 2 / 3 * np.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 = np.array([np.cos(next_latitude_rad) * np.cos(next_longitude_rad),
-                                                 np.cos(next_latitude_rad) * np.sin(next_longitude_rad),
-                                                 np.sin(next_latitude_rad)])
-            next_local_north_vector = np.array([-np.sin(next_latitude_rad) * np.cos(next_longitude_rad),
-                                                -np.sin(next_latitude_rad) * np.sin(next_longitude_rad),
-                                                np.cos(next_latitude_rad)])
-            next_local_east_vector = np.array([-np.sin(next_longitude_rad),
-                                               np.cos(next_longitude_rad), 0])
-
-            next_xarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_x, next_local_normal_vector))
-            next_yarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_y, next_local_normal_vector))
-            next_xarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_x, next_local_north_vector),
-                                               np.dot(next_unit_vector_x, next_local_east_vector))
-            next_yarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_y, next_local_north_vector),
-                                               np.dot(next_unit_vector_y, next_local_east_vector))
-
-            latitude = np.rad2deg(next_latitude_rad)
-            longitude = np.rad2deg(next_longitude_rad)
-            elevation = next_elevation
-            xarm_azimuth = np.rad2deg(next_xarm_azimuth_rad)
-            yarm_azimuth = np.rad2deg(next_yarm_azimuth_rad)
-            xarm_tilt = next_xarm_tilt_rad
-            yarm_tilt = next_yarm_tilt_rad
+
+            latitude, longitude, elevation, xarm_azimuth, yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
+            
+    def _get_next_vertex_parameters(self, current_index):
+        """
+        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.
+        
+        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(np.array([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 = np.cross(unit_vector_x, unit_vector_y)
+        rotation_vector /= np.linalg.norm(rotation_vector)
+        rotation_angle = 2 / 3 * np.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 = np.array([np.cos(next_latitude_rad) * np.cos(next_longitude_rad),
+                                             np.cos(next_latitude_rad) * np.sin(next_longitude_rad),
+                                             np.sin(next_latitude_rad)])
+        next_local_north_vector = np.array([-np.sin(next_latitude_rad) * np.cos(next_longitude_rad),
+                                            -np.sin(next_latitude_rad) * np.sin(next_longitude_rad),
+                                            np.cos(next_latitude_rad)])
+        next_local_east_vector = np.array([-np.sin(next_longitude_rad),
+                                           np.cos(next_longitude_rad), 0])
+
+        next_xarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_x, next_local_normal_vector))
+        next_yarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_y, next_local_normal_vector))
+        next_xarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_x, next_local_north_vector),
+                                           np.dot(next_unit_vector_x, next_local_east_vector))
+        next_yarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_y, next_local_north_vector),
+                                           np.dot(next_unit_vector_y, next_local_east_vector))
+
+        return (np.rad2deg(next_latitude_rad),
+               np.rad2deg(next_longitude_rad),
+               next_elevation,
+               np.rad2deg(next_xarm_azimuth_rad),
+               np.rad2deg(next_yarm_azimuth_rad),
+               next_xarm_tilt_rad,
+               next_yarm_tilt_rad)
 
 def get_empty_interferometer(name):
     """

From e7bc7a687d4ced70464b43b319d7f16342584892 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Thu, 17 Sep 2026 16:36:50 +0200
Subject: [PATCH 08/15] style: pre-commit checks

---
 bilby/gw/detector/networks.py | 29 ++++++++++++++++-------------
 bilby/gw/utils.py             | 11 +++++++----
 test/gw/utils_test.py         |  3 ++-
 3 files changed, 25 insertions(+), 18 deletions(-)

diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index b1dda598a..f386c034d 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -399,22 +399,24 @@ def __init__(
                 )
             )
 
+            latitude, longitude, elevation, xarm_azimuth,
+            yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
 
-            latitude, longitude, elevation, xarm_azimuth, yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
-            
     def _get_next_vertex_parameters(self, current_index):
         """
         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.
+        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.
-        
+
         Returns
         =======
         latitude: float
@@ -468,12 +470,13 @@ def _get_next_vertex_parameters(self, current_index):
                                            np.dot(next_unit_vector_y, next_local_east_vector))
 
         return (np.rad2deg(next_latitude_rad),
-               np.rad2deg(next_longitude_rad),
-               next_elevation,
-               np.rad2deg(next_xarm_azimuth_rad),
-               np.rad2deg(next_yarm_azimuth_rad),
-               next_xarm_tilt_rad,
-               next_yarm_tilt_rad)
+                np.rad2deg(next_longitude_rad),
+                next_elevation,
+                np.rad2deg(next_xarm_azimuth_rad),
+                np.rad2deg(next_yarm_azimuth_rad),
+                next_xarm_tilt_rad,
+                next_yarm_tilt_rad)
+
 
 def get_empty_interferometer(name):
     """
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index f76f867ff..99223631c 100644
--- a/bilby/gw/utils.py
+++ b/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)
@@ -55,6 +55,7 @@ def psd_from_freq_series(freq_data, df, *, xp=None):
     """
     return asd_from_freq_series(freq_data, df, xp=xp) ** 2
 
+
 @xp_wrap
 def get_vertex_position_geocentric(position, _longitude=None, _elevation=None, *, xp=None):
     """
@@ -66,7 +67,7 @@ def get_vertex_position_geocentric(position, _longitude=None, _elevation=None, *
     Parameters
     ==========
     position: array_like
-        3D representation of the ellipsoidal vertex position (latitude [rad], longitude [rad], elevation [m]) 
+        3D representation of the ellipsoidal vertex position (latitude [rad], longitude [rad], elevation [m])
 
     Returns
     =======
@@ -75,11 +76,12 @@ def get_vertex_position_geocentric(position, _longitude=None, _elevation=None, *
     """
 
     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.")
+        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:
+    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.")
@@ -93,6 +95,7 @@ def get_vertex_position_geocentric(position, _longitude=None, _elevation=None, *
     z_comp = ((semi_minor_axis / semi_major_axis)**2 * radius + elevation) * xp.sin(latitude)
     return xp.asarray([x_comp, y_comp, z_comp])
 
+
 @xp_wrap
 def get_vertex_position_ellipsoid(position, *, xp=None):
     """
diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index a8d09e899..45c2555e1 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -281,7 +281,8 @@ def test_get_vertex_position_geocentric_ellipsoid_conversion(self):
         latitude_true = np.deg2rad(46 + 27. / 60 + 18.528 / 3600)
         longitude_true = np.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
         elevation_true = 142.554
-        vertex = gwutils.get_vertex_position_geocentric(self.xp.asarray([latitude_true, longitude_true, elevation_true]))
+        vertex = gwutils.get_vertex_position_geocentric(
+            self.xp.asarray([latitude_true, longitude_true, elevation_true]))
         latitude, longitude, elevation = gwutils.get_vertex_position_ellipsoid(vertex)
         self.assertAlmostEqual(latitude, latitude_true, 5)
         self.assertAlmostEqual(longitude, longitude_true, 5)

From 3ab8e06d0773442afcb3e4594872e7edb5e88bc9 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Thu, 17 Sep 2026 18:02:15 +0200
Subject: [PATCH 09/15] fix: syntax error for line continuation

---
 bilby/gw/detector/networks.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index f386c034d..cf7a39c77 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -399,8 +399,8 @@ def __init__(
                 )
             )
 
-            latitude, longitude, elevation, xarm_azimuth,
-            yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
+            latitude, longitude, elevation, xarm_azimuth, \
+                yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
 
     def _get_next_vertex_parameters(self, current_index):
         """

From 90816c6ef7d00aaa4922c08ab040fa16f6eccd4a Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Wed, 30 Sep 2026 15:26:45 +0200
Subject: [PATCH 10/15] array_backend: removed numpy references in triangular
 ifo code

To infer the correct array backend in `InterferometerGeometry.vertex`,
the property `array_backend` was added in line with the implementations
in `Interferometer` and `InterferometerList`. This also makes the
`Interferometer.array_backend` method independent of implementation
details in `InterferometerGeometry`.
---
 bilby/gw/detector/geometry.py       | 10 +++++-
 bilby/gw/detector/interferometer.py |  2 +-
 bilby/gw/detector/networks.py       | 53 +++++++++++++++--------------
 bilby/gw/utils.py                   | 14 ++++----
 4 files changed, 44 insertions(+), 35 deletions(-)

diff --git a/bilby/gw/detector/geometry.py b/bilby/gw/detector/geometry.py
index 0d7167956..e900111d9 100644
--- a/bilby/gw/detector/geometry.py
+++ b/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,8 +208,11 @@ 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(np.array([self._latitude, self._longitude,
+            self._vertex = gwutils.get_vertex_position_geocentric(xp.array([self._latitude, self._longitude,
                                                                   self.elevation]))
             self._vertex_updated = True
         return self._vertex
@@ -314,3 +318,7 @@ def set_array_backend(self, xp):
         self.yarm_azimuth = xp.asarray(self.yarm_azimuth)
         self.xarm_tilt = xp.asarray(self.xarm_tilt)
         self.yarm_tilt = xp.asarray(self.yarm_tilt)
+
+    @property
+    def array_backend(self):
+        return array_module(self.length)
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index 5f181539c..42c275875 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -945,4 +945,4 @@ def set_array_backend(self, xp):
 
     @property
     def array_backend(self):
-        return array_module(self.geometry.length)
+        return self.geometry.array_backend
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index cf7a39c77..fc0d8989a 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -1,9 +1,9 @@
 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
@@ -402,7 +402,8 @@ def __init__(
             latitude, longitude, elevation, xarm_azimuth, \
                 yarm_azimuth, xarm_tilt, yarm_tilt = self._get_next_vertex_parameters(ii)
 
-    def _get_next_vertex_parameters(self, current_index):
+    @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
@@ -439,41 +440,41 @@ def _get_next_vertex_parameters(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(np.array([current_ifo.latitude_radians,
+        vertex_geocentric = get_vertex_position_geocentric(xp.array([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 = np.cross(unit_vector_x, unit_vector_y)
-        rotation_vector /= np.linalg.norm(rotation_vector)
-        rotation_angle = 2 / 3 * np.pi
+        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 = np.array([np.cos(next_latitude_rad) * np.cos(next_longitude_rad),
-                                             np.cos(next_latitude_rad) * np.sin(next_longitude_rad),
-                                             np.sin(next_latitude_rad)])
-        next_local_north_vector = np.array([-np.sin(next_latitude_rad) * np.cos(next_longitude_rad),
-                                            -np.sin(next_latitude_rad) * np.sin(next_longitude_rad),
-                                            np.cos(next_latitude_rad)])
-        next_local_east_vector = np.array([-np.sin(next_longitude_rad),
-                                           np.cos(next_longitude_rad), 0])
-
-        next_xarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_x, next_local_normal_vector))
-        next_yarm_tilt_rad = np.arcsin(np.dot(next_unit_vector_y, next_local_normal_vector))
-        next_xarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_x, next_local_north_vector),
-                                           np.dot(next_unit_vector_x, next_local_east_vector))
-        next_yarm_azimuth_rad = np.arctan2(np.dot(next_unit_vector_y, next_local_north_vector),
-                                           np.dot(next_unit_vector_y, next_local_east_vector))
-
-        return (np.rad2deg(next_latitude_rad),
-                np.rad2deg(next_longitude_rad),
+        next_local_normal_vector = xp.array([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.array([-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.array([-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,
-                np.rad2deg(next_xarm_azimuth_rad),
-                np.rad2deg(next_yarm_azimuth_rad),
+                xp.rad2deg(next_xarm_azimuth_rad),
+                xp.rad2deg(next_yarm_azimuth_rad),
                 next_xarm_tilt_rad,
                 next_yarm_tilt_rad)
 
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 99223631c..44e8c4712 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -118,20 +118,20 @@ def get_vertex_position_ellipsoid(position, *, xp=None):
     semi_major_axis = 6378137  # for ellipsoid model of Earth, in m
     semi_minor_axis = 6356752.314  # in m
 
-    r = np.sqrt(position[0] ** 2 + position[1] ** 2)
+    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 = (np.sqrt(D) - Q) ** (1 / 3) - (np.sqrt(D) + Q) ** (1 / 3)
+    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 * (np.sqrt(E ** 2 + v) + E)
-    t = np.sqrt(G ** 2 + (F - v * G) / (2 * G - E)) - G
-    latitude = np.sign(position[2]) * np.arctan((semi_major_axis * (1 - t ** 2)) / (2 * semi_minor_axis * t))
-    longitude = np.arctan2(position[1], position[0])
-    elevation = (r - semi_major_axis * t) * np.cos(latitude) + (position[2] - semi_minor_axis) * np.sin(latitude)
+    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])
 

From 60b2c462197ed27dd8f19e4c00dd745801373d4e Mon Sep 17 00:00:00 2001
From: Milan Wils <58658688+milanwils@users.noreply.github.com>
Date: Wed, 30 Sep 2026 17:04:57 +0200
Subject: [PATCH 11/15] Apply suggestion from @ColmTalbot about remaining np
 function

Co-authored-by: Colm Talbot <talbotcolm@gmail.com>
---
 bilby/gw/detector/interferometer.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index 42c275875..c53b2ad46 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -586,7 +586,9 @@ def vertex_position_geocentric(self):
         =======
         array_like: A 3D array representation of the vertex
         """
-        return gwutils.get_vertex_position_geocentric(np.array([self.geometry.latitude_radians,
+        xp = self.array_backend
+
+        return gwutils.get_vertex_position_geocentric(xp.array([self.geometry.latitude_radians,
                                                                 self.geometry.longitude_radians,
                                                                 self.geometry.elevation]))
 

From 684734779371a56aa7ea0057b510c1714d239212 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Wed, 30 Sep 2026 17:07:08 +0200
Subject: [PATCH 12/15] changed two remaining np calls to xp

---
 test/gw/utils_test.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index 45c2555e1..5a7cd07ec 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -278,8 +278,8 @@ def test_lalsim_SimInspiralChooseFDWaveform(self):
             )
 
     def test_get_vertex_position_geocentric_ellipsoid_conversion(self):
-        latitude_true = np.deg2rad(46 + 27. / 60 + 18.528 / 3600)
-        longitude_true = np.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
+        latitude_true = self.xp.deg2rad(46 + 27. / 60 + 18.528 / 3600)
+        longitude_true = self.xp.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
         elevation_true = 142.554
         vertex = gwutils.get_vertex_position_geocentric(
             self.xp.asarray([latitude_true, longitude_true, elevation_true]))

From 265331850c4dcec73aea4504f79a34ac95042058 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Wed, 30 Sep 2026 17:38:24 +0200
Subject: [PATCH 13/15] refactor: replace set_array_backend methods with
 @array_backend.setter property

---
 bilby/gw/detector/geometry.py                         | 11 ++++++-----
 bilby/gw/detector/interferometer.py                   |  9 +++++----
 bilby/gw/detector/networks.py                         |  9 +++++----
 bilby/gw/detector/psd.py                              |  8 +++++++-
 bilby/gw/likelihood/base.py                           |  6 +++---
 docs/array_api.rst                                    |  4 ++--
 .../injection_examples/jax_fast_tutorial.py           |  2 +-
 test/gw/detector/geometry_test.py                     |  6 +++---
 test/gw/likelihood/base_test.py                       |  2 +-
 test/gw/likelihood/basic_test.py                      |  2 +-
 test/gw/likelihood/multiband_test.py                  |  2 +-
 test/gw/likelihood/roq_test.py                        |  8 ++++----
 test/gw/utils_test.py                                 |  4 ++--
 13 files changed, 41 insertions(+), 32 deletions(-)

diff --git a/bilby/gw/detector/geometry.py b/bilby/gw/detector/geometry.py
index e900111d9..d5fe3e04c 100644
--- a/bilby/gw/detector/geometry.py
+++ b/bilby/gw/detector/geometry.py
@@ -309,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)
@@ -318,7 +323,3 @@ def set_array_backend(self, xp):
         self.yarm_azimuth = xp.asarray(self.yarm_azimuth)
         self.xarm_tilt = xp.asarray(self.xarm_tilt)
         self.yarm_tilt = xp.asarray(self.yarm_tilt)
-
-    @property
-    def array_backend(self):
-        return array_module(self.length)
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index c53b2ad46..bd0075ef0 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -941,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 self.geometry.array_backend
+
+    @array_backend.setter
+    def array_backend(self, xp):
+        self.geometry.array_backend = xp
+        self.power_spectral_density.array_backend = xp
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index fc0d8989a..2b010d204 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -345,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__(
diff --git a/bilby/gw/detector/psd.py b/bilby/gw/detector/psd.py
index aab53b828..864dbd380 100644
--- a/bilby/gw/detector/psd.py
+++ b/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
diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py
index 0030716ae..3f69a10a9 100644
--- a/bilby/gw/likelihood/base.py
+++ b/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):
         """
diff --git a/docs/array_api.rst b/docs/array_api.rst
index 02cb5c2c2..7b914e525 100644
--- a/docs/array_api.rst
+++ b/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,
diff --git a/examples/gw_examples/injection_examples/jax_fast_tutorial.py b/examples/gw_examples/injection_examples/jax_fast_tutorial.py
index 22a64c935..bfb7f1b08 100644
--- a/examples/gw_examples/injection_examples/jax_fast_tutorial.py
+++ b/examples/gw_examples/injection_examples/jax_fast_tutorial.py
@@ -181,7 +181,7 @@ def main():
         parameters=injection_parameters,
         raise_error=False,
     )
-    ifos.set_array_backend(jnp)
+    ifos.array_backend = jnp
 
     # Initialise the likelihood by passing in the interferometer data (ifos) and
     # the waveform generator
diff --git a/test/gw/detector/geometry_test.py b/test/gw/detector/geometry_test.py
index 7340a5f8d..a7c27bb59 100644
--- a/test/gw/detector/geometry_test.py
+++ b/test/gw/detector/geometry_test.py
@@ -30,7 +30,7 @@ def setUp(self):
             xarm_tilt=self.xarm_tilt,
             yarm_tilt=self.yarm_tilt,
         )
-        self.geometry.set_array_backend(self.xp)
+        self.geometry.array_backend = self.xp
 
     def tearDown(self):
         del self.length
@@ -196,7 +196,7 @@ def test_unit_vector_along_arm_x(self):
         self.geometry.latitude = 0
         self.geometry.xarm_tilt = 0
         self.geometry.xarm_azimuth = 0
-        self.geometry.set_array_backend(self.xp)
+        self.geometry.array_backend = self.xp
         arm = self.geometry.unit_vector_along_arm("x")
         self.assertTrue(np.allclose(arm, np.array([0, 1, 0])))
         self.assertEqual(aac.get_namespace(arm), self.xp)
@@ -206,7 +206,7 @@ def test_unit_vector_along_arm_y(self):
         self.geometry.latitude = 0
         self.geometry.yarm_tilt = 0
         self.geometry.yarm_azimuth = 90
-        self.geometry.set_array_backend(self.xp)
+        self.geometry.array_backend = self.xp
         arm = self.geometry.unit_vector_along_arm("y")
         self.assertTrue(np.allclose(arm, np.array([0, 0, 1])))
         self.assertEqual(aac.get_namespace(arm), self.xp)
diff --git a/test/gw/likelihood/base_test.py b/test/gw/likelihood/base_test.py
index 21d50cd64..ff0e7c130 100644
--- a/test/gw/likelihood/base_test.py
+++ b/test/gw/likelihood/base_test.py
@@ -36,7 +36,7 @@ def setUp(self):
         self.interferometers.set_strain_data_from_power_spectral_densities(
             sampling_frequency=self.sampling_frequency, duration=self.duration
         )
-        self.interferometers.set_array_backend(self.xp)
+        self.interferometers.array_backend = self.xp
         wfg = bilby.gw.waveform_generator.WaveformGenerator(
             duration=self.duration,
             sampling_frequency=self.sampling_frequency,
diff --git a/test/gw/likelihood/basic_test.py b/test/gw/likelihood/basic_test.py
index 10751a271..516ecbc32 100644
--- a/test/gw/likelihood/basic_test.py
+++ b/test/gw/likelihood/basic_test.py
@@ -34,7 +34,7 @@ def setUp(self):
         self.interferometers.set_strain_data_from_power_spectral_densities(
             sampling_frequency=self.xp.asarray(2048.0), duration=self.xp.asarray(4.0)
         )
-        self.interferometers.set_array_backend(self.xp)
+        self.interferometers.array_backend = self.xp
         base_wfg = bilby.gw.waveform_generator.GWSignalWaveformGenerator(
             duration=self.xp.asarray(4.0),
             sampling_frequency=self.xp.asarray(2048.0),
diff --git a/test/gw/likelihood/multiband_test.py b/test/gw/likelihood/multiband_test.py
index 2e0330462..33954d8bc 100644
--- a/test/gw/likelihood/multiband_test.py
+++ b/test/gw/likelihood/multiband_test.py
@@ -47,7 +47,7 @@ def setUp(self):
         )
         for ifo in self.ifos:
             ifo.minimum_frequency = self.fmin
-        self.ifos.set_array_backend(self.xp)
+        self.ifos.array_backend = self.xp
 
         spline_calibration_nodes = 10
         self.calibration_parameters = {}
diff --git a/test/gw/likelihood/roq_test.py b/test/gw/likelihood/roq_test.py
index 34fa5504f..5e46e0e59 100644
--- a/test/gw/likelihood/roq_test.py
+++ b/test/gw/likelihood/roq_test.py
@@ -65,7 +65,7 @@ def setUp(self):
         ifos.set_strain_data_from_power_spectral_densities(
             sampling_frequency=self.sampling_frequency, duration=self.duration
         )
-        ifos.set_array_backend(self.xp)
+        ifos.array_backend = self.xp
 
         self.priors = bilby.gw.prior.BBHPriorDict()
         self.priors.pop("mass_1")
@@ -423,7 +423,7 @@ def test_fails_with_frequency_duration_mismatch(
             duration=self.xp.asarray(duration),
             start_time=self.xp.asarray(self.injection_parameters["geocent_time"] - duration + 1)
         )
-        interferometers.set_array_backend(self.xp)
+        interferometers.array_backend = self.xp
         for ifo in interferometers:
             ifo.minimum_frequency = self.xp.asarray(minimum_frequency)
             ifo.maximum_frequency = self.xp.asarray(maximum_frequency)
@@ -457,7 +457,7 @@ def test_fails_with_prior_mismatch(self, basis, chirp_mass_min, chirp_mass_max):
             duration=self.xp.asarray(self.duration),
             start_time=self.xp.asarray(self.injection_parameters["geocent_time"] - self.duration + 1)
         )
-        interferometers.set_array_backend(self.xp)
+        interferometers.array_backend = self.xp
         for ifo in interferometers:
             ifo.minimum_frequency = self.minimum_frequency
         search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
@@ -596,7 +596,7 @@ def assertLess_likelihood_errors(
         self.priors["chirp_mass"].maximum = mc_max
 
         interferometers = bilby.gw.detector.InterferometerList(["H1", "L1"])
-        interferometers.set_array_backend(self.xp)
+        interferometers.array_backend = self.xp
         for ifo in interferometers:
             if minimum_frequency is None:
                 ifo.minimum_frequency = self.minimum_frequency
diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index 5a7cd07ec..7a9d6f541 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -333,7 +333,7 @@ def test_conversion_single(self) -> None:
         zenith = self.xp.asarray(sample["zenith"])
         azimuth = self.xp.asarray(sample["azimuth"])
         time = self.xp.asarray(sample["time"])
-        self.ifos.set_array_backend(self.xp)
+        self.ifos.array_backend = self.xp
         ra, dec = bilby.gw.utils.zenith_azimuth_to_ra_dec(
             zenith, azimuth, time, self.ifos
         )
@@ -344,7 +344,7 @@ def test_conversion_gives_correct_prior(self) -> None:
         zeniths = self.xp.asarray(self.samples["zenith"])
         azimuths = self.xp.asarray(self.samples["azimuth"])
         times = self.xp.asarray(self.samples["time"])
-        self.ifos.set_array_backend(self.xp)
+        self.ifos.array_backend = self.xp
         ras, decs = bilby.gw.utils.zenith_azimuth_to_ra_dec(
             zeniths, azimuths, times, self.ifos
         )

From b72bc0f52e62a5f4778a4c630b5472d9a8cc59f7 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Mon, 5 Oct 2026 10:20:05 +0200
Subject: [PATCH 14/15] fix: changed xp.array to xp.asarray for torch
 compatibility

---
 bilby/gw/detector/geometry.py       |  2 +-
 bilby/gw/detector/interferometer.py |  6 +++---
 bilby/gw/detector/networks.py       | 22 +++++++++++-----------
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/bilby/gw/detector/geometry.py b/bilby/gw/detector/geometry.py
index d5fe3e04c..f2461761e 100644
--- a/bilby/gw/detector/geometry.py
+++ b/bilby/gw/detector/geometry.py
@@ -212,7 +212,7 @@ def vertex(self):
         xp = self.array_backend
 
         if not self._vertex_updated:
-            self._vertex = gwutils.get_vertex_position_geocentric(xp.array([self._latitude, self._longitude,
+            self._vertex = gwutils.get_vertex_position_geocentric(xp.asarray([self._latitude, self._longitude,
                                                                   self.elevation]))
             self._vertex_updated = True
         return self._vertex
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index bd0075ef0..4954c17f8 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -588,9 +588,9 @@ def vertex_position_geocentric(self):
         """
         xp = self.array_backend
 
-        return gwutils.get_vertex_position_geocentric(xp.array([self.geometry.latitude_radians,
-                                                                self.geometry.longitude_radians,
-                                                                self.geometry.elevation]))
+        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):
         """
diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py
index 2b010d204..a5e65259a 100644
--- a/bilby/gw/detector/networks.py
+++ b/bilby/gw/detector/networks.py
@@ -441,9 +441,9 @@ def _get_next_vertex_parameters(self, current_index, *, xp):
         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.array([current_ifo.latitude_radians,
-                                                                     current_ifo.longitude_radians,
-                                                                     current_ifo.elevation]))
+        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
@@ -455,14 +455,14 @@ def _get_next_vertex_parameters(self, current_index, *, xp):
         next_unit_vector_x = rotation.apply(unit_vector_x)
         next_unit_vector_y = rotation.apply(unit_vector_y)
 
-        next_local_normal_vector = xp.array([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.array([-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.array([-xp.sin(next_longitude_rad),
-                                           xp.cos(next_longitude_rad), 0])
+        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))

From 8ac76beebd6f85242b7734d71bf4da1e8bab9678 Mon Sep 17 00:00:00 2001
From: Milan Wils <milan.wils@kuleuven.be>
Date: Tue, 6 Oct 2026 13:51:29 +0200
Subject: [PATCH 15/15] fix: gw/utils tests failing for pytorch

One test was failing in `gw/utils.py` because `torch.deg2rad` only
accepts Tensor and not a float. This is fixed by turning the argument
into an array using `xp.asarray`.

Additionally, the test used assertAlmostEqual which requires the builtin
function `__round__` which is not defined for `Tensor` objects.
Therefore, the tensors are cast into float before comparing.
---
 test/gw/utils_test.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py
index 7a9d6f541..e63bdc88a 100644
--- a/test/gw/utils_test.py
+++ b/test/gw/utils_test.py
@@ -278,15 +278,15 @@ def test_lalsim_SimInspiralChooseFDWaveform(self):
             )
 
     def test_get_vertex_position_geocentric_ellipsoid_conversion(self):
-        latitude_true = self.xp.deg2rad(46 + 27. / 60 + 18.528 / 3600)
-        longitude_true = self.xp.deg2rad(-(119 + 24. / 60 + 27.5657 / 3600))
+        latitude_true = self.xp.deg2rad(self.xp.asarray(46 + 27. / 60 + 18.528 / 3600))
+        longitude_true = self.xp.deg2rad(self.xp.asarray(-(119 + 24. / 60 + 27.5657 / 3600)))
         elevation_true = 142.554
         vertex = gwutils.get_vertex_position_geocentric(
             self.xp.asarray([latitude_true, longitude_true, elevation_true]))
         latitude, longitude, elevation = gwutils.get_vertex_position_ellipsoid(vertex)
-        self.assertAlmostEqual(latitude, latitude_true, 5)
-        self.assertAlmostEqual(longitude, longitude_true, 5)
-        self.assertAlmostEqual(elevation, elevation_true, 5)
+        self.assertAlmostEqual(float(latitude), float(latitude_true), 5)
+        self.assertAlmostEqual(float(longitude), float(longitude_true), 5)
+        self.assertAlmostEqual(float(elevation), float(elevation_true), 5)
 
     def test_safe_cast_mode_to_int(self):
         # Valid cases