From d371f910fbac7433f1d9eb17b7ebb3ae6322e20d Mon Sep 17 00:00:00 2001 From: Andreas Specovius Date: Mon, 15 Jun 2020 14:25:28 +0200 Subject: [PATCH] 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. --- include/GSkyRegionRect.hpp | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- pyext/GSkyRegionRect.i | 24 +++++++++++++++++------- src/sky/GSkyRegionRect.cpp | 45 +++++++++++++++++++++++---------------------- 3 files changed, 124 insertions(+), 31 deletions(-) diff --git a/include/GSkyRegionRect.hpp b/include/GSkyRegionRect.hpp index 8016e77..d57eb2c 100644 --- a/include/GSkyRegionRect.hpp +++ b/include/GSkyRegionRect.hpp @@ -33,6 +33,7 @@ #include "GMath.hpp" #include "GSkyDir.hpp" #include "GSkyRegion.hpp" +#include "GTools.hpp" /***********************************************************************//** @@ -53,8 +54,14 @@ class GSkyRegionRect : public GSkyRegion { public: // Constructors and destructors GSkyRegionRect(void); - GSkyRegionRect(const GSkyDir& centre, const double& w, const double& h); - GSkyRegionRect(const double& ra, const double& dec, const double& w, const double& h); + GSkyRegionRect(const GSkyDir& centre + ,const double& w, const double& h + ,const double& posang_deg + ); + GSkyRegionRect(const double& ra, const double& dec + ,const double& w, const double& h + ,const double& posang_deg + ); explicit GSkyRegionRect(const std::string& line); GSkyRegionRect(const GSkyRegionRect& region); virtual ~GSkyRegionRect(void); @@ -66,6 +73,10 @@ public: void clear(void); GSkyRegionRect* clone(void) const; std::string classname(void) const; + double posang(void) const; + void posang(const double& posang); + double posang_deg(void) const; + void posang_deg(const double& posang); double width(void) const; void width(const double& width); double height(void) const; @@ -93,6 +104,7 @@ protected: GSkyDir m_centre; //!< Centre or reference point of the region double m_halfwidth; //!< Half width of the region [deg] double m_halfheight; //!< Half height of the region [deg] + double m_posang; //!< Position angle, counterclockwise from North [radians] }; @@ -109,6 +121,76 @@ std::string GSkyRegionRect::classname(void) const /***********************************************************************//** + * @brief Return region position angle (in radians) + * + * @return Region position angle [radians]. + * + * Returns the region position angle in radians. The position angle is + * counted counterclockwise from North. + ***************************************************************************/ +inline +double GSkyRegionRect::posang(void) const +{ + // Return + return (m_posang); +} + + +/***********************************************************************//** + * @brief Return region position angle (in degrees) + * + * @return Region position angle [degrees]. + * + * Returns the region position angle in degrees. The position angle is + * counted counterclockwise from North. + ***************************************************************************/ +inline +double GSkyRegionRect::posang_deg(void) const +{ + // Return + return (m_posang * gammalib::rad2deg); +} + + +/***********************************************************************//** + * @brief Set position angle of rectangular region + * + * @param[in] posang Position angle [radians]. + * + * Sets the position angle of the rectangular sky region. The position angle + * is counted counterclockwise from North. + ***************************************************************************/ +inline +void GSkyRegionRect::posang(const double& posang) +{ + // Set the position angle + m_posang = posang; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set position angle of rectangular region + * + * @param[in] posang Position angle [degrees]. + * + * Sets the position angle of the rectangular sky region. The position angle + * is counted counterclockwise from North. + ***************************************************************************/ +inline +void GSkyRegionRect::posang_deg(const double& posang) +{ + // Set the position angle + m_posang = posang * gammalib::deg2rad; + + // Return + return; +} + + +/***********************************************************************//** * @brief Return region width extension (in degrees) * * @return Region width [deg]. diff --git a/pyext/GSkyRegionRect.i b/pyext/GSkyRegionRect.i index 9c93d7b..d30633f 100644 --- a/pyext/GSkyRegionRect.i +++ b/pyext/GSkyRegionRect.i @@ -1,7 +1,7 @@ /*************************************************************************** - * GSkyRegionRect.i - Sky region class * + * GSkyRegionRect.i - Sky region class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2019 by Andreas Specovius * + * copyright (C) 2019-2020 by Andreas Specovius * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -38,8 +38,14 @@ class GSkyRegionRect : public GSkyRegion { public: // Constructors and destructors GSkyRegionRect(void); - GSkyRegionRect(const GSkyDir& centre, const double& w, const double& h); - GSkyRegionRect(const double& ra, const double& dec, const double& w, const double& h); + GSkyRegionRect(const GSkyDir& centre + ,const double& w, const double& h + ,const double& posang_deg + ); + GSkyRegionRect(const double& ra, const double& dec + ,const double& w, const double& h + ,const double& posang_deg + ); explicit GSkyRegionRect(const std::string& line); GSkyRegionRect(const GSkyRegionRect& region); virtual ~GSkyRegionRect(void); @@ -48,6 +54,10 @@ public: void clear(void); GSkyRegionRect* clone(void) const; std::string classname(void) const; + double posang(void) const; + void posang(const double& posang); + double posang_deg(void) const; + void posang_deg(const double& posang); double width(void) const; void width(const double& width); double height(void) const; @@ -74,10 +84,10 @@ public: } %pythoncode { def __getstate__(self): - state = (self.centre(), self.width(), self.height(), self.name()) + state = (self.centre(), self.width(), self.height(), self.posang(), self.name()) return state def __setstate__(self, state): - self.__init__(state[0], state[1], state[2]) - self.name(state[3]) + self.__init__(state[0], state[1], state[2], state[3]) + self.name(state[4]) } }; diff --git a/src/sky/GSkyRegionRect.cpp b/src/sky/GSkyRegionRect.cpp index f3b2765..f0dbfd8 100644 --- a/src/sky/GSkyRegionRect.cpp +++ b/src/sky/GSkyRegionRect.cpp @@ -74,7 +74,7 @@ GSkyRegionRect::GSkyRegionRect(void) : GSkyRegion() * @param[in] h Region height [deg]. ***************************************************************************/ GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w, - const double& h) : + const double& h, const double& posang_deg) : GSkyRegion() { // Initialise members @@ -84,6 +84,7 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w, this->centre(centre); this->width(w); this->height(h); + this->posang_deg(posang_deg); // Compute solid angle compute_solid_angle(); @@ -99,9 +100,11 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w, * @param[in] ra Right Ascension of region centre [deg]. * @param[in] dec Declination of region centre [deg]. * @param[in] radius Region radius [deg]. + * @param[in] radius Region radius [deg]. ***************************************************************************/ GSkyRegionRect::GSkyRegionRect(const double& ra, const double& dec, - const double& w, const double& h) : + const double& w, const double& h, + const double& posang_deg) : GSkyRegion() { // Initialise members @@ -111,6 +114,7 @@ GSkyRegionRect::GSkyRegionRect(const double& ra, const double& dec, this->centre(ra, dec); this->width(w); this->height(h); + this->posang_deg(posang_deg); // Compute solid angle compute_solid_angle(); @@ -344,7 +348,7 @@ void GSkyRegionRect::read(const std::string& line) // Get the substring of the important values unsigned pos = region_def.find("box("); unsigned end = region_def.find(")"); - std::string regionstring = region_def.substr(pos+7, end); + std::string regionstring = region_def.substr(pos+4, end); regionstring.erase(regionstring.find(")"), 1); // Get the values of the region x,y,and radius @@ -362,14 +366,6 @@ void GSkyRegionRect::read(const std::string& line) double h = gammalib::todouble(values[3]); double posang = gammalib::todouble(values[4]); - // Assert posang is zero - if (std::abs(posang)>0.1) { - std::string msg = - "Box position angle is not zero! Handling of rotated rectangular" - " sky regions is generally possible but not implemented."; - throw GException::invalid_value(G_READ, msg); - } - // Get radius units if (gammalib::contains(values[2], "'")) { w /= 60.0; @@ -383,6 +379,12 @@ void GSkyRegionRect::read(const std::string& line) else if (gammalib::contains(values[3], "\"")) { h /= 3600.0; } + if (gammalib::contains(values[4], "'")) { + posang /= 60.0; + } + else if (gammalib::contains(values[4], "\"")) { + posang /= 3600.0; + } // Initialise centre direction GSkyDir centre = GSkyDir(); @@ -407,6 +409,7 @@ void GSkyRegionRect::read(const std::string& line) this->centre(centre); this->width(w); this->height(h); + this->posang_deg(posang); // Compute solid angle compute_solid_angle(); @@ -455,7 +458,7 @@ std::string GSkyRegionRect::write(void) const result.append(","); result.append(gammalib::str(2*m_halfheight)); result.append(","); - result.append(gammalib::str(0)); + result.append(gammalib::str(m_posang*gammalib::rad2deg)); result.append(")"); // Optionally add region name @@ -495,6 +498,8 @@ std::string GSkyRegionRect::print(const GChatter& chatter) const result.append(gammalib::str(2*m_halfwidth)+" deg"); result.append("\n"+gammalib::parformat("Height")); result.append(gammalib::str(2*m_halfheight)+" deg"); + result.append("\n"+gammalib::parformat("PA")); + result.append(gammalib::str(m_posang*gammalib::rad2deg)+" deg"); } // endif: chatter was not silent @@ -597,6 +602,7 @@ void GSkyRegionRect::init_members(void) m_centre = GSkyDir(); m_halfwidth = 0.0; m_halfheight = 0.0; + m_posang = 0.0; m_type = "Rect"; //Return @@ -612,9 +618,10 @@ void GSkyRegionRect::init_members(void) void GSkyRegionRect::copy_members(const GSkyRegionRect& region) { // Copy attributes - m_centre = region.m_centre; + m_centre = region.m_centre; m_halfwidth = region.m_halfwidth; m_halfheight = region.m_halfheight; + m_posang = region.m_posang; // Return return; @@ -632,19 +639,13 @@ void GSkyRegionRect::free_members(void) /***********************************************************************//** - * @brief Compute solid angle - * - * As width and height are opening angles the rectangular boundaries are - * not aligned to spherical coordinate lines, the box keeps it opening angles. - * Hence the solidangle will vary. + * @brief Compute solid angle [sr] ***************************************************************************/ void GSkyRegionRect::compute_solid_angle(void) { - // Compute width scaling factor from box center - double scale = std::cos(m_centre.dec()); - // Compute solid angle - m_solid = (2*m_halfwidth)/scale * (2*m_halfheight); + m_solid = (2*m_halfwidth) * (2*m_halfheight) * \ + (gammalib::deg2rad * gammalib::deg2rad); // Return return; -- libgit2 0.21.4