Created a sub method to compute the necessary reflected regions separation angle…

…. As the necessary separation angle depends on the rotation angle but also on the rectangle extension (overlaps may easily occur) the computation is not so easy: I implemented a rough (slightly overestimating) algorithm which chooses the maximum opening angle between the corners of the rectangle with respect to the pointing position. Extended implementation of _reflected_regions() method to allow the placement of rectangular off regions, which is essentially finally placing a rectangle instead of a circle.

Created a sub method to compute the necessary reflected regions separation angle…
…. As the necessary separation angle depends on the rotation angle but also on the rectangle extension (overlaps may easily occur) the computation is not so easy: I implemented a rough (slightly overestimating) algorithm which chooses the maximum opening angle between the corners of the rectangle with respect to the pointing position. Extended implementation of _reflected_regions() method to allow the placement of rectangular off regions, which is essentially finally placing a rectangle instead of a circle.
Andreas Specovius · Jürgen Knödlseder
1 parent a9168621
Showing 1 changed file with 121 additions and 3 deletions
cscripts/csphagen.py
@@ -363,6 +363,116 @@ class csphagen(ctools.csobservation):
         # Return
         return
  
+    def _compute_region_separation(self, pnt_dir):
+        """ Compute the separation angle for reflected off regions in radians.
+
+        Returns
+        -------
+        angle : `float`
+            Separation angle of two off regions (rad)
+        """
+        # Initialise the result
+        separation = -1
+
+        # Compute offset of reflected regions to pointing position
+        offset = pnt_dir.dist_deg(self._src_dir)
+
+        if self._srcshape == 'CIRCLE':
+
+            # Compute apparent diameter of the circle
+            separation = 2.0 * self._rad / offset
+
+        elif self._srcshape == 'RECT':
+
+            # Compute the opening angle towards combinations of rectangle corners.
+            # This method overestimates the real need of space between the
+            # rectangles, so the method may be optimised to gain more off regions!
+            # Anyway, it is assured that the off regions will never overlap.
+
+            # Get the corners
+            cs = [self._src_reg[0].get_corner(icorner) for icorner in range(4)]
+
+            # Auxillary, local method for opening angle computation
+            def __compute_posang(pnt_dir,a,b):
+                """ Compute the difference in position angle to the pointing in
+                    degrees.
+                """
+                posang_a = pnt_dir.posang_deg(a) %360
+                posang_b = pnt_dir.posang_deg(b) %360
+                posang = abs(posang_a - posang_b)
+                return posang
+
+            # Compute the 6 opening angles
+            combinations = [[0,1], [0,2], [0,3], [1,2], [1,3], [2,3]]
+            angles = [__compute_posang(pnt_dir, cs[i], cs[j]) for i,j in combinations]
+
+            # The desired separation is the maximum opening angle
+            separation = max(angles) *gammalib.deg2rad
+
+
+            ## Different approach:
+            ## Find the separation angle for which two rectangles to not overlap
+            ## Issue: Scanning not optimal.
+            # posang    = pnt_dir.posang_deg(self._src_dir)
+            # dphi_incr = 0.1
+            # dphi      = dphi_incr
+            # print("Src ctr direction: ra={} dec={}".format(self._src_dir.ra_deg(), self._src_dir.dec_deg()))
+            # print(self._src_reg)
+            # while dphi<360:
+
+            #     ctr_dir = pnt_dir.clone()
+            #     ctr_dir.rotate_deg(posang + dphi, offset)
+            #     print("Query ctr direction: ra={} dec={}".format(ctr_dir.ra_deg(), ctr_dir.dec_deg()))
+
+            #     # Create test rectangle
+            #     reg = gammalib.GSkyRegionRect(ctr_dir
+            #                                  ,self._reg_width
+            #                                  ,self._reg_height
+            #                                  ,self._reg_posang + dphi
+            #                                  )
+            #     print(reg)
+
+            #     # Check test rectangle corner containment in reference rectangle
+            #     is_corner_contained = self._src_reg.overlaps(reg)
+            #     # is_corner_contained = False
+            #     # for icorner in range(4):
+            #     #     corner_dir = reg.get_corner(icorner)
+            #     #     is_corner_contained |= self._src_reg.contains(corner_dir)
+
+            #     # Check if any corner was contained
+            #     if not is_corner_contained:
+            #         diameter = dphi
+            #         break
+
+            #     # Increment dphi
+            #     dphi += dphi_incr
+
+
+            ## Very different approach: Compute the width along the circular
+            ## region placement path.
+            ## Issue: Long regions will overlap in the inner part.
+            # def l(w,h,phi):
+            #     """ Super complicated method to compute the height along dec-axis
+            #         of a rotated rectangle.
+            #     """
+            #     phi = np.deg2rad(phi) % (2*np.pi)
+            #     phi_c = np.arctan2(w,h)
+            #     if phi>=(0.5*np.pi) and (phi<np.pi):
+            #         phi = np.pi-phi
+            #     elif (phi>=np.pi) and (phi<(1.5*np.pi)):
+            #         phi -=np.pi
+            #     elif (phi>=(1.5*np.pi)) and (phi<(2*np.pi)):
+            #         phi = 2*np.pi - phi
+
+            #     if (0<=phi) and (phi<phi_c):
+            #         result = h/2/np.cos(phi)
+            #     else:
+            #         result = w/2/np.sin(phi)
+            #     return 2*result
+
+        # Return
+        return separation
+
     def _reflected_regions(self, obs):
         """
         Calculate list of reflected regions for a single observation (pointing)
@@ -402,7 +512,7 @@ class csphagen(ctools.csobservation):
         # ... otherwise
         else:
             posang = pnt_dir.posang_deg(self._src_dir)
-            if self._srcshape == 'CIRCLE':
+            if (self._srcshape == 'CIRCLE') or (self._srcshape == 'RECT'):
  
                 # Determine number of background regions to skip
                 N_skip  = self['bkgregskip'].integer()
@@ -412,7 +522,7 @@ class csphagen(ctools.csobservation):
                 # Compute the angular separation of reflected regions wrt
                 # camera center. The factor 1.05 ensures background regions
                 # do not overlap due to numerical precision issues
-                alpha = 1.05 * 2.0 * self._rad / offset
+                alpha = 1.05 * self._compute_region_separation(pnt_dir)
  
                 # Compute number of reflected regions by dividing the angular
                 # separation by 2 pi.
@@ -437,7 +547,15 @@ class csphagen(ctools.csobservation):
                     while dphi <= dphi_max:
                         ctr_dir = pnt_dir.clone()
                         ctr_dir.rotate_deg(posang + dphi, offset)
-                        region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
+                        if self._srcshape == 'CIRCLE':
+                            region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
+                        elif self._srcshape == 'RECT':
+                            # Adjust the posang of the rectangle correspondingly
+                            region = gammalib.GSkyRegionRect(ctr_dir
+                                                            ,self._reg_width
+                                                            ,self._reg_height
+                                                            ,self._reg_posang + dphi
+                                                            )
                         if self._has_exclusion:
                             if self._excl_reg.overlaps(region):