Skip to content
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

Jump to bottom

add in altitude support and simplify adding detector function #4884

Merged
merged 3 commits into from
Oct 2, 2024
Merged

add in altitude support and simplify adding detector function #4884

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
464bcad
add in altitude support and simplify adding detector function
ahnitz Sep 19, 2024
c516405
update
ahnitz Oct 1, 2024
ed43e04
update unittest
ahnitz Oct 1, 2024
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
102 changes: 56 additions & 46 deletions pycbc/detector.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,8 @@ def get_available_lal_detectors():

def add_detector_on_earth(name, longitude, latitude,
yangle=0, xangle=None, height=0,
xlength=4000, ylength=4000):
xlength=4000, ylength=4000,
xaltitude=0, yaltitude=0):
""" Add a new detector on the earth

Parameters
Expand All @@ -98,6 +99,10 @@ def add_detector_on_earth(name, longitude, latitude,
xangle: float
Azimuthal angle of the x-arm (angle drawn from point north). If not set
we assume a right angle detector following the right-hand rule.
xaltitude: float
The altitude angle of the x-arm measured from the local horizon.
yaltitude: float
The altitude angle of the y-arm measured from the local horizon.
height: float
The height in meters of the detector above the standard
reference ellipsoidal earth
Expand All @@ -106,27 +111,23 @@ def add_detector_on_earth(name, longitude, latitude,
# assume right angle detector if no separate xarm direction given
xangle = yangle + np.pi / 2.0

# Rotation matrix to move detector to correct orientation
rm1 = rotation_matrix(longitude * units.rad, 'z')
rm2 = rotation_matrix((np.pi / 2.0 - latitude) * units.rad, 'y')
rm = np.matmul(rm2, rm1)

# baseline response of a single arm pointed in the -X direction
resp = np.array([[-1, 0, 0], [0, 0, 0], [0, 0, 0]])
rm2 = rotation_matrix(-longitude * units.rad, 'z')
rm1 = rotation_matrix(-1.0 * (np.pi / 2.0 - latitude) * units.rad, 'y')
# Calculate response in earth centered coordinates
# by rotation of response in coordinates aligned
# with the detector arms
resps = []
vecs = []
for angle in [yangle, xangle]:
a, b = cos(2 * angle), sin(2 * angle)
resp = np.array([[-a, b, 0], [b, a, 0], [0, 0, 0]])

for angle, azi in [(yangle, yaltitude), (xangle, xaltitude)]:
rm0 = rotation_matrix(angle * units.rad, 'z')
rmN = rotation_matrix(-azi * units.rad, 'y')
rm = rm2 @ rm1 @ rm0 @ rmN
# apply rotation
resp = np.matmul(resp, rm)
resp = np.matmul(rm.T, resp) / 4.0
resps.append(resp)

vec = np.matmul(rm.T, np.array([-np.cos(angle), np.sin(angle), 0]))
vecs.append(vec)
resps.append(rm @ resp @ rm.T / 2.0)
vecs.append(rm @ np.array([-1, 0, 0]))
Copy link
Contributor

Choose a reason for hiding this comment

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

So there's changes in here. I think to include the "altitude angles" of the two arms. Have you checked/can you show that this is unchanged in the limit that the altitude angles are both 0.

For current instruments (H1,L1,V1,K1), will this make any actual change?

Copy link
Member Author

Choose a reason for hiding this comment

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

@spxiwh That's actually what I added as a unit test. The key thing is that this sets up the response matrix. I numerically compare the response matrix here against the cached ones in lal. There is no change to the response there (except the E* case which has a 10^-4 level change due to the issue I raised in slack).

The altitudes are near zero for the H1/ L1 / V1 comparison, but not exactly so (nor should they be). In fact if you don't account for the altitude the response matrix would also be different at the O(10^-4) level like E* is.


full_resp = (resps[0] - resps[1])
loc = coordinates.EarthLocation.from_geodetic(longitude * units.rad,
Expand All @@ -142,8 +143,8 @@ def add_detector_on_earth(name, longitude, latitude,
'yangle': yangle,
'xangle': xangle,
'height': height,
'xaltitude': 0.0,
'yaltitude': 0.0,
'xaltitude': xaltitude,
'yaltitude': yaltitude,
'ylength': ylength,
'xlength': xlength,
}
Expand Down Expand Up @@ -191,6 +192,22 @@ def load_detector_config(config_files):
method(det.upper(), *args, **kwds)


# prepopulate using detectors hardcoded into lalsuite
for pref, name in get_available_lal_detectors():
lalsim = pycbc.libutils.import_optional('lalsimulation')
lal_det = lalsim.DetectorPrefixToLALDetector(pref).frDetector
add_detector_on_earth(pref,
lal_det.vertexLongitudeRadians,
lal_det.vertexLatitudeRadians,
height=lal_det.vertexElevation,
xangle=lal_det.xArmAzimuthRadians,
yangle=lal_det.yArmAzimuthRadians,
xlength=lal_det.xArmMidpoint * 2,
ylength=lal_det.yArmMidpoint * 2,
xaltitude=lal_det.xArmAltitudeRadians,
yaltitude=lal_det.yArmAltitudeRadians,
)

# autoload detector config files
if 'PYCBC_DETECTOR_CONFIG' in os.environ:
load_detector_config(os.environ['PYCBC_DETECTOR_CONFIG'].split(':'))
Expand Down Expand Up @@ -218,11 +235,6 @@ def __init__(self, detector_name, reference_time=1126259462.0):
self.info = _ground_detectors[detector_name]
self.response = self.info['response']
self.location = self.info['location']
elif detector_name in lal_detectors:
lalsim = pycbc.libutils.import_optional('lalsimulation')
self._lal = lalsim.DetectorPrefixToLALDetector(self.name)
self.response = self._lal.response
self.location = self._lal.location
else:
raise ValueError("Unkown detector {}".format(detector_name))

Expand All @@ -247,29 +259,27 @@ def set_gmst_reference(self):

def lal(self):
""" Return lal data type detector instance """
if hasattr(self, '_lal'):
return self._lal
else:
import lal
d = lal.FrDetector()
d.vertexLongitudeRadians = self.longitude
d.vertexLatitudeRadians = self.latitude
d.vertexElevation = self.info['height']
d.xArmAzimuthRadians = self.info['xangle']
d.yArmAzimuthRadians = self.info['yangle']
d.xArmAltitudeRadians = self.info['yaltitude']
d.xArmAltitudeRadians = self.info['xaltitude']

# This is somewhat abused by lalsimulation at the moment
# to determine a filter kernel size. We set this only so that
# value gets a similar number of samples as other detectors
# it is used for nothing else
d.yArmMidpoint = 4000.0

x = lal.Detector()
r = lal.CreateDetector(x, d, lal.LALDETECTORTYPE_IFODIFF)
self._lal = r
return r
import lal
d = lal.FrDetector()
d.vertexLongitudeRadians = self.longitude
d.vertexLatitudeRadians = self.latitude
d.vertexElevation = self.info['height']
d.xArmAzimuthRadians = self.info['xangle']
d.yArmAzimuthRadians = self.info['yangle']
d.xArmAltitudeRadians = self.info['xaltitude']
d.yArmAltitudeRadians = self.info['yaltitude']

# This is somewhat abused by lalsimulation at the moment
# to determine a filter kernel size. We set this only so that
# value gets a similar number of samples as other detectors
# it is used for nothing else
d.yArmMidpoint = self.info['ylength'] / 2.0
d.xArmMidpoint = self.info['xlength'] / 2.0
Copy link
Contributor

Choose a reason for hiding this comment

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

Why is this being changed?

Copy link
Member Author

Choose a reason for hiding this comment

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

@spxiwh The lal code doesn't actually use these numbers for much, but it does have them set properly to the right values. This change means that they will be set in the same way as the internal cached detectors in lal in case lal every does use them (for example if someone were to add high-frequency corrections to lalsuite it might matter). Before, I was not bothering to really set them because in the code itself the only thing it used them for was was a guess to a filter length.


x = lal.Detector()
r = lal.CreateDetector(x, d, lal.LALDETECTORTYPE_IFODIFF)
self._lal = r
return r

def gmst_estimate(self, gps_time):
if self.reference_time is None:
Expand Down
18 changes: 18 additions & 0 deletions test/test_detector.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,24 @@ def test_light_time(self):
t2 = d1.light_travel_time_to_detector(d2)
self.assertAlmostEqual(t1, t2, 7)

def test_custom_detector(self):
det.add_detector_on_earth("TEST", 1.3, 0.5, yangle=0, xaltitude=0.01)
d = det.Detector("TEST")

# Check that we can call the new detector response
fp, fc = d.antenna_pattern(1.5, 1.0, 0, 1000000000)

# Check it interacts with existing detectors
d2 = det.Detector("H1")
t1 = d2.light_travel_time_to_detector(d)

def test_response_matrix(self):
import lal
for ifo in ['H1', 'L1', 'V1', 'K1', 'I1']:
ref_resp = lal.cached_detector_by_prefix[ifo].response
resp = det.Detector(ifo).response
self.assertAlmostEqual((ref_resp - resp).max(), 0, places=6)

def test_antenna_pattern(self):
vals = list(zip(self.ra, self.dec, self.pol, self.time))
for ifo in self.d:
Expand Down
Loading