Eirik Hatlen / gammalib

Commit d371f910fbac7433f1d9eb17b7ebb3ae6322e20d

Authored by Andreas Specovius
Committed by Jürgen Knödlseder
1 parent 4638c52b

Added position angle :m_posang:, stored in radians. Can be read from string line… 

… (ds9 line) and written to ds9 string. Can be set via radians and degrees, the latter is implemented for the constructor. Changed solid angle computation algorithm.
Inline Side-by-side
Showing 3 changed files with 124 additions and 31 deletions
include/GSkyRegionRect.hpp
  View file @d371f91
... ... @@ -33,6 +33,7 @@
33 33 #include "GMath.hpp"
34 34 #include "GSkyDir.hpp"
35 35 #include "GSkyRegion.hpp"
  36 +#include "GTools.hpp"
36 37  
37 38  
38 39 /***********************************************************************//**
... ... @@ -53,8 +54,14 @@ class GSkyRegionRect : public GSkyRegion {
53 54 public:
54 55 // Constructors and destructors
55 56 GSkyRegionRect(void);
56   - GSkyRegionRect(const GSkyDir& centre, const double& w, const double& h);
57   - GSkyRegionRect(const double& ra, const double& dec, const double& w, const double& h);
  57 + GSkyRegionRect(const GSkyDir& centre
  58 + ,const double& w, const double& h
  59 + ,const double& posang_deg
  60 + );
  61 + GSkyRegionRect(const double& ra, const double& dec
  62 + ,const double& w, const double& h
  63 + ,const double& posang_deg
  64 + );
58 65 explicit GSkyRegionRect(const std::string& line);
59 66 GSkyRegionRect(const GSkyRegionRect& region);
60 67 virtual ~GSkyRegionRect(void);
... ... @@ -66,6 +73,10 @@ public:
66 73 void clear(void);
67 74 GSkyRegionRect* clone(void) const;
68 75 std::string classname(void) const;
  76 + double posang(void) const;
  77 + void posang(const double& posang);
  78 + double posang_deg(void) const;
  79 + void posang_deg(const double& posang);
69 80 double width(void) const;
70 81 void width(const double& width);
71 82 double height(void) const;
... ... @@ -93,6 +104,7 @@ protected:
93 104 GSkyDir m_centre; //!< Centre or reference point of the region
94 105 double m_halfwidth; //!< Half width of the region [deg]
95 106 double m_halfheight; //!< Half height of the region [deg]
  107 + double m_posang; //!< Position angle, counterclockwise from North [radians]
96 108 };
97 109  
98 110  
... ... @@ -109,6 +121,76 @@ std::string GSkyRegionRect::classname(void) const
109 121  
110 122  
111 123 /***********************************************************************//**
  124 + * @brief Return region position angle (in radians)
  125 + *
  126 + * @return Region position angle [radians].
  127 + *
  128 + * Returns the region position angle in radians. The position angle is
  129 + * counted counterclockwise from North.
  130 + ***************************************************************************/
  131 +inline
  132 +double GSkyRegionRect::posang(void) const
  133 +{
  134 + // Return
  135 + return (m_posang);
  136 +}
  137 +
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Return region position angle (in degrees)
  141 + *
  142 + * @return Region position angle [degrees].
  143 + *
  144 + * Returns the region position angle in degrees. The position angle is
  145 + * counted counterclockwise from North.
  146 + ***************************************************************************/
  147 +inline
  148 +double GSkyRegionRect::posang_deg(void) const
  149 +{
  150 + // Return
  151 + return (m_posang * gammalib::rad2deg);
  152 +}
  153 +
  154 +
  155 +/***********************************************************************//**
  156 + * @brief Set position angle of rectangular region
  157 + *
  158 + * @param[in] posang Position angle [radians].
  159 + *
  160 + * Sets the position angle of the rectangular sky region. The position angle
  161 + * is counted counterclockwise from North.
  162 + ***************************************************************************/
  163 +inline
  164 +void GSkyRegionRect::posang(const double& posang)
  165 +{
  166 + // Set the position angle
  167 + m_posang = posang;
  168 +
  169 + // Return
  170 + return;
  171 +}
  172 +
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Set position angle of rectangular region
  176 + *
  177 + * @param[in] posang Position angle [degrees].
  178 + *
  179 + * Sets the position angle of the rectangular sky region. The position angle
  180 + * is counted counterclockwise from North.
  181 + ***************************************************************************/
  182 +inline
  183 +void GSkyRegionRect::posang_deg(const double& posang)
  184 +{
  185 + // Set the position angle
  186 + m_posang = posang * gammalib::deg2rad;
  187 +
  188 + // Return
  189 + return;
  190 +}
  191 +
  192 +
  193 +/***********************************************************************//**
112 194 * @brief Return region width extension (in degrees)
113 195 *
114 196 * @return Region width [deg].
... ...
pyext/GSkyRegionRect.i
  View file @d371f91
1 1 /***************************************************************************
2   - * GSkyRegionRect.i - Sky region class *
  2 + * GSkyRegionRect.i - Sky region class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2019 by Andreas Specovius *
  4 + * copyright (C) 2019-2020 by Andreas Specovius *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -38,8 +38,14 @@ class GSkyRegionRect : public GSkyRegion {
38 38 public:
39 39 // Constructors and destructors
40 40 GSkyRegionRect(void);
41   - GSkyRegionRect(const GSkyDir& centre, const double& w, const double& h);
42   - GSkyRegionRect(const double& ra, const double& dec, const double& w, const double& h);
  41 + GSkyRegionRect(const GSkyDir& centre
  42 + ,const double& w, const double& h
  43 + ,const double& posang_deg
  44 + );
  45 + GSkyRegionRect(const double& ra, const double& dec
  46 + ,const double& w, const double& h
  47 + ,const double& posang_deg
  48 + );
43 49 explicit GSkyRegionRect(const std::string& line);
44 50 GSkyRegionRect(const GSkyRegionRect& region);
45 51 virtual ~GSkyRegionRect(void);
... ... @@ -48,6 +54,10 @@ public:
48 54 void clear(void);
49 55 GSkyRegionRect* clone(void) const;
50 56 std::string classname(void) const;
  57 + double posang(void) const;
  58 + void posang(const double& posang);
  59 + double posang_deg(void) const;
  60 + void posang_deg(const double& posang);
51 61 double width(void) const;
52 62 void width(const double& width);
53 63 double height(void) const;
... ... @@ -74,10 +84,10 @@ public:
74 84 }
75 85 %pythoncode {
76 86 def __getstate__(self):
77   - state = (self.centre(), self.width(), self.height(), self.name())
  87 + state = (self.centre(), self.width(), self.height(), self.posang(), self.name())
78 88 return state
79 89 def __setstate__(self, state):
80   - self.__init__(state[0], state[1], state[2])
81   - self.name(state[3])
  90 + self.__init__(state[0], state[1], state[2], state[3])
  91 + self.name(state[4])
82 92 }
83 93 };
... ...
src/sky/GSkyRegionRect.cpp
  View file @d371f91
... ... @@ -74,7 +74,7 @@ GSkyRegionRect::GSkyRegionRect(void) : GSkyRegion()
74 74 * @param[in] h Region height [deg].
75 75 ***************************************************************************/
76 76 GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w,
77   - const double& h) :
  77 + const double& h, const double& posang_deg) :
78 78 GSkyRegion()
79 79 {
80 80 // Initialise members
... ... @@ -84,6 +84,7 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir&amp; centre, const double&amp; w,
84 84 this->centre(centre);
85 85 this->width(w);
86 86 this->height(h);
  87 + this->posang_deg(posang_deg);
87 88  
88 89 // Compute solid angle
89 90 compute_solid_angle();
... ... @@ -99,9 +100,11 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir&amp; centre, const double&amp; w,
99 100 * @param[in] ra Right Ascension of region centre [deg].
100 101 * @param[in] dec Declination of region centre [deg].
101 102 * @param[in] radius Region radius [deg].
  103 + * @param[in] radius Region radius [deg].
102 104 ***************************************************************************/
103 105 GSkyRegionRect::GSkyRegionRect(const double& ra, const double& dec,
104   - const double& w, const double& h) :
  106 + const double& w, const double& h,
  107 + const double& posang_deg) :
105 108 GSkyRegion()
106 109 {
107 110 // Initialise members
... ... @@ -111,6 +114,7 @@ GSkyRegionRect::GSkyRegionRect(const double&amp; ra, const double&amp; dec,
111 114 this->centre(ra, dec);
112 115 this->width(w);
113 116 this->height(h);
  117 + this->posang_deg(posang_deg);
114 118  
115 119 // Compute solid angle
116 120 compute_solid_angle();
... ... @@ -344,7 +348,7 @@ void GSkyRegionRect::read(const std::string&amp; line)
344 348 // Get the substring of the important values
345 349 unsigned pos = region_def.find("box(");
346 350 unsigned end = region_def.find(")");
347   - std::string regionstring = region_def.substr(pos+7, end);
  351 + std::string regionstring = region_def.substr(pos+4, end);
348 352 regionstring.erase(regionstring.find(")"), 1);
349 353  
350 354 // Get the values of the region x,y,and radius
... ... @@ -362,14 +366,6 @@ void GSkyRegionRect::read(const std::string&amp; line)
362 366 double h = gammalib::todouble(values[3]);
363 367 double posang = gammalib::todouble(values[4]);
364 368  
365   - // Assert posang is zero
366   - if (std::abs(posang)>0.1) {
367   - std::string msg =
368   - "Box position angle is not zero! Handling of rotated rectangular"
369   - " sky regions is generally possible but not implemented.";
370   - throw GException::invalid_value(G_READ, msg);
371   - }
372   -
373 369 // Get radius units
374 370 if (gammalib::contains(values[2], "'")) {
375 371 w /= 60.0;
... ... @@ -383,6 +379,12 @@ void GSkyRegionRect::read(const std::string&amp; line)
383 379 else if (gammalib::contains(values[3], "\"")) {
384 380 h /= 3600.0;
385 381 }
  382 + if (gammalib::contains(values[4], "'")) {
  383 + posang /= 60.0;
  384 + }
  385 + else if (gammalib::contains(values[4], "\"")) {
  386 + posang /= 3600.0;
  387 + }
386 388  
387 389 // Initialise centre direction
388 390 GSkyDir centre = GSkyDir();
... ... @@ -407,6 +409,7 @@ void GSkyRegionRect::read(const std::string&amp; line)
407 409 this->centre(centre);
408 410 this->width(w);
409 411 this->height(h);
  412 + this->posang_deg(posang);
410 413  
411 414 // Compute solid angle
412 415 compute_solid_angle();
... ... @@ -455,7 +458,7 @@ std::string GSkyRegionRect::write(void) const
455 458 result.append(",");
456 459 result.append(gammalib::str(2*m_halfheight));
457 460 result.append(",");
458   - result.append(gammalib::str(0));
  461 + result.append(gammalib::str(m_posang*gammalib::rad2deg));
459 462 result.append(")");
460 463  
461 464 // Optionally add region name
... ... @@ -495,6 +498,8 @@ std::string GSkyRegionRect::print(const GChatter&amp; chatter) const
495 498 result.append(gammalib::str(2*m_halfwidth)+" deg");
496 499 result.append("\n"+gammalib::parformat("Height"));
497 500 result.append(gammalib::str(2*m_halfheight)+" deg");
  501 + result.append("\n"+gammalib::parformat("PA"));
  502 + result.append(gammalib::str(m_posang*gammalib::rad2deg)+" deg");
498 503  
499 504 } // endif: chatter was not silent
500 505  
... ... @@ -597,6 +602,7 @@ void GSkyRegionRect::init_members(void)
597 602 m_centre = GSkyDir();
598 603 m_halfwidth = 0.0;
599 604 m_halfheight = 0.0;
  605 + m_posang = 0.0;
600 606 m_type = "Rect";
601 607  
602 608 //Return
... ... @@ -612,9 +618,10 @@ void GSkyRegionRect::init_members(void)
612 618 void GSkyRegionRect::copy_members(const GSkyRegionRect& region)
613 619 {
614 620 // Copy attributes
615   - m_centre = region.m_centre;
  621 + m_centre = region.m_centre;
616 622 m_halfwidth = region.m_halfwidth;
617 623 m_halfheight = region.m_halfheight;
  624 + m_posang = region.m_posang;
618 625  
619 626 // Return
620 627 return;
... ... @@ -632,19 +639,13 @@ void GSkyRegionRect::free_members(void)
632 639  
633 640  
634 641 /***********************************************************************//**
635   - * @brief Compute solid angle
636   - *
637   - * As width and height are opening angles the rectangular boundaries are
638   - * not aligned to spherical coordinate lines, the box keeps it opening angles.
639   - * Hence the solidangle will vary.
  642 + * @brief Compute solid angle [sr]
640 643 ***************************************************************************/
641 644 void GSkyRegionRect::compute_solid_angle(void)
642 645 {
643   - // Compute width scaling factor from box center
644   - double scale = std::cos(m_centre.dec());
645   -
646 646 // Compute solid angle
647   - m_solid = (2*m_halfwidth)/scale * (2*m_halfheight);
  647 + m_solid = (2*m_halfwidth) * (2*m_halfheight) * \
  648 + (gammalib::deg2rad * gammalib::deg2rad);
648 649  
649 650 // Return
650 651 return;
... ...