From a5f498e208e3b50c645873b635d55e61ec6611b1 Mon Sep 17 00:00:00 2001 From: Milan Wils Date: Mon, 13 Jan 2025 13:40:36 +0100 Subject: [PATCH 1/4] 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 Date: Thu, 16 Jan 2025 11:32:50 +0100 Subject: [PATCH 2/4] 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 Date: Thu, 16 Jan 2025 15:22:33 +0100 Subject: [PATCH 3/4] 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 Date: Fri, 17 Jan 2025 17:36:41 +0100 Subject: [PATCH 4/4] 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):