diff --git a/include/GSkyRegionMap.hpp b/include/GSkyRegionMap.hpp index b3a7779..63b3baf 100644 --- a/include/GSkyRegionMap.hpp +++ b/include/GSkyRegionMap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyRegionMap.hpp - Sky region map class * * ----------------------------------------------------------------------- * - * copyright (C) 2017-2018 by Pierrick Martin * + * copyright (C) 2017-2020 by Pierrick Martin * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * diff --git a/include/GSkyRegionRect.hpp b/include/GSkyRegionRect.hpp index e47d5b9..632d5b9 100644 --- a/include/GSkyRegionRect.hpp +++ b/include/GSkyRegionRect.hpp @@ -30,10 +30,12 @@ /* __ Includes ___________________________________________________________ */ #include #include -#include "GMath.hpp" +//#include "GMath.hpp" #include "GSkyDir.hpp" #include "GSkyRegion.hpp" -#include "GTools.hpp" +//#include "GTools.hpp" + +/* __ Forward declarations _______________________________________________ */ /***********************************************************************//** @@ -63,14 +65,15 @@ class GSkyRegionRect : public GSkyRegion { public: // Constructors and destructors GSkyRegionRect(void); - 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 - ); + GSkyRegionRect(const GSkyDir& centre, + const double& width, + const double& height, + const double& posang); + GSkyRegionRect(const double& ra, + const double& dec, + const double& width, + const double& height, + const double& posang); explicit GSkyRegionRect(const std::string& line); GSkyRegionRect(const GSkyRegionRect& region); virtual ~GSkyRegionRect(void); @@ -79,33 +82,30 @@ public: GSkyRegionRect& operator=(const GSkyRegionRect& region); // Implemented methods - 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; - void height(const double& height); - const GSkyDir& centre(void) const; - void centre(const GSkyDir& centre); - void centre(const double& ra,const double& dec); - double ra(void) const; - double dec(void) const; - void read(const std::string& line); - std::string write(void) const; - bool contains(const GSkyDir& dir) const; - bool contains_local(const GSkyDir& locdir) const; - bool contains(const GSkyRegion& reg) const; - bool overlaps(const GSkyRegion& reg) const; - std::string print(const GChatter& chatter = NORMAL) const; - - GSkyDir transform_to_local(const GSkyDir& skydir) const; - GSkyDir transform_to_global(const GSkyDir& locdir) const; - GSkyDir get_corner(const int& index) const; + void clear(void); + GSkyRegionRect* clone(void) const; + std::string classname(void) const; + const GSkyDir& centre(void) const; + void centre(const GSkyDir& centre); + void centre(const double& ra,const double& dec); + double ra(void) const; + double dec(void) const; + const double& width(void) const; + void width(const double& width); + const double& height(void) const; + void height(const double& height); + const double& posang(void) const; + void posang(const double& posang); + void read(const std::string& line); + std::string write(void) const; + bool contains(const GSkyDir& dir) const; + bool contains_local(const GSkyDir& locdir) const; + bool contains(const GSkyRegion& reg) const; + bool overlaps(const GSkyRegion& reg) const; + GSkyDir transform_to_local(const GSkyDir& skydir) const; + GSkyDir transform_to_global(const GSkyDir& locdir) const; + GSkyDir get_corner(const int& index) const; + std::string print(const GChatter& chatter = NORMAL) const; protected: // Protected methods @@ -116,9 +116,9 @@ protected: // Protected members 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] + double m_width; //!< Width of the region (degrees) + double m_height; //!< Height of the region (degrees) + double m_posang; //!< Position angle, counterclockwise from North (degrees) }; @@ -135,144 +135,136 @@ std::string GSkyRegionRect::classname(void) const /***********************************************************************//** - * @brief Return region position angle (in radians) + * @brief Return rectangular region centre * - * @return Region position angle [radians]. + * @return Region centre. * - * Returns the region position angle in radians. The position angle is - * counted counterclockwise from North. + * Returns the region centre. ***************************************************************************/ inline -double GSkyRegionRect::posang(void) const +const GSkyDir& GSkyRegionRect::centre(void) const { - // Return - return (m_posang); + return (m_centre); } /***********************************************************************//** - * @brief Return region position angle (in degrees) + * @brief Set rectangular region centre * - * @return Region position angle [degrees]. + * @param[in] dir Region centre. * - * Returns the region position angle in degrees. The position angle is - * counted counterclockwise from North. + * Sets the centre of the rectangular region to the specified sky direction. ***************************************************************************/ inline -double GSkyRegionRect::posang_deg(void) const +void GSkyRegionRect::centre(const GSkyDir& dir) { - // Return - return (m_posang * gammalib::rad2deg); + m_centre = dir; + return; } /***********************************************************************//** - * @brief Return region width extension (in degrees) + * @brief Set rectangular region centre Right Ascension and Declincation * - * @return Region width [deg]. + * @param[in] ra Right Ascension (degrees). + * @param[in] dec Declination (degrees). * - * Returns the region width extension in degrees. + * Sets the centre of the rectangular region to the specified Right Ascension + * and Declination. ***************************************************************************/ inline -double GSkyRegionRect::width(void) const +void GSkyRegionRect::centre(const double& ra, const double& dec) { - double width = 2*m_halfwidth; - return (width); + m_centre.radec_deg(ra, dec); + return; } /***********************************************************************//** - * @brief Return region height extension (in degrees) + * @brief Return rectangular region centre Right Ascension * - * @return Region height [deg]. + * @return Region centre Right Ascension (degrees). * - * Returns the region height extension in degrees. + * Returns the region centre Right Ascension in degrees. ***************************************************************************/ inline -double GSkyRegionRect::height(void) const +double GSkyRegionRect::ra(void) const { - double height = 2*m_halfheight; - return (height); + return (m_centre.ra_deg()); } /***********************************************************************//** - * @brief Return rectangular region centre + * @brief Return rectangular region centre Declination * - * @return Region centre. + * @return Region centre Declination (degrees). * - * Returns the region centre. + * Returns the region centre Declination in degrees. ***************************************************************************/ inline -const GSkyDir& GSkyRegionRect::centre(void) const +double GSkyRegionRect::dec(void) const { - return (m_centre); + return (m_centre.dec_deg()); } /***********************************************************************//** - * @brief Set rectangular region centre + * @brief Return region width extension (in degrees) * - * @param[in] dir Region centre. + * @return Region width (degrees). * - * Sets the centre of the rectangular region to the specified sky direction. + * Returns the region width extension in degrees. ***************************************************************************/ inline -void GSkyRegionRect::centre(const GSkyDir& dir) +const double& GSkyRegionRect::width(void) const { - // Set centre - m_centre = dir; - - // Return - return; + return (m_width); } /***********************************************************************//** - * @brief Set rectangular region centre Right Ascension and Declincation + * @brief Return region height extension (in degrees) * - * @param[in] ra Right Ascension [deg]. - * @param[in] dec Declination [deg]. + * @return Region height (degrees). * - * Sets the centre of the rectangular region to the specified Right Ascension - * and Declination. + * Returns the region height extension in degrees. ***************************************************************************/ inline -void GSkyRegionRect::centre(const double& ra, const double& dec) +const double& GSkyRegionRect::height(void) const { - // Set centre values - m_centre.radec_deg(ra,dec); - - // Return - return; + return (m_height); } /***********************************************************************//** - * @brief Return rectangular region centre Right Ascension + * @brief Return region position angle (in degrees) * - * @return Region centre Right Ascension [deg]. + * @return Region position angle (degrees). * - * Returns the region centre Right Ascension in degrees. + * Returns the region position angle in degrees. The position angle is + * counted counterclockwise from North. ***************************************************************************/ inline -double GSkyRegionRect::ra(void) const +const double& GSkyRegionRect::posang(void) const { - return (m_centre.ra_deg()); + return (m_posang); } /***********************************************************************//** - * @brief Return rectangular region centre Declination + * @brief Set position angle of rectangular region * - * @return Region centre Declination [deg]. + * @param[in] posang Position angle (degrees). * - * Returns the region centre Declination in degrees. + * Sets the position angle of the rectangular sky region. The position angle + * is counted counterclockwise from celestial North and is aligned to the + * rectangle height. ***************************************************************************/ inline -double GSkyRegionRect::dec(void) const +void GSkyRegionRect::posang(const double& posang) { - return (m_centre.dec_deg()); + m_posang = posang; + return; } #endif /* GSKYREGIONRECT_HPP */ diff --git a/pyext/GSkyRegionRect.i b/pyext/GSkyRegionRect.i index 9b298e1..4a549a6 100644 --- a/pyext/GSkyRegionRect.i +++ b/pyext/GSkyRegionRect.i @@ -35,48 +35,47 @@ * @brief Interface for the rectangular sky region class ***************************************************************************/ class GSkyRegionRect : public GSkyRegion { + public: // Constructors and destructors GSkyRegionRect(void); - 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 - ); + GSkyRegionRect(const GSkyDir& centre, + const double& width, + const double& height, + const double& posang); + GSkyRegionRect(const double& ra, + const double& dec, + const double& width, + const double& height, + const double& posang); explicit GSkyRegionRect(const std::string& line); GSkyRegionRect(const GSkyRegionRect& region); virtual ~GSkyRegionRect(void); // Implemented methods - 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; - void height(const double& height); - const GSkyDir& centre(void) const; - void centre(const GSkyDir& centre); - void centre(const double& ra,const double& dec); - double ra(void) const; - double dec(void) const; - void read(const std::string& line); - std::string write(void) const; - bool contains(const GSkyDir& dir) const; - bool contains_local(const GSkyDir& locdir) const; - bool contains(const GSkyRegion& reg) const; - bool overlaps(const GSkyRegion& reg) const; - - GSkyDir transform_to_local(const GSkyDir& skydir) const; - GSkyDir transform_to_global(const GSkyDir& localdir) const; - GSkyDir get_corner(const int& index) const; + void clear(void); + GSkyRegionRect* clone(void) const; + std::string classname(void) const; + const GSkyDir& centre(void) const; + void centre(const GSkyDir& centre); + void centre(const double& ra,const double& dec); + double ra(void) const; + double dec(void) const; + const double& width(void) const; + void width(const double& width); + const double& height(void) const; + void height(const double& height); + const double& posang(void) const; + void posang(const double& posang); + void read(const std::string& line); + std::string write(void) const; + bool contains(const GSkyDir& dir) const; + bool contains_local(const GSkyDir& locdir) const; + bool contains(const GSkyRegion& reg) const; + bool overlaps(const GSkyRegion& reg) const; + GSkyDir transform_to_local(const GSkyDir& skydir) const; + GSkyDir transform_to_global(const GSkyDir& locdir) const; + GSkyDir get_corner(const int& index) const; }; diff --git a/src/sky/GSkyRegionMap.cpp b/src/sky/GSkyRegionMap.cpp index d684f34..eddb256 100644 --- a/src/sky/GSkyRegionMap.cpp +++ b/src/sky/GSkyRegionMap.cpp @@ -657,24 +657,24 @@ void GSkyRegionMap::set_region_circle(const GSkyRegionCircle* circle) * @param[in] rect Region rectangle * * Sets a region map from a region rectangle. - * A skymap is allocated with a resolution goal of 0.01 degrees (same for both - * axes). - * The particular resolution of either axes depends on the w/h of the RoI of the - * rotated rectangle. - * The map extent is chosen dynamically to match the RoI of the rotated rectangle. - * Three sparse pixels are added at either side of the map. - * The sky map is created in celestial coordinates and TAN projection. + * + * A skymap is allocated with a resolution goal of 0.01 degrees (same for + * both axes). The particular resolution of either axes depends on the width + * over height of the RoI of the rotated rectangle. The map extent is chosen + * dynamically to match the RoI of the rotated rectangle. Three sparse pixels + * are added at either side of the map. The sky map is created in celestial + * coordinates and TAN projection. ***************************************************************************/ void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect) { // Define desired resolution (in degrees) - const double resolution_goal = 0.01; + const double resolution_goal = 0.01; // Try to estimate extent of the (rotated) rectangle: // Loop over the rectangle corners to get the extension - double dxmax=0; - double dymax=0; - for(int icorner=0; icorner<4; ++icorner) { + double dxmax = 0.0; + double dymax = 0.0; + for (int icorner = 0; icorner < 4; ++icorner) { // Get skydir of current corner GSkyDir corner = rect->get_corner(icorner); @@ -691,22 +691,19 @@ void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect) if (dy > dymax) { dymax = dy; } - if (dx > dxmax) { dxmax = dx; } - } + + } // endfor: looped over all corners // Compute width of map RoI in degrees - double map_width = 2*dxmax *gammalib::rad2deg; - double map_height = 2*dymax *gammalib::rad2deg; - // std::cout<<"Generate map for rectangle:\n"<<(*rect)<ra(), rect->dec(), -resolution_x, resolution_y, - nbins_x, nbins_y - ); + nbins_x, nbins_y); - // Set sky map pixels + // Set sky map pixels that are contained in rectangle to 1 for (int ipix = 0; ipix < map_fine.npix(); ++ipix) { - - // If pixel is contained in rectangle then set map pixel to 1 if (rect->contains(map_fine.inx2dir(ipix))) { map_fine(ipix) = 1.0; } - - } // endfor: looped over sky map pixels + } // Set sky map as region map this->map(map_fine); // Return return; -} \ No newline at end of file +} diff --git a/src/sky/GSkyRegionRect.cpp b/src/sky/GSkyRegionRect.cpp index 23fa64b..6f34591 100644 --- a/src/sky/GSkyRegionRect.cpp +++ b/src/sky/GSkyRegionRect.cpp @@ -28,18 +28,19 @@ #ifdef HAVE_CONFIG_H #include #endif +//#include "GMath.hpp" +#include "GTools.hpp" #include "GSkyRegionCircle.hpp" #include "GSkyRegionRect.hpp" #include "GSkyRegionMap.hpp" -#include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_WIDTH "GSkyRegionRect::width(double&)" -#define G_HEIGHT "GSkyRegionRect::height(double&)" -#define G_READ "GSkyRegionRect::read(std::string&)" -#define G_CONTAINS "GSkyRegionRect::contains(GSkyRegion&)" -#define G_OVERLAPS "GSkyRegionRect::overlaps(GSkyRegion&)" -#define G_GET_CORNER "GSkyRegionRect::get_corner(int&)" +#define G_WIDTH "GSkyRegionRect::width(double&)" +#define G_HEIGHT "GSkyRegionRect::height(double&)" +#define G_READ "GSkyRegionRect::read(std::string&)" +#define G_CONTAINS "GSkyRegionRect::contains(GSkyRegion&)" +#define G_OVERLAPS "GSkyRegionRect::overlaps(GSkyRegion&)" +#define G_GET_CORNER "GSkyRegionRect::get_corner(int&)" /* __ Macros _____________________________________________________________ */ @@ -70,24 +71,26 @@ GSkyRegionRect::GSkyRegionRect(void) : GSkyRegion() /***********************************************************************//** - * @brief Direction constructor + * @brief Sky direction constructor * - * @param[in] centre Centre sky direction. - * @param[in] w Region width [deg]. - * @param[in] h Region height [deg]. + * @param[in] centre Centre of rectangle. + * @param[in] width Region width (degrees). + * @param[in] height Region height (degrees). + * @param[in] posang Position angle (degrees). ***************************************************************************/ -GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w, - const double& h, const double& posang_deg) : - GSkyRegion() +GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, + const double& width, + const double& height, + const double& posang) : GSkyRegion() { // Initialise members init_members(); // Set members this->centre(centre); - this->width(w); - this->height(h); - this->posang_deg(posang_deg); + this->width(width); + this->height(height); + this->posang(posang); // Compute solid angle compute_solid_angle(); @@ -98,26 +101,28 @@ GSkyRegionRect::GSkyRegionRect(const GSkyDir& centre, const double& w, /***********************************************************************//** - * @brief Direction constructor + * @brief Right Ascension and Declination constructor * * @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]. + * @param[in] width Region width (degrees). + * @param[in] height Region height (degrees). + * @param[in] posang Position angle (degrees). ***************************************************************************/ -GSkyRegionRect::GSkyRegionRect(const double& ra, const double& dec, - const double& w, const double& h, - const double& posang_deg) : - GSkyRegion() +GSkyRegionRect::GSkyRegionRect(const double& ra, + const double& dec, + const double& width, + const double& height, + const double& posang) : GSkyRegion() { // Initialise members init_members(); // Set members this->centre(ra, dec); - this->width(w); - this->height(h); - this->posang_deg(posang_deg); + this->width(width); + this->height(height); + this->posang(posang); // Compute solid angle compute_solid_angle(); @@ -153,8 +158,7 @@ GSkyRegionRect::GSkyRegionRect(const std::string& line) : GSkyRegion() * * @param[in] region Rectangular sky region. ***************************************************************************/ -GSkyRegionRect::GSkyRegionRect(const GSkyRegionRect& region) : - GSkyRegion(region) +GSkyRegionRect::GSkyRegionRect(const GSkyRegionRect& region) : GSkyRegion(region) { // Initialise members init_members(); @@ -253,67 +257,29 @@ GSkyRegionRect* GSkyRegionRect::clone(void) const /***********************************************************************//** - * @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 cel. North and is aligned to the rect - * height. - ***************************************************************************/ -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 cel. North and is aligned to the rect - * height. - ***************************************************************************/ -void GSkyRegionRect::posang_deg(const double& posang) -{ - // Set the position angle - this->posang(posang * gammalib::deg2rad); - - // Return - return; -} - - -/***********************************************************************//** * @brief Set width of rectangular region * - * @param[in] width Width [deg]. + * @param[in] width Width (degrees). * * @exception GException::invalid_argument * Width value is less than 0. * - * Sets the width of the rectangular sky region. Only non-negative radii are - * allowed. + * Sets the @p width of the rectangular sky region. Only non-negative widths + * are allowed. ***************************************************************************/ void GSkyRegionRect::width(const double& width) { + // Throw an exception if the width is less than zero if (width < 0.0) { std::string msg = - "A width of "+gammalib::str(width)+" degrees has been" - " specified for a rectangular sky region.\n" - "The width of a rectangular sky region can't be less than 0" - " degrees."; + "A negative width of "+gammalib::str(width)+" degrees has been " + "specified for a rectangular sky region. Please specify a " + "non-negative width."; throw GException::invalid_argument(G_WIDTH, msg); } // Set radius value - m_halfwidth = 0.5*width; + m_width = width; // Compute the solid angle compute_solid_angle(); @@ -326,27 +292,29 @@ void GSkyRegionRect::width(const double& width) /***********************************************************************//** * @brief Set height of rectangular region * - * @param[in] height Height [deg]. + * @param[in] height Height (degrees). * * @exception GException::invalid_argument * Height value is less than 0. * - * Sets the height of the rectangular sky region. Only non-negative radii are - * allowed. For zero posang, the height axis is pointing to celestial north. + * Sets the @p height of the rectangular sky region. Only non-negative + * heights are allowed. + * + * Note that for a position angle of zero, the height axis is pointing to + * celestial north. ***************************************************************************/ void GSkyRegionRect::height(const double& height) { if (height < 0.0) { std::string msg = - "A height of "+gammalib::str(height)+" degrees has been" - " specified for a rectangular sky region.\n" - "The height of a rectangular sky region can't be less than 0" - " degrees."; + "A negative height of "+gammalib::str(height)+" degrees has been " + "specified for a rectangular sky region. Please specify a " + "non-negative height."; throw GException::invalid_argument(G_HEIGHT, msg); } // Set radius value - m_halfheight = 0.5*height; + m_height = height; // Compute the solid angle compute_solid_angle(); @@ -377,9 +345,8 @@ void GSkyRegionRect::read(const std::string& line) // Finding the box if (region_def.find("box") == std::string::npos) { std::string msg = - "Unable to find the key word \"box\" in provided string" - " \""+line+"\".\n" - "The \"box\" key word is mandatory."; + "Unable to find the key word \"box\" in provided string " + "\""+line+"\". The \"box\" key word is mandatory."; throw GException::invalid_value(G_READ, msg); } @@ -396,29 +363,29 @@ void GSkyRegionRect::read(const std::string& line) std::vector values = gammalib::split(regionstring,","); if (values.size() != 5) { std::string msg = - "Invalid number of "+gammalib::str(values.size())+" arguments" - " after the \"box\" key word in provided string \""+line+"\".\n" + "Invalid number of "+gammalib::str(values.size())+" arguments " + "after the \"box\" key word in provided string \""+line+"\". " "Exactly 5 arguments are expected."; throw GException::invalid_value(G_READ, msg); } double x = gammalib::todouble(values[0]); double y = gammalib::todouble(values[1]); - double w = gammalib::todouble(values[2]); - double h = gammalib::todouble(values[3]); + double width = gammalib::todouble(values[2]); + double height = gammalib::todouble(values[3]); double posang = gammalib::todouble(values[4]); // Get radius units if (gammalib::contains(values[2], "'")) { - w /= 60.0; + width /= 60.0; } else if (gammalib::contains(values[2], "\"")) { - w /= 3600.0; + width /= 3600.0; } if (gammalib::contains(values[3], "'")) { - h /= 60.0; + height /= 60.0; } else if (gammalib::contains(values[3], "\"")) { - h /= 3600.0; + height /= 3600.0; } if (gammalib::contains(values[4], "'")) { posang /= 60.0; @@ -439,32 +406,31 @@ void GSkyRegionRect::read(const std::string& line) } else { std::string msg = - "Unsupported coordinate system \""+system+"\" in provided string" - " \""+line+"\".\n" - "Only the following coordinate systems are supported: " - "\"fk5\", \"icrs\" and \"galactic\"."; + "Unsupported coordinate system \""+system+"\" in provided string " + "\""+line+"\". Only the following coordinate systems are " + "supported: \"fk5\", \"icrs\" and \"galactic\"."; throw GException::invalid_value(G_READ, msg); } // Set members this->centre(centre); - this->width(w); - this->height(h); - this->posang_deg(posang); + this->width(width); + this->height(height); + this->posang(posang); // Compute solid angle compute_solid_angle(); // Check if there is a given name for the region and set it std::vectorcomments = gammalib::split(comment, " "); - for (int i = 0; i < comments.size(); i++) { + for (int i = 0; i < comments.size(); ++i) { if (gammalib::contains(comments[i], "text")) { std::vector attributes = gammalib::split(comments[i], "="); if (attributes.size() < 2) { std::string msg = - "Invalid character sequence encountered in provided" - " string \""+line+"\".\n" - "An attribute of the type \"text=Name\" is expected."; + "Invalid character sequence encountered in provided " + "string \""+line+"\". An attribute of the type " + "\"text=Name\" is expected."; throw GException::invalid_value(G_READ, msg); } m_name = attributes[1]; @@ -495,11 +461,11 @@ std::string GSkyRegionRect::write(void) const result.append(","); result.append(gammalib::str(m_centre.dec_deg())); result.append(","); - result.append(gammalib::str(2*m_halfwidth)); + result.append(gammalib::str(m_width)); result.append(","); - result.append(gammalib::str(2*m_halfheight)); + result.append(gammalib::str(m_height)); result.append(","); - result.append(gammalib::str(m_posang*gammalib::rad2deg)); + result.append(gammalib::str(m_posang)); result.append(")"); // Optionally add region name @@ -514,41 +480,6 @@ std::string GSkyRegionRect::write(void) const /***********************************************************************//** - * @brief Print rectangular region - * - * @param[in] chatter Chattiness - * @return String containing region information. - ***************************************************************************/ -std::string GSkyRegionRect::print(const GChatter& chatter) const -{ - // Initialise result string - std::string result;; - - // Continue only if chatter is not silent - if (chatter != SILENT) { - - // Append header - result.append("=== GSkyRegionRect ==="); - - // Append sky region information - result.append("\n"+gammalib::parformat("Right Ascension of centre")); - result.append(gammalib::str(m_centre.ra_deg())+" deg"); - result.append("\n"+gammalib::parformat("Declination of centre")); - result.append(gammalib::str(m_centre.dec_deg())+" deg"); - result.append("\n"+gammalib::parformat("Width")); - 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 - - // Return result - return result; -} - -/***********************************************************************//** * @brief Checks if sky direction lies within region * * @param[in] dir Sky direction. @@ -572,7 +503,7 @@ bool GSkyRegionRect::contains(const GSkyDir& dir) const /***********************************************************************//** * @brief Checks if local direction lies within region * - * @param[in] dir Direction in local coordinate system. + * @param[in] locdir Sky direction in local coordinate system. * * A local direction lies within a region when its distance to the region * centre is not larger than the region extension in both axes directions. @@ -583,10 +514,10 @@ bool GSkyRegionRect::contains_local(const GSkyDir& locdir) const bool dir_is_in = false; // Check containment: declination axis - if ((std::abs(locdir.dec_deg()) - 1.0e-10) <= m_halfheight) { + if ((std::abs(locdir.dec_deg()) - 1.0e-10) <= 0.5 * m_height) { // Check containment: right ascension axis - if ((std::abs(locdir.ra_deg()) - 1.0e-10) <= m_halfwidth) { + if ((std::abs(locdir.ra_deg()) - 1.0e-10) <= 0.5 * m_width) { // Sky direction is inside the rectangle dir_is_in = true; @@ -625,10 +556,10 @@ bool GSkyRegionRect::contains(const GSkyRegion& reg) const const double circ_rad = regcirc->radius(); // Check extension of circle along declination axis - if ((std::abs(local_centre.dec_deg()) + circ_rad) <= m_halfheight) { + if ((std::abs(local_centre.dec_deg()) + circ_rad) <= 0.5 * m_height) { // Check extension of circle along right ascension axis - if ((std::abs(local_centre.ra_deg()) + circ_rad) <= m_halfwidth) { + if ((std::abs(local_centre.ra_deg()) + circ_rad) <= 0.5 * m_width) { // Circle is fully contained in this rectangle is_fully_inside = true; @@ -644,7 +575,7 @@ bool GSkyRegionRect::contains(const GSkyRegion& reg) const dynamic_cast(®); // Loop over the four corners of the rectangle :regrect: - for(int icorner=0; icorner<4; ++icorner) { + for (int icorner = 0; icorner < 4; ++icorner) { // Get skydir of current corner GSkyDir corner = regrect->get_corner(icorner); @@ -654,7 +585,7 @@ bool GSkyRegionRect::contains(const GSkyRegion& reg) const break; } // Did we arrive at the last corner? - else if (icorner==3) { + else if (icorner == 3) { // As we reached this point, the three prior corners were inside // and the last corner now is also inside: @@ -705,10 +636,10 @@ bool GSkyRegionRect::overlaps(const GSkyRegion& reg) const const double circ_rad = regcirc->radius(); // Check extension of circle along declination axis - if (std::abs(local_centre.dec_deg()) <= (m_halfheight + circ_rad)) { + if (std::abs(local_centre.dec_deg()) <= (0.5 * m_height + circ_rad)) { // Check extension of circle along right ascension axis - if (std::abs(local_centre.ra_deg()) <= (m_halfwidth + circ_rad)) { + if (std::abs(local_centre.ra_deg()) <= (0.5 * m_width + circ_rad)) { // Circle is fully contained in this rectangle is_overlapping = true; @@ -722,7 +653,7 @@ bool GSkyRegionRect::overlaps(const GSkyRegion& reg) const const GSkyRegionRect* regrect = dynamic_cast(®); - // Dirty cludge: compare vs map + // Dirty kludge: compare vs map GSkyRegionMap regmap = GSkyRegionMap(regrect); is_overlapping = regmap.overlaps(*this); @@ -746,17 +677,17 @@ bool GSkyRegionRect::overlaps(const GSkyRegion& reg) const * @param[in] skydir Sky direction in global coordinate system. * @return Sky direction in local cartesian region object coordinates. * - * Transform the sky direction :skydir: to the local cartesian coordinate system. - * The origin of the local coordinate system is fixed to the center of the - * rectangle and aligned in +ra (width) and +dec (height). - * Please note, that a GSkyDir object is returned, even if the local coordsys + * Transform the sky direction @p skydir to the local cartesian coordinate + * system. The origin of the local coordinate system is fixed to the center + * of the rectangle and aligned in +ra (width) and +dec (height). Please + * note, that a GSkyDir object is returned, even if the local coordsys is * is cartesian. ***************************************************************************/ GSkyDir GSkyRegionRect::transform_to_local(const GSkyDir& skydir) const { // Get distance and polar angle relative to rectangle centre double dist = m_centre.dist(skydir); - double pa = m_centre.posang(skydir) - m_posang; + double pa = m_centre.posang(skydir) - m_posang * gammalib::deg2rad; // Compute corresponding x and y coordinate from polar coords double x = dist * std::sin(pa); @@ -777,7 +708,7 @@ GSkyDir GSkyRegionRect::transform_to_local(const GSkyDir& skydir) const * @param[in] locdir Direction in local coordinate system. * @return Sky direction in global coordinate system. * - * Transform the local direction :locdir: to the global coordinate system. + * Transform the local direction @p locdir to the global coordinate system. * The origin of the local coordinate system is fixed to the center of the * rectangle and aligned in +ra (width) and +dec (height). ***************************************************************************/ @@ -793,11 +724,8 @@ GSkyDir GSkyRegionRect::transform_to_global(const GSkyDir& locdir) const // Compute dist and position angle in polar coords double dist = std::sqrt(x*x + y*y); - double pa = gammalib::pihalf - std::atan2(y, x) + m_posang; - - // // Compute distance and position angle for locdir in global coordinates - // double dist = centre_local.dist(locdir); - // double pa = centre_local.posang(locdir) + m_posang; + double pa = gammalib::pihalf - std::atan2(y, x) + + m_posang * gammalib::deg2rad; // Create output global sky direction GSkyDir transformed(m_centre); @@ -822,14 +750,14 @@ GSkyDir GSkyRegionRect::transform_to_global(const GSkyDir& locdir) const GSkyDir GSkyRegionRect::get_corner(const int& index) const { // Assert index is in [0,3] - if ((index<0) || (index>3)) { - throw GException::out_of_range(G_GET_CORNER, index, 0, 3); + if ((index < 0) || (index > 3)) { + throw GException::out_of_range(G_GET_CORNER, "Corner index", index, 3); } // Compute the offset to the corners // Index counts CCW starting from +RA,+DEC: +R+D, +R-D, -R-D, -R+D - double dx = m_halfwidth * (1 - 2*int(index >= 2)); - double dy = m_halfheight * (1 - 2*int((index==1) || (index==2))); + double dx = 0.5 * m_width * (1 - 2*int(index >= 2)); + double dy = 0.5 * m_height * (1 - 2*int((index==1) || (index==2))); // Define local corner sky direction GSkyDir corner; @@ -843,6 +771,41 @@ GSkyDir GSkyRegionRect::get_corner(const int& index) const } +/***********************************************************************//** + * @brief Print rectangular region + * + * @param[in] chatter Chattiness + * @return String containing region information. + ***************************************************************************/ +std::string GSkyRegionRect::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result;; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GSkyRegionRect ==="); + + // Append sky region information + result.append("\n"+gammalib::parformat("Right Ascension of centre")); + result.append(gammalib::str(ra())+" deg"); + result.append("\n"+gammalib::parformat("Declination of centre")); + result.append(gammalib::str(dec())+" deg"); + result.append("\n"+gammalib::parformat("Width")); + result.append(gammalib::str(width())+" deg"); + result.append("\n"+gammalib::parformat("Height")); + result.append(gammalib::str(height())+" deg"); + result.append("\n"+gammalib::parformat("PA")); + result.append(gammalib::str(posang())+" deg"); + + } // endif: chatter was not silent + + // Return result + return result; +} + /*========================================================================== = = @@ -856,11 +819,11 @@ GSkyDir GSkyRegionRect::get_corner(const int& index) const void GSkyRegionRect::init_members(void) { // Initialise members - m_centre = GSkyDir(); - m_halfwidth = 0.0; - m_halfheight = 0.0; - m_posang = 0.0; m_type = "Rect"; + m_centre.clear(); + m_width = 0.0; + m_height = 0.0; + m_posang = 0.0; //Return return; @@ -875,10 +838,10 @@ void GSkyRegionRect::init_members(void) void GSkyRegionRect::copy_members(const GSkyRegionRect& region) { // Copy attributes - m_centre = region.m_centre; - m_halfwidth = region.m_halfwidth; - m_halfheight = region.m_halfheight; - m_posang = region.m_posang; + m_centre = region.m_centre; + m_width = region.m_width; + m_height = region.m_height; + m_posang = region.m_posang; // Return return; @@ -896,13 +859,12 @@ void GSkyRegionRect::free_members(void) /***********************************************************************//** - * @brief Compute solid angle [sr] + * @brief Compute solid angle (sr) ***************************************************************************/ void GSkyRegionRect::compute_solid_angle(void) { // Compute solid angle - m_solid = (2*m_halfwidth) * (2*m_halfheight) * \ - (gammalib::deg2rad * gammalib::deg2rad); + m_solid = m_width * m_height * gammalib::deg2rad * gammalib::deg2rad; // Return return;