GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GHorizDir.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GHorizDir.cpp - Horizontal direction class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014 by Karl Kosack *
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 GHorizDir.cpp
23  * @brief Horizontal direction class implementation
24  * @author Karl Kosack
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GHorizDir.hpp"
36 #include "GMatrix.hpp"
37 #include "GVector.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 
47 /* __ Prototypes _________________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Constructor
58  ***************************************************************************/
60 {
61  // Initialise members
62  init_members();
63 
64  // Return
65  return;
66 }
67 
68 
69 /***********************************************************************//**
70  * @brief Copy constructor
71  *
72  * @param[in] dir Horizontal direction.
73  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Copy members
80  copy_members(dir);
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief Destructor
89  ***************************************************************************/
91 {
92  // Free members
93  free_members();
94 
95  // Return
96  return;
97 }
98 
99 
100 /*==========================================================================
101  = =
102  = Operators =
103  = =
104  ==========================================================================*/
105 
106 /***********************************************************************//**
107  * @brief Assignment operator
108  *
109  * @param[in] dir Horizontal direction.
110  * @return Horizontal direction.
111  ***************************************************************************/
113 {
114  // Execute only if object is not identical
115  if (this != &dir) {
116 
117  // Free members
118  free_members();
119 
120  // Initialise members
121  init_members();
122 
123  // Copy members
124  copy_members(dir);
125 
126  } // endif: object was not identical
127 
128  // Return this object
129  return *this;
130 }
131 
132 
133 /*==========================================================================
134  = =
135  = Public methods =
136  = =
137  ==========================================================================*/
138 
139 /***********************************************************************//**
140  * @brief Clear horizontal direction
141  ***************************************************************************/
143 {
144  // Free members
145  free_members();
146 
147  // Initialise members
148  init_members();
149 
150  // Return
151  return;
152 }
153 
154 
155 /***********************************************************************//**
156  * @brief Clone horizontal direction
157  *
158  * @return Pointer to deep copy of horizontal direction.
159  ***************************************************************************/
161 {
162  // Clone horizon direction
163  return new GHorizDir(*this);
164 }
165 
166 
167 /***********************************************************************//**
168  * @brief Set horizontal direction (radians)
169  *
170  * @param[in] alt Altitude in radians.
171  * @param[in] az Azimuth in radians.
172  *
173  * Sets Altitude and Azimuth in radians.
174  ***************************************************************************/
175 void GHorizDir::altaz(const double& alt, const double& az)
176 {
177 
178  // Set direction
179  m_alt = alt;
180  m_az = az;
181 
182  // Return
183  return;
184 }
185 
186 
187 /***********************************************************************//**
188  * @brief Set horizontal direction (degrees)
189  *
190  * @param[in] alt Altitude in radians.
191  * @param[in] az Azimuth in radians.
192  *
193  * Sets Altitude and Azimuth in degrees.
194  ***************************************************************************/
195 void GHorizDir::altaz_deg(const double& alt, const double& az)
196 {
197 
198  // Set direction
199  m_az = az * gammalib::deg2rad;
200  m_alt = alt * gammalib::deg2rad;
201 
202  // Return
203  return;
204 }
205 
206 
207 /***********************************************************************//**
208  * @brief Set horizontal direction from 3D vector
209  *
210  * @param[in] vector 3D vector.
211  *
212  * Convert a 3-dimensional vector in into a horizontal direction. The
213  * transformation is given by
214  *
215  * \f[
216  * \alpha = \arctan \left( \frac{x_1}{x_0} \right)
217  * \f]
218  *
219  * \f[
220  * \delta = \arcsin x_2
221  * \f]
222  ***************************************************************************/
223 void GHorizDir::celvector(const GVector& vector)
224 {
225  // Convert vector into horizontal direction
226  m_alt = std::asin(vector[2]);
227  m_az = std::atan2(vector[1], vector[0]);
228 
229  // Return
230  return;
231 }
232 
233 
234 /***********************************************************************//**
235  * @brief Rotate horizontal direction by zenith and azimuth angle
236  *
237  * @param[in] phi Azimuth angle (deg).
238  * @param[in] theta Zenith angle (deg).
239  *
240  * Rotate horizontal direction by a zenith and azimuth angle given in the
241  * system of the horizontal direction and aligned in TBD.
242  * The azimuth angle is counted counter clockwise from TBD.
243  ***************************************************************************/
244 void GHorizDir::rotate_deg(const double& phi, const double& theta)
245 {
246 
247  // Allocate Euler and rotation matrices
248  GMatrix ry;
249  GMatrix rz;
250 
251  // Set up rotation matrix to rotate from native coordinates to
252  // celestial coordinates
253  ry.eulery(m_alt * gammalib::rad2deg - 90.0);
255  GMatrix rot = (ry * rz).transpose();
256 
257  // Set up native coordinate vector
258  double phi_rad = phi * gammalib::deg2rad;
259  double theta_rad = theta * gammalib::deg2rad;
260  double cos_phi = std::cos(phi_rad);
261  double sin_phi = std::sin(phi_rad);
262  double cos_theta = std::cos(theta_rad);
263  double sin_theta = std::sin(theta_rad);
264  GVector native(-cos_phi*sin_theta, sin_phi*sin_theta, cos_theta);
265 
266  // Rotate vector into horizontal coordinates
267  GVector dir = rot * native;
268 
269  // Convert vector into horizontal direction
270  celvector(dir);
271 
272  // Return
273  return;
274 }
275 
276 
277 /***********************************************************************//**
278  * @brief Return horizontal direction as 3D vector
279  *
280  * @return Horizontal direction as 3D vector.
281  ***************************************************************************/
283 {
284  // Compute 3D vector
285  double cosaz = std::cos(m_az);
286  double sinaz = std::sin(m_az);
287  double cosalt = std::cos(m_alt);
288  double sinalt = std::sin(m_alt);
289  GVector vector(cosalt*cosaz, cosalt*sinaz, sinalt);
290 
291  // Return vector
292  return vector;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Compute angular distance to horizontal direction in radians
298  *
299  * @param[in] dir Horizontal direction.
300  * @return Angular distance in radians.
301  *
302  * Computes the angular distance to a specified horizontal direction in
303  * radians.
304  ***************************************************************************/
305 double GHorizDir::dist(const GHorizDir& dir) const
306 {
307  // Compute distance
308  double cosdis = std::sin(m_alt) * std::sin(dir.alt()) +
309  std::cos(m_alt) * std::cos(dir.alt()) *
310  std::cos(dir.az() - m_az);
311 
312  // Compute distance (use argument save GTools function)
313  double dist = gammalib::acos(cosdis);
314 
315  // Return distance
316  return dist;
317 }
318 
319 
320 /***********************************************************************//**
321  * @brief Compute angular distance to horizontal direction in degrees
322  *
323  * @param[in] dir Horizontal direction.
324  * @return Angular distance in degrees.
325  *
326  * Computes the angular distance to a specified horizontal direction in
327  * degrees.
328  ***************************************************************************/
329 double GHorizDir::dist_deg(const GHorizDir& dir) const
330 {
331  // Return distance in degrees
332  return (dist(dir) * gammalib::rad2deg);
333 }
334 
335 
336 
337 /***********************************************************************//**
338  * @brief Print horizontal direction information
339  *
340  * @param[in] chatter Chattiness (defaults to NORMAL).
341  * @return String containing horizontal direction information.
342  ***************************************************************************/
343 std::string GHorizDir::print(const GChatter& chatter) const
344 {
345  // Initialise result string
346  std::string result;
347 
348  // Continue only if chatter is not silent
349  if (chatter != SILENT) {
350 
351  // Put coordinates in string
352  result = "(Az,Alt)=("+gammalib::str(m_az*gammalib::rad2deg) + "," +
353  gammalib::str(m_alt*gammalib::rad2deg)+")";
354 
355  } // endif: chatter was not silent
356 
357  // Return result
358  return result;
359 }
360 
361 
362 /*==========================================================================
363  = =
364  = Private methods =
365  = =
366  ==========================================================================*/
367 
368 /***********************************************************************//**
369  * @brief Initialise class members
370  ***************************************************************************/
372 {
373  // Initialise members
374  m_az = 0.0;
375  m_alt = 0.0;
376 
377  // Return
378  return;
379 }
380 
381 
382 /***********************************************************************//**
383  * @brief Copy class members
384  *
385  * @param[in] dir Horizontal direction.
386  ***************************************************************************/
388 {
389  // Copy members
390  m_az = dir.m_az;
391  m_alt = dir.m_alt;
392 
393  // Return
394  return;
395 }
396 
397 
398 /***********************************************************************//**
399  * @brief Delete class members
400  ***************************************************************************/
402 {
403  // Return
404  return;
405 }
406 
407 
408 /*==========================================================================
409  = =
410  = Friends =
411  = =
412  ==========================================================================*/
413 
414 /***********************************************************************//**
415  * @brief Equality operator
416  *
417  * @param[in] a First horizontal direction.
418  * @param[in] b Second horizontal direction.
419  *
420  * Compare two horizontal directions. If the coordinate is at the pole,
421  * the azimuth value is irrelevant.
422  *
423  * @todo: really should test for being within some tolerance here
424  ***************************************************************************/
425 bool operator==(const GHorizDir &a, const GHorizDir &b)
426 {
427  // Initialise result
428  bool equal = false;
429 
430  // Do check
431  if (std::abs(std::abs(a.m_alt) - 90.0) < 1e-10) {
432  equal = (a.m_alt == b.m_alt);
433  }
434  else {
435  equal = (a.m_alt == b.m_alt && a.m_az == b.m_az);
436  }
437 
438  // Return equality
439  return equal;
440 }
441 
442 
443 /***********************************************************************//**
444  * @brief Non equality operator
445  *
446  * @param[in] a First horizontal direction.
447  * @param[in] b Second horizontal direction.
448  ***************************************************************************/
449 bool operator!=(const GHorizDir &a, const GHorizDir &b)
450 {
451  // Return non equality
452  return (!(a==b));
453 }
const double & alt(void) const
Return altitude angle in radians.
Definition: GHorizDir.hpp:140
void rotate_deg(const double &phi, const double &theta)
Rotate horizontal direction by zenith and azimuth angle.
Definition: GHorizDir.cpp:244
void clear(void)
Clear horizontal direction.
Definition: GHorizDir.cpp:142
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
Generic matrix class definition.
void free_members(void)
Delete class members.
Definition: GHorizDir.cpp:401
double m_az
azimuth in radians
Definition: GHorizDir.hpp:94
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Horizontal (Alt/Az) direction class.
Definition: GHorizDir.hpp:52
Gammalib tools definition.
GHorizDir & operator=(const GHorizDir &dir)
Assignment operator.
Definition: GHorizDir.cpp:112
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
const double deg2rad
Definition: GMath.hpp:43
virtual ~GHorizDir(void)
Destructor.
Definition: GHorizDir.cpp:90
GHorizDir(void)
Constructor.
Definition: GHorizDir.cpp:59
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
GChatter
Definition: GTypemaps.hpp:33
void copy_members(const GHorizDir &dir)
Copy class members.
Definition: GHorizDir.cpp:387
double dist(const GHorizDir &dir) const
Compute angular distance to horizontal direction in radians.
Definition: GHorizDir.cpp:305
Vector class interface definition.
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
std::string print(const GChatter &chatter=NORMAL) const
Print horizontal direction information.
Definition: GHorizDir.cpp:343
GHorizDir * clone(void) const
Clone horizontal direction.
Definition: GHorizDir.cpp:160
void altaz_deg(const double &alt, const double &az)
Set horizontal direction (degrees)
Definition: GHorizDir.cpp:195
void init_members(void)
Initialise class members.
Definition: GHorizDir.cpp:371
void altaz(const double &alt, const double &az)
Set horizontal direction (radians)
Definition: GHorizDir.cpp:175
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
double dist_deg(const GHorizDir &dir) const
Compute angular distance to horizontal direction in degrees.
Definition: GHorizDir.cpp:329
Exception handler interface definition.
GVector celvector(void) const
Return horizontal direction as 3D vector.
Definition: GHorizDir.cpp:282
double m_alt
altitude in radians
Definition: GHorizDir.hpp:93
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
const double & az(void) const
Return azimuth angle in radians.
Definition: GHorizDir.hpp:164
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
const double rad2deg
Definition: GMath.hpp:44
bool operator==(const GEnergy &a, const GEnergy &b)
Energy equality operator friend.
Definition: GEnergy.hpp:297
bool operator!=(const GEbounds &a, const GEbounds &b)
Energy boundaries inequality operator friend.
Definition: GEbounds.hpp:207
Mathematical function definitions.
Horizontal direction class interface definition.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413