From 20133af489caec0efc4be818c501c97a8fa678a1 Mon Sep 17 00:00:00 2001 From: Andreas Specovius Date: Tue, 16 Jun 2020 18:31:40 +0200 Subject: [PATCH] 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. --- cscripts/csphagen.py | 124 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 121 insertions(+), 3 deletions(-) diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py index 282fd5a..2b7ee40 100755 --- a/cscripts/csphagen.py +++ b/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) 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