Skip to content

Align Einstein Telescope geometry with LAL #894

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
milanwils wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Align Einstein Telescope geometry with LAL #894

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a5f498e
naming ET detectors anti-clockwise
milan-wils Jan 13, 2025
ca46c1e
Implemented transform ECEF to WGS84
milan-wils Jan 16, 2025
3fa49f0
added test for ECEF to WGS84 conversion
milan-wils Jan 16, 2025
058d588
Changed triangular config in line with LAL definition of ET
milan-wils Jan 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions bilby/gw/detector/detectors/ET.interferometer
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -12,3 +13,4 @@ elevation = 51.884
xarm_azimuth = 70.5674
yarm_azimuth = 130.5674
shape = 'Triangle'
clockwise = False
69 changes: 46 additions & 23 deletions bilby/gw/detector/networks.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -310,6 +312,7 @@ def __init__(
yarm_azimuth,
xarm_tilt=0.0,
yarm_tilt=0.0,
clockwise=True,
Copy link
Copy Markdown
Collaborator

@ColmTalbot ColmTalbot Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would personally vote for just making a breaking change here (and waiting for a major release) to avoid adding a new parameter that doesn't really mean anything (the set of three Interferometers are just a permutation of each other with the different directions).

Copy link
Copy Markdown
Author

@milanwils milanwils Feb 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine by me, if people really want they can still override the detector names

):
super(TriangularInterferometer, self).__init__([])
self.name = name
Expand All @@ -322,6 +325,7 @@ def __init__(
maximum_frequency = [maximum_frequency] * 3

for ii in range(3):

self.append(
Interferometer(
"{}{}".format(name, ii + 1),
Expand All @@ -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")
Copy link
Copy Markdown
Collaborator

@ColmTalbot ColmTalbot Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this is not the change for the current code, but given that this is clearly better thought out than the current method, can the logic for updating the geometry be moved into a standalone function (with a nice docstring)?

Copy link
Copy Markdown
Author

@milanwils milanwils Feb 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps we could add setters to the InterferometerGeometry class for x, y and vertex? For the lat, long, elevation, tilt and azimuth these setters are separate but it does not seem sensible to me that someone can change only one Cartesian coordinate.

All the update logic could then be moved into the geometry class and only the calculation of the new corner points (which is triangle specific) remains in the 'TriangularInterferometer` class.

My main concern with this is approach is that one would need to be careful to update the vertex position before changing the unit vector. Otherwise the calculation of tilt and azimuth is wrong. This could be avoided by making the unit vectors the "base information" and deriving the tilt and azimuth from it.

This also raises the question of whether the constructor should be adapted to allow users to define the detector directly from Cartesian parameters. I guess you could argue that this is actually the only way you should be able to set these parameters as interferometers are constant.

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):
Expand Down
43 changes: 43 additions & 0 deletions bilby/gw/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Copy link
Copy Markdown
Collaborator

@ColmTalbot ColmTalbot Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like this function should take and recieve the fame type, either (float, float, float) -> (float, float, float) or array -> array. I don't really mind which, but it would be good to be consistent.

I now see that the same thing happens in get_vertex_position_geocentric. sigh, I'm not sure if it is worth changing the other function.

We could change the other function to something like

def get_vertex_position_geocentric(position, _longitude=None, _elevation=None):
    if longitude is not None and not isarray(position):
        WARN("... syntax changed in Bilby 3.0.0 ...")
        position = np.array([position, _longitude, _elevation])
    ...

This would need a wider set of eyes/opinions.

Copy link
Copy Markdown
Author

@milanwils milanwils Feb 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I went through the exact same thought process when writing the function. This seems like a good compromise.

"""
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
Expand Down
10 changes: 10 additions & 0 deletions test/gw/utils_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):

Expand Down