GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPointing.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPointing.cpp - CTA pointing class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2019 by Juergen Knoedlseder *
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 GCTAPointing.cpp
23  * @brief CTA pointing class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GFilename.hpp"
33 #include "GMatrix.hpp"
34 #include "GSkyDir.hpp"
35 #include "GCTAPointing.hpp"
36 #include "GCTAInstDir.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_READ_XML "GCTAPointing::read(GXmlElement&)"
40 #define G_WRITE_XML "GCTAPointing::write(GXmlElement&)"
41 
42 /* __ Macros _____________________________________________________________ */
43 
44 /* __ Coding definitions _________________________________________________ */
45 //#define G_USE_VECTORS
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  ***************************************************************************/
61 {
62  // Initialise members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief Sky direction constructor
72  *
73  * @param[in] dir Sky direction.
74  *
75  * Construct CTA pointing from sky direction.
76  ***************************************************************************/
78 {
79  // Initialise members
80  init_members();
81 
82  // Assign sky direction
83  this->dir(dir);
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief XML constructor
92  *
93  * @param[in] xml XML element.
94  *
95  * Construct CTA pointing from XML element.
96  ***************************************************************************/
98 {
99  // Initialise members
100  init_members();
101 
102  // Read information from XML element
103  read(xml);
104 
105  // Return
106  return;
107 }
108 
109 
110 /***********************************************************************//**
111  * @brief Copy constructor
112  *
113  * @param[in] pnt CTA pointing.
114  ***************************************************************************/
116 {
117  // Initialise members
118  init_members();
119 
120  // Copy members
121  copy_members(pnt);
122 
123  // Return
124  return;
125 }
126 
127 
128 /***********************************************************************//**
129  * @brief Destructor
130  ***************************************************************************/
132 {
133  // Free members
134  free_members();
135 
136  // Return
137  return;
138 }
139 
140 
141 /*==========================================================================
142  = =
143  = Operators =
144  = =
145  ==========================================================================*/
146 
147 /***********************************************************************//**
148  * @brief Assignment operator
149  *
150  * @param[in] pnt CTA pointing.
151  * @return CTA pointing.
152  ***************************************************************************/
154 {
155  // Execute only if object is not identical
156  if (this != &pnt) {
157 
158  // Free members
159  free_members();
160 
161  // Initialise members
162  init_members();
163 
164  // Copy members
165  copy_members(pnt);
166 
167  } // endif: object was not identical
168 
169  // Return this object
170  return *this;
171 }
172 
173 
174 /*==========================================================================
175  = =
176  = Public methods =
177  = =
178  ==========================================================================*/
179 
180 /***********************************************************************//**
181  * @brief Clear CTA pointing
182  ***************************************************************************/
184 {
185  // Free members
186  free_members();
187 
188  // Initialise private members
189  init_members();
190 
191  // Return
192  return;
193 }
194 
195 
196 /***********************************************************************//**
197  * @brief Clone CTA pointing
198  *
199  * @return Poiter to deep copy of CTA pointing.
200  ***************************************************************************/
202 {
203  return new GCTAPointing(*this);
204 }
205 
206 
207 /***********************************************************************//**
208  * @brief Set pointing direction
209  *
210  * @param[in] dir Sky direction of pointing.
211  *
212  * Set the pointing direction to the specified @p sky direction.
213  ***************************************************************************/
214 void GCTAPointing::dir(const GSkyDir& dir)
215 {
216  // Set sky direction
217  m_dir = dir;
218 
219  // Invalidate cache
220  m_has_cache = false;
221 
222  // Signal that pointing is valid
223  m_valid = true;
224 
225  // Return
226  return;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Get instrument direction from sky direction
232  *
233  * @param[in] skydir Sky direction.
234  * @return CTA instrument direction.
235  *
236  * Returns instrument direction as function of a sky direction
237  ***************************************************************************/
239 {
240  #if defined(G_USE_VECTORS)
241  // Compute rotation matrix
242  update();
243 
244  // Get celestial vector from sky coordinate
245  GVector celvector = skydir.celvector();
246 
247  // Transform to instrument system
248  GVector inst = m_Rback.transpose() * celvector;
249 
250  // Initialise instrument coordinates
251  double detx(0.0);
252  double dety(0.0);
253 
254  // Get offset from cel vector, i.e. distance from camera center
255  double theta = std::acos(inst[2]);
256 
257  // Check if theta and phi are defined
258  if (theta > 0.0 ) {
259  double phi = std::asin(inst[1] / std::sin(theta));
260  detx = theta * std::cos(phi);
261  dety = theta * std::sin(phi);
262  }
263  #else
264  double theta = m_dir.dist(skydir);
265  double phi = m_dir.posang(skydir);
266  double detx(0.0);
267  double dety(0.0);
268  if (theta > 0.0 ) {
269  detx = theta * std::cos(phi);
270  dety = theta * std::sin(phi);
271  }
272  #endif
273 
274  // Initialise instrument direction
275  GCTAInstDir inst_direction(skydir);
276  inst_direction.detx(detx);
277  inst_direction.dety(dety);
278 
279  // Return
280  return inst_direction;
281 }
282 
283 
284 /***********************************************************************//**
285  * @brief Get sky direction direction from instrument direction
286  *
287  * @param[in] instdir Instrument direction.
288  * @return Sky direction.
289  *
290  * Returns sky direction as function of an instrument direction
291  ***************************************************************************/
293 {
294  // Compute rotation matrix
295  update();
296 
297  // Retrieve instrument coordinates
298  double inst_x = instdir.detx();
299  double inst_y = instdir.dety();
300 
301  // Convert to polar coordinates
302  double theta = std::sqrt(inst_x * inst_x + inst_y * inst_y);
303  double phi = std::atan2(inst_y, inst_x);
304  double sin_phi = std::sin(phi);
305  double cos_phi = std::cos(phi);
306  double sin_theta = std::sin(theta);
307  double cos_theta = std::cos(theta);
308 
309  // Build vector from polar coordinates
310  GVector native(-cos_phi*sin_theta, sin_phi*sin_theta, cos_theta);
311 
312  // Rotate from instrument system into sky system
313  GVector skyvector = m_Rback * native;
314  GSkyDir sky;
315  sky.celvector(skyvector);
316 
317  // Return
318  return sky;
319 
320 }
321 
322 
323 /***********************************************************************//**
324  * @brief Return rotation matrix
325  *
326  * @return Rotation matrix.
327  ***************************************************************************/
328 const GMatrix& GCTAPointing::rot(void) const
329 {
330  // Update cache
331  update();
332 
333  // Return rotation matrix
334  return m_Rback;
335 }
336 
337 
338 /***********************************************************************//**
339  * @brief Read pointing from XML element
340  *
341  * @param[in] xml XML element.
342  *
343  * @exception GException::invalid_value
344  * Invalid XML format encountered.
345  *
346  * Read pointing parameter from an XML element. The format of the pointing
347  * parameter is
348  *
349  * <parameter name="Pointing" ra="83.0" dec="22.1"/>
350  *
351  ***************************************************************************/
353 {
354  // Clear pointing
355  clear();
356 
357  // Get pointing parameter
358  const GXmlElement* par = gammalib::xml_get_par(G_READ_XML, xml, "Pointing");
359 
360  // Extract position attributes
361  if (par->has_attribute("ra") && par->has_attribute("dec")) {
362  double ra = gammalib::todouble(par->attribute("ra"));
363  double dec = gammalib::todouble(par->attribute("dec"));
364  m_dir.radec_deg(ra,dec);
365  }
366  else if (par->has_attribute("lon") && par->has_attribute("lat")) {
367  double lon = gammalib::todouble(par->attribute("lon"));
368  double lat = gammalib::todouble(par->attribute("lat"));
369  m_dir.lb_deg(lon,lat);
370  }
371  else {
372  std::string msg = "Attributes \"ra\" and \"dec\" or \"lon\" and"
373  " \"lat\"not found in XML parameter \"Pointing\"."
374  " Please verify the XML format.";
376  }
377 
378  // Signal that pointing is valid
379  m_valid = true;
380 
381  // Return
382  return;
383 }
384 
385 
386 /***********************************************************************//**
387  * @brief Write pointing information into XML element
388  *
389  * @param[in] xml XML element.
390  *
391  * Writes pointing parameter into an XML element. The format of the pointing
392  * parameter is
393  *
394  * <parameter name="Pointing" ra="83.0" dec="22.1"/>
395  *
396  ***************************************************************************/
398 {
399  // Get parameter
400  GXmlElement* par = gammalib::xml_need_par(G_WRITE_XML, xml, "Pointing");
401 
402  // Write attributes
403  par->attribute("ra", gammalib::str(m_dir.ra_deg()));
404  par->attribute("dec", gammalib::str(m_dir.dec_deg()));
405 
406  // Return
407  return;
408 }
409 
410 
411 /***********************************************************************//**
412  * @brief Print CTA pointing information
413  *
414  * @param[in] chatter Chattiness.
415  * @return String containing pointing information.
416  ***************************************************************************/
417 std::string GCTAPointing::print(const GChatter& chatter) const
418 {
419  // Initialise result string
420  std::string result;
421 
422  // Continue only if chatter is not silent
423  if (chatter != SILENT) {
424 
425  // Append header
426  result.append("=== GCTAPointing ===");
427 
428  // Append information
429  result.append("\n"+gammalib::parformat("Pointing direction"));
430  result.append(this->dir().print());
431 
432  } // endif: chatter was not silent
433 
434  // Return result
435  return result;
436 }
437 
438 
439 /*==========================================================================
440  = =
441  = Private methods =
442  = =
443  ==========================================================================*/
444 
445 /***********************************************************************//**
446  * @brief Initialise class members
447  ***************************************************************************/
449 {
450  // Initialise members
451  m_dir.clear();
452  m_valid = false;
453  m_zenith = 0.0;
454  m_azimuth = 0.0;
455 
456  // Initialise cache
457  m_has_cache = false;
458  m_Rback.clear();
459 
460  // Return
461  return;
462 }
463 
464 
465 /***********************************************************************//**
466  * @brief Copy class members
467  *
468  * @param[in] pnt CTA pointing.
469  ***************************************************************************/
471 {
472  // Copy members
473  m_dir = pnt.m_dir;
474  m_valid = pnt.m_valid;
475  m_zenith = pnt.m_zenith;
476  m_azimuth = pnt.m_azimuth;
477 
478  // Copy cache
479  m_has_cache = pnt.m_has_cache;
480  m_Rback = pnt.m_Rback;
481 
482  // Return
483  return;
484 }
485 
486 
487 /***********************************************************************//**
488  * @brief Delete class members
489  ***************************************************************************/
491 {
492  // Return
493  return;
494 }
495 
496 
497 /***********************************************************************//**
498  * @brief Update coordinate transformation cache
499  ***************************************************************************/
500 void GCTAPointing::update(void) const
501 {
502  // Update coordinate transformation cache only if required
503  if (!m_has_cache) {
504 
505  // Set up Euler matrices
506  GMatrix Ry;
507  GMatrix Rz;
508  Ry.eulery(m_dir.dec_deg() - 90.0);
509  Rz.eulerz(-m_dir.ra_deg());
510 
511  // Compute rotation matrix
512  m_Rback = (Ry * Rz).transpose();
513 
514  // Signal that we have a valid transformation cache
515  m_has_cache = true;
516 
517  } // endif: Update of cache was required
518 
519  // Return
520  return;
521 }
void clear(void)
Clear CTA pointing.
GMatrix m_Rback
Rotation matrix.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
void detx(const double &x)
Set DETX coordinate (in radians)
GSkyDir skydir(const GCTAInstDir &instdir) const
Get sky direction direction from instrument direction.
Sky direction class interface definition.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
Generic matrix class definition.
XML element node class.
Definition: GXmlElement.hpp:47
CTA pointing class interface definition.
virtual ~GCTAPointing(void)
Destructor.
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Gammalib tools definition.
void free_members(void)
Delete class members.
void read(const GXmlElement &xml)
Read pointing from XML element.
void write(GXmlElement &xml) const
Write pointing information into XML element.
#define G_READ_XML
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
virtual void clear(void)
Clear matrix.
Definition: GMatrix.cpp:511
bool m_has_cache
Has transformation cache.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
double m_zenith
Pointing zenith angle (deg)
const GMatrix & rot(void) const
Return rotation matrix.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GMatrix transpose(void) const
Return transposed matrix.
Definition: GMatrix.cpp:847
double m_azimuth
Pointing azimuth angle (deg)
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
void celvector(const GVector &vector)
Set sky direction from 3D vector in celestial coordinates.
Definition: GSkyDir.cpp:291
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1511
CTA pointing class.
GCTAPointing * clone(void) const
Clone CTA pointing.
CTA instrument direction class interface definition.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
GChatter
Definition: GTypemaps.hpp:33
virtual GCTAPointing & operator=(const GCTAPointing &pnt)
Assignment operator.
#define G_WRITE_XML
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:142
bool m_valid
Validity flag.
std::string print(const GChatter &chatter=NORMAL) const
Print CTA pointing information.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:256
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
void copy_members(const GCTAPointing &pnt)
Copy class members.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
const GSkyDir & dir(void) const
Return pointing sky direction.
void dety(const double &y)
Set DETY coordinate (in radians)
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1020
double posang(const GSkyDir &dir) const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1022
Sky direction class.
Definition: GSkyDir.hpp:62
void update(void) const
Update coordinate transformation cache.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1560
Filename class interface definition.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:803
GSkyDir m_dir
Pointing direction in sky coordinates.
GCTAPointing(void)
Void constructor.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:411
void init_members(void)
Initialise class members.