diff --git a/bilby/gw/detector/detectors/ET.interferometer b/bilby/gw/detector/detectors/ET.interferometer index e89384095..79b81d1d0 100644 --- a/bilby/gw/detector/detectors/ET.interferometer +++ b/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 @@ -12,3 +13,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..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): @@ -310,6 +312,7 @@ def __init__( yarm_azimuth, xarm_tilt=0.0, yarm_tilt=0.0, + clockwise=True, ): super(TriangularInterferometer, self).__init__([]) self.name = name @@ -322,6 +325,7 @@ def __init__( maximum_frequency = [maximum_frequency] * 3 for ii in range(3): + self.append( Interferometer( "{}{}".format(name, ii + 1), @@ -339,29 +343,48 @@ def __init__( ) ) - xarm_azimuth += 240 - yarm_azimuth += 240 - - 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): 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 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):