add in altitude support and simplify adding detector function (#4884)

* add in altitude support and simplify adding detector function

* update

* update unittest

pycbc/detector.py

@@ -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
@@ -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
@@ -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]))
 
     full_resp = (resps[0] - resps[1])
     loc = coordinates.EarthLocation.from_geodetic(longitude * units.rad,
@@ -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,
                               }
@@ -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(':'))
@@ -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))
 
@@ -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
+
+        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:

test/test_detector.py

@@ -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: