GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSkyRegionCircle.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSkyRegionCircle.cpp - Circular sky region class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2016 by Michael Mayer *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GSkyRegionCircle.hpp
23  * @brief Circular sky region implementation
24  * @author Michael Mayer
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GSkyRegionCircle.hpp"
32 #include "GTools.hpp"
33 
34 /* __ Method name definitions ____________________________________________ */
35 #define G_RADIUS "GSkyRegionCircle::radius(double&)"
36 #define G_READ "GSkyRegionCircle::read(std::string&)"
37 #define G_CONTAINS "GSkyRegionCircle::contains(GSkyRegion&)"
38 #define G_OVERLAPS "GSkyRegionCircle::overlaps(GSkyRegion&)"
39 
40 /* __ Macros _____________________________________________________________ */
41 
42 /* __ Coding definitions _________________________________________________ */
43 
44 /* __ Debug definitions __________________________________________________ */
45 
46 /* __ Prototype __________________________________________________________ */
47 
48 
49 /*==========================================================================
50  = =
51  = Constructors/destructors =
52  = =
53  ==========================================================================*/
54 
55 /***********************************************************************//**
56  * @brief Void constructor
57  ***************************************************************************/
59 {
60  // Initialise members
61  init_members();
62 
63  // Return
64  return;
65 }
66 
67 
68 /***********************************************************************//**
69  * @brief Direction constructor
70  *
71  * @param[in] centre Centre sky direction.
72  * @param[in] radius Region radius [deg].
73  ***************************************************************************/
74 GSkyRegionCircle::GSkyRegionCircle(const GSkyDir& centre, const double& radius) :
75  GSkyRegion()
76 {
77  // Initialise members
78  init_members();
79 
80  // Set members
81  this->centre(centre);
82  this->radius(radius);
83 
84  // Compute solid angle
86 
87  // Return
88  return;
89 }
90 
91 
92 /***********************************************************************//**
93  * @brief Direction constructor
94  *
95  * @param[in] ra Right Ascension of region centre [deg].
96  * @param[in] dec Declination of region centre [deg].
97  * @param[in] radius Region radius [deg].
98  ***************************************************************************/
99 GSkyRegionCircle::GSkyRegionCircle(const double& ra, const double& dec,
100  const double& radius) : GSkyRegion()
101 {
102  // Initialise members
103  init_members();
104 
105  // Set members
106  this->centre(ra, dec);
107  this->radius(radius);
108 
109  // Compute solid angle
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief String constructor
119  *
120  * @param[in] line DS9 region file line.
121  *
122  * Constructs region from a DS9 region file line.
123  ***************************************************************************/
124 GSkyRegionCircle::GSkyRegionCircle(const std::string& line) : GSkyRegion()
125 {
126  // Initialise members
127  init_members();
128 
129  // Read information from DS9 region file line
130  read(line);
131 
132  // Return
133  return;
134 
135 }
136 
137 
138 /***********************************************************************//**
139  * @brief Copy constructor
140  *
141  * @param[in] region Circular sky region.
142  ***************************************************************************/
144  GSkyRegion(region)
145 {
146  // Initialise members
147  init_members();
148 
149  // Copy members
150  copy_members(region);
151 
152  // Return
153  return;
154 }
155 
156 
157 /***********************************************************************//**
158  * @brief Destructor
159  ***************************************************************************/
161 {
162  // Free members
163  free_members();
164 
165  // Return
166  return;
167 }
168 
169 
170 /*==========================================================================
171  = =
172  = Operators =
173  = =
174  ==========================================================================*/
175 
176 /***********************************************************************//**
177  * @brief Assignment operator
178  *
179  * @param[in] region Circular sky region.
180  * @return Circular sky region.
181  ***************************************************************************/
183 {
184  // Execute only if object is not identical
185  if (this != &region) {
186 
187  // Copy base class members
188  this->GSkyRegion::operator=(region);
189 
190  // Free members
191  free_members();
192 
193  // Initialise private members for clean destruction
194  init_members();
195 
196  // Copy members
197  copy_members(region);
198 
199  } // endif: object was not identical
200 
201  // Return this object
202  return *this;
203 }
204 
205 
206 /*==========================================================================
207  = =
208  = Public methods =
209  = =
210  ==========================================================================*/
211 
212 /***********************************************************************//**
213  * @brief Clear instance
214  ***************************************************************************/
216 {
217  // Free members
218  free_members();
219  this->GSkyRegion::free_members();
220 
221  // Initialise private members
222  this->GSkyRegion::init_members();
223  init_members();
224 
225  // Return
226  return;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Clone circular sky region
232  *
233  * @return Deep copy to circular sky region.
234  ***************************************************************************/
236 {
237  // Clone circular sky region
238  return new GSkyRegionCircle(*this);
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Set radius of circular region
244  *
245  * @param[in] radius Radius [deg].
246  *
247  * @exception GException::invalid_argument
248  * Radius value is less than 0.
249  *
250  * Sets the radius of the circular sky region. Only non-negative radii are
251  * allowed.
252  ***************************************************************************/
253 void GSkyRegionCircle::radius(const double& radius)
254 {
255  if (radius < 0.0) {
256  std::string msg =
257  "A radius of "+gammalib::str(radius)+" degrees has been"
258  " specified for a circular sky region.\n"
259  "The radius of a circular sky region can't be less than 0"
260  " degrees.";
262  }
263 
264  // Set radius value
265  m_radius = radius;
266 
267  // Compute the solid angle
269 
270  // Return
271  return;
272 }
273 
274 
275 /***********************************************************************//**
276  * @brief Read region from DS9 string
277  *
278  * @param[in] line String in DS9 format.
279  *
280  * @exception GException::invalid_value
281  * Invalid value found in DS9 format string.
282  ***************************************************************************/
283 void GSkyRegionCircle::read(const std::string& line)
284 {
285  // Clear the current instance
286  clear();
287 
288  // Split the string into 2 parts seperated by #
289  std::vector<std::string> substrings = gammalib::split(line,"#");
290  std::string region_def = (!substrings.empty()) ? substrings[0] : "";
291  std::string comment = (substrings.size() > 1) ? substrings[1] : "";
292 
293  // Finding the circle
294  if (region_def.find("circle") == std::string::npos) {
295  std::string msg =
296  "Unable to find the key word \"circle\" in provided string"
297  " \""+line+"\".\n"
298  "The \"circle\" key word is mandatory.";
299  throw GException::invalid_value(G_READ, msg);
300  }
301 
302  // Get the the coordinate system of the values
303  std::string system = gammalib::split(region_def, ";")[0];
304 
305  // Get the substring of the important values
306  unsigned pos = region_def.find("circle(");
307  unsigned end = region_def.find(")");
308  std::string circlestring = region_def.substr(pos+7, end);
309  circlestring.erase(circlestring.find(")"), 1);
310 
311  // Get the values of the region x,y,and radius
312  std::vector<std::string> values = gammalib::split(circlestring,",");
313  if (values.size() != 3) {
314  std::string msg =
315  "Invalid number of "+gammalib::str(values.size())+" arguments"
316  " after the \"circle\" key word in provided string \""+line+"\".\n"
317  "Exactly 3 arguments are expected.";
318  throw GException::invalid_value(G_READ, msg);
319  }
320  double x = gammalib::todouble(values[0]);
321  double y = gammalib::todouble(values[1]);
322  double radius = gammalib::todouble(values[2]);
323 
324  // Get radius units
325  if (gammalib::contains(values[2], "'")) {
326  radius /= 60.0;
327  }
328  else if (gammalib::contains(values[2], "\"")) {
329  radius /= 3600.0;
330  }
331 
332  // Initialise centre direction
333  GSkyDir centre = GSkyDir();
334 
335  // Set the values in the correct system
336  if (system == "fk5" || system == "icrs") {
337  centre.radec_deg(x,y);
338  }
339  else if (system == "galactic") {
340  centre.lb_deg(x,y);
341  }
342  else {
343  std::string msg =
344  "Unsupported coordinate system \""+system+"\" in provided string"
345  " \""+line+"\".\n"
346  "Only the following coordinate systems are supported: "
347  "\"fk5\", \"icrs\" and \"galactic\".";
348  throw GException::invalid_value(G_READ, msg);
349  }
350 
351  // Set members
352  this->centre(centre);
353  this->radius(radius);
354 
355  // Compute solid angle
357 
358  // Check if there is a given name for the region and set it
359  std::vector<std::string>comments = gammalib::split(comment, " ");
360  for (int i = 0; i < comments.size(); i++) {
361  if (gammalib::contains(comments[i], "text")) {
362  std::vector<std::string> attributes = gammalib::split(comments[i], "=");
363  if (attributes.size() < 2) {
364  std::string msg =
365  "Invalid character sequence encountered in provided"
366  " string \""+line+"\".\n"
367  "An attribute of the type \"text=Name\" is expected.";
368  throw GException::invalid_value(G_READ, msg);
369  }
370  m_name = attributes[1];
371  }
372  }
373 
374  // Return
375  return;
376 }
377 
378 
379 /***********************************************************************//**
380  * @brief Write region into a string
381  *
382  * @return String to be written in a DS9 region file
383  *
384  * Writes a DS9 region into a string. The region name is only written if it
385  * is defined.
386  ***************************************************************************/
387 std::string GSkyRegionCircle::write(void) const
388 {
389  // Allocate string
390  std::string result;
391 
392  // Set string
393  result.append("fk5;circle(");
394  result.append(gammalib::str(m_centre.ra_deg()));
395  result.append(",");
396  result.append(gammalib::str(m_centre.dec_deg()));
397  result.append(",");
398  result.append(gammalib::str(m_radius));
399  result.append(")");
400 
401  // Optionally add region name
402  if (m_name.length() > 0) {
403  result.append(" # text=");
404  result.append(m_name);
405  }
406 
407  // Return string
408  return result;
409 }
410 
411 
412 /***********************************************************************//**
413  * @brief Print circular region
414  *
415  * @param[in] chatter Chattiness
416  * @return String containing region information.
417  ***************************************************************************/
418 std::string GSkyRegionCircle::print(const GChatter& chatter) const
419 {
420  // Initialise result string
421  std::string result;;
422 
423  // Continue only if chatter is not silent
424  if (chatter != SILENT) {
425 
426  // Append header
427  result.append("=== GSkyRegionCircle ===");
428 
429  // Append sky circle information
430  result.append("\n"+gammalib::parformat("Right Ascension of centre"));
431  result.append(gammalib::str(m_centre.ra_deg())+" deg");
432  result.append("\n"+gammalib::parformat("Declination of centre"));
433  result.append(gammalib::str(m_centre.dec_deg())+" deg");
434  result.append("\n"+gammalib::parformat("Radius"));
435  result.append(gammalib::str(m_radius)+" deg");
436 
437  } // endif: chatter was not silent
438 
439  // Return result
440  return result;
441 }
442 
443 /***********************************************************************//**
444  * @brief Checks if sky direction lies within region
445  *
446  * @param[in] dir Sky direction.
447  *
448  * A sky direction lies within a region when its distance to the region
449  * centre is not larger than the region radius.
450  ***************************************************************************/
451 bool GSkyRegionCircle::contains(const GSkyDir& dir) const
452 {
453  // Initialise return value
454  bool dir_is_in = false;
455 
456  // calculate distance from dir to centre
457  double distance = dir.dist_deg(m_centre);
458 
459  // Check if distance < radius
460  if (distance <= radius()) {
461  dir_is_in = true;
462  }
463 
464  // Return bool
465  return dir_is_in;
466 }
467 
468 /***********************************************************************//**
469  * @brief Checks if region is fully contained within this region
470  *
471  * @param[in] reg Sky region.
472  *
473  * @exception GException::feature_not_implemented
474  * Regions differ in type.
475  ***************************************************************************/
477 {
478  // Initialise return value
479  bool fully_inside = false;
480 
481  // If other region is circle use a simple way to calculate
482  if (reg.type() == "Circle") {
483 
484  // Create circular region from reg
485  const GSkyRegionCircle* regcirc =
486  dynamic_cast<const GSkyRegionCircle*>(&reg);
487 
488  // Calculate angular distance between the centres
489  double ang_dist = m_centre.dist_deg(regcirc->centre());
490 
491  // Check if the region is contained in this
492  if ((ang_dist + regcirc->radius()) <= m_radius) {
493  fully_inside = true;
494  }
495 
496  }
497 
498  // ... otherwise throw an exception
499  else {
501  "Cannot compare two different region types yet");
502  }
503 
504  // Return value
505  return fully_inside;
506 }
507 
508 
509 /***********************************************************************//**
510  * @brief Checks if region is overlapping with this region
511  *
512  * @param[in] reg Sky region.
513  *
514  * @exception GException::feature_not_implemented
515  * Regions differ in type.
516  ***************************************************************************/
518 {
519  // Initialise return value
520  bool overlap = false;
521 
522  // If other region is circle use a simple way to calculate
523  if (reg.type() == "Circle") {
524 
525  // Create circular region from reg
526  const GSkyRegionCircle* regcirc =
527  dynamic_cast<const GSkyRegionCircle*>(&reg);
528 
529  // Calculate angular distance between the centres
530  double ang_dist = m_centre.dist_deg(regcirc->centre());
531 
532  // Check if the distance is smaller than the sum of both radii
533  if (ang_dist <= (m_radius + regcirc->radius())) {
534  overlap = true;
535  }
536 
537  }
538 
539  // ... otherwise throw an exception
540  else {
542  "Cannot compare two different region types yet");
543  }
544 
545  // Return value
546  return overlap;
547 }
548 
549 
550 
551 /*==========================================================================
552  = =
553  = Private methods =
554  = =
555  ==========================================================================*/
556 
557 /***********************************************************************//**
558  * @brief Initialise class members
559  ***************************************************************************/
561 {
562  // Initialise members
563  m_centre = GSkyDir();
564  m_radius = 0.0;
565  m_type = "Circle";
566 
567  //Return
568  return;
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Copy class members
574  *
575  * @param[in] region Circular sky region.
576  ***************************************************************************/
578 {
579  // Copy attributes
580  m_centre = region.m_centre;
581  m_radius = region.m_radius;
582 
583  // Return
584  return;
585 }
586 
587 
588 /***********************************************************************//**
589  * @brief Delete class members
590  ***************************************************************************/
592 {
593  // Return
594  return;
595 }
596 
597 
598 /***********************************************************************//**
599  * @brief Compute solid angle
600  ***************************************************************************/
602 {
603  // Convert radius from degrees to radians
604  double radius_rad = m_radius * gammalib::deg2rad;
605 
606  // Compute solid angle
607  m_solid = gammalib::twopi * (1.0 - std::cos(radius_rad));
608 
609  // Return
610  return;
611 }
void free_members(void)
Delete class members.
std::string print(const GChatter &chatter=NORMAL) const
Print circular region.
bool contains(const std::string &str, const std::string &substring)
Checks if a substring is in a string.
Definition: GTools.cpp:1221
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:280
double m_solid
Solid angle subtended by the region (sr)
Definition: GSkyRegion.hpp:96
bool overlaps(const GSkyRegion &reg) const
Checks if region is overlapping with this region.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
GSkyDir m_centre
Centre or reference point of the region.
void free_members(void)
Delete class members.
Definition: GSkyRegion.cpp:168
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
void init_members(void)
Initialise class members.
Definition: GSkyRegion.cpp:136
const std::string & type(void) const
Return region type.
Definition: GSkyRegion.hpp:137
GSkyRegionCircle(void)
Void constructor.
#define G_RADIUS
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:862
Gammalib tools definition.
void read(const std::string &line)
Read region from DS9 string.
Interface for the circular sky region class.
std::string m_type
Type of the region (circle, rectangle,...)
Definition: GSkyRegion.hpp:94
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:997
const double & radius(void) const
Return circular region radius (in degrees)
#define G_READ
virtual ~GSkyRegionCircle(void)
Destructor.
const double deg2rad
Definition: GMath.hpp:43
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
std::string write(void) const
Write region into a string.
std::string m_name
Name of the region.
Definition: GSkyRegion.hpp:95
double m_radius
Radius of circular the region [deg].
void compute_solid_angle(void)
Compute solid angle.
virtual GSkyRegion & operator=(const GSkyRegion &region)
Assignment operator.
Definition: GSkyRegion.cpp:106
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
GChatter
Definition: GTypemaps.hpp:33
#define G_OVERLAPS
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
#define G_CONTAINS
const GSkyDir & centre(void) const
Return circular region centre.
void copy_members(const GSkyRegionCircle &region)
Copy class members.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:256
bool contains(const GSkyDir &dir) const
Checks if sky direction lies within region.
const double twopi
Definition: GMath.hpp:36
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
GSkyRegionCircle & operator=(const GSkyRegionCircle &region)
Assignment operator.
Sky direction class.
Definition: GSkyDir.hpp:62
Circular sky region class interface definition.
GSkyRegionCircle * clone(void) const
Clone circular sky region.
void clear(void)
Clear instance.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:805
void init_members(void)
Initialise class members.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413