GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSkyDir.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSkyDir.cpp - Sky direction class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-2020 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 GSkyDir.cpp
23  * @brief Sky direction class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdlib>
32 #include <cmath>
33 #include "GException.hpp"
34 #include "GTools.hpp"
35 #include "GSkyDir.hpp"
36 #include "GTime.hpp"
37 #include "GMatrix.hpp"
38 #include "GVector.hpp"
39 
40 /* __ Method name definitions ____________________________________________ */
41 #define G_POSANG "GSkyDir::posang(GSkyDir&, std::string&)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 /* __ Prototypes _________________________________________________________ */
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Constructor
59  ***************************************************************************/
61 {
62  // Initialise members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief Vector constructor
72  *
73  * @param[in] vector Sky vector.
74  *
75  * Constructs a sky direction from a 3-element celestial vector. See the
76  * celvector() method for more information.
77  ***************************************************************************/
78 GSkyDir::GSkyDir(const GVector& vector)
79 {
80  // Initialise members
81  init_members();
82 
83  // Set sky direction from vector
84  celvector(vector);
85 
86  // Return
87  return;
88 }
89 
90 
91 /***********************************************************************//**
92  * @brief Copy constructor
93  *
94  * @param[in] dir Sky direction.
95  ***************************************************************************/
97 {
98  // Initialise members
99  init_members();
100 
101  // Copy members
102  copy_members(dir);
103 
104  // Return
105  return;
106 }
107 
108 
109 /***********************************************************************//**
110  * @brief Destructor
111  ***************************************************************************/
113 {
114  // Free members
115  free_members();
116 
117  // Return
118  return;
119 }
120 
121 
122 /*==========================================================================
123  = =
124  = Operators =
125  = =
126  ==========================================================================*/
127 
128 /***********************************************************************//**
129  * @brief Assignment operator
130  *
131  * @param[in] dir Sky direction.
132  * @return Sky direction.
133  ***************************************************************************/
135 {
136  // Execute only if object is not identical
137  if (this != &dir) {
138 
139  // Free members
140  free_members();
141 
142  // Initialise members
143  init_members();
144 
145  // Copy members
146  copy_members(dir);
147 
148  } // endif: object was not identical
149 
150  // Return this object
151  return *this;
152 }
153 
154 
155 /*==========================================================================
156  = =
157  = Public methods =
158  = =
159  ==========================================================================*/
160 
161 /***********************************************************************//**
162  * @brief Clear sky direction
163  ***************************************************************************/
164 void GSkyDir::clear(void)
165 {
166  // Free members
167  free_members();
168 
169  // Initialise members
170  init_members();
171 
172  // Return
173  return;
174 }
175 
176 
177 /***********************************************************************//**
178  * @brief Clone sky direction
179  *
180  * @return Pointer to deep copy of sky direction.
181  ***************************************************************************/
183 {
184  // Clone sky direction
185  return new GSkyDir(*this);
186 }
187 
188 
189 /***********************************************************************//**
190  * @brief Set equatorial sky direction (radians)
191  *
192  * @param[in] ra Right Ascension in radians.
193  * @param[in] dec Declination in radians.
194  *
195  * Sets Right Ascension and Declination in radians.
196  ***************************************************************************/
197 void GSkyDir::radec(const double& ra, const double& dec)
198 {
199  // Set attributes
200  m_has_lb = false;
201  m_has_radec = true;
202  #if defined(G_SINCOS_CACHE)
203  m_has_lb_cache = false;
204  m_has_radec_cache = false;
205  #endif
206 
207  // Set direction
208  m_ra = ra;
209  m_dec = dec;
210 
211  // Return
212  return;
213 }
214 
215 
216 /***********************************************************************//**
217  * @brief Set equatorial sky direction (degrees)
218  *
219  * @param[in] ra Right Ascension in degrees.
220  * @param[in] dec Declination in degrees.
221  *
222  * Sets Right Ascension and Declination in degrees.
223  ***************************************************************************/
224 void GSkyDir::radec_deg(const double& ra, const double& dec)
225 {
226  // Set attributes
227  m_has_lb = false;
228  m_has_radec = true;
229  #if defined(G_SINCOS_CACHE)
230  m_has_lb_cache = false;
231  m_has_radec_cache = false;
232  #endif
233 
234  // Set direction
235  m_ra = ra * gammalib::deg2rad;
236  m_dec = dec * gammalib::deg2rad;
237 
238  // Return
239  return;
240 }
241 
242 
243 /***********************************************************************//**
244  * @brief Set galactic sky direction (radians)
245  *
246  * @param[in] l Galactic longitude in radians.
247  * @param[in] b Galactic latitude in radians.
248  *
249  * Sets Galactic longitude and latitude in radians.
250  ***************************************************************************/
251 void GSkyDir::lb(const double& l, const double& b)
252 {
253  // Set attributes
254  m_has_lb = true;
255  m_has_radec = false;
256  #if defined(G_SINCOS_CACHE)
257  m_has_lb_cache = false;
258  m_has_radec_cache = false;
259  #endif
260 
261  // Set direction
262  m_l = l;
263  m_b = b;
264 
265  // Return
266  return;
267 }
268 
269 
270 /***********************************************************************//**
271  * @brief Set galactic sky direction (degrees)
272  *
273  * @param[in] l Galactic longitude in degrees.
274  * @param[in] b Galactic latitude in degrees.
275  *
276  * Sets Galactic longitude and latitude in degrees.
277  ***************************************************************************/
278 void GSkyDir::lb_deg(const double& l, const double& b)
279 {
280  // Set attributes
281  m_has_lb = true;
282  m_has_radec = false;
283  #if defined(G_SINCOS_CACHE)
284  m_has_lb_cache = false;
285  m_has_radec_cache = false;
286  #endif
287 
288  // Set direction
289  m_l = l * gammalib::deg2rad;
290  m_b = b * gammalib::deg2rad;
291 
292  // Return
293  return;
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Set sky direction from 3D vector in celestial coordinates
299  *
300  * @param[in] vector 3D vector.
301  *
302  * Convert a 3-dimensional vector in celestial coordinates into a sky
303  * direction. The transformation is given by
304  *
305  * \f[
306  * \alpha = \arctan \left( \frac{x_1}{x_0} \right)
307  * \f]
308  *
309  * \f[
310  * \delta = \arcsin x_2
311  * \f]
312  ***************************************************************************/
313 void GSkyDir::celvector(const GVector& vector)
314 {
315  // Set attributes
316  m_has_lb = false;
317  m_has_radec = true;
318  #if defined(G_SINCOS_CACHE)
319  m_has_lb_cache = false;
320  m_has_radec_cache = false;
321  #endif
322 
323  // Convert vector into sky position
324  m_dec = std::asin(vector[2]);
325  m_ra = std::atan2(vector[1], vector[0]);
326 
327  // Return
328  return;
329 }
330 
331 
332 /***********************************************************************//**
333  * @brief Rotate sky direction by zenith and azimuth angle
334  *
335  * @param[in] phi Azimuth angle (radians).
336  * @param[in] theta Zenith angle (radians).
337  *
338  * Rotate sky direction by a zenith and azimuth angle given in the system
339  * of the sky direction and aligned in celestial coordinates.
340  * The azimuth angle is counted counter clockwise from celestial north
341  * (this is identical to the astronomical definition of a position angle).
342  ***************************************************************************/
343 void GSkyDir::rotate(const double& phi, const double& theta)
344 {
345  // If we have no equatorial coordinates then get them now
346  if (!m_has_radec && m_has_lb) {
347  gal2equ();
348  }
349 
350  // Allocate Euler and rotation matrices
351  GMatrix ry;
352  GMatrix rz;
353 
354  // Set up rotation matrix to rotate from native coordinates to
355  // celestial coordinates
356  ry.eulery(m_dec * gammalib::rad2deg - 90.0);
358  GMatrix rot = (ry * rz).transpose();
359 
360  // Set up native coordinate vector
361  double cos_phi = std::cos(phi);
362  double sin_phi = std::sin(phi);
363  double cos_theta = std::cos(theta);
364  double sin_theta = std::sin(theta);
365  GVector native(-cos_phi*sin_theta, sin_phi*sin_theta, cos_theta);
366 
367  // Rotate vector into celestial coordinates
368  GVector dir = rot * native;
369 
370  // Convert vector into sky position
371  celvector(dir);
372 
373  // Return
374  return;
375 }
376 
377 
378 /***********************************************************************//**
379  * @brief Rotate sky direction by zenith and azimuth angle
380  *
381  * @param[in] phi Azimuth angle (deg).
382  * @param[in] theta Zenith angle (deg).
383  *
384  * Rotate sky direction by a zenith and azimuth angle given in the system
385  * of the sky direction and aligned in celestial coordinates.
386  * The azimuth angle is counted counter clockwise from celestial north
387  * (this is identical to the astronomical definition of a position angle).
388  ***************************************************************************/
389 void GSkyDir::rotate_deg(const double& phi, const double& theta)
390 {
391  // Convert phi and theta in radians
392  double phi_rad = phi * gammalib::deg2rad;
393  double theta_rad = theta * gammalib::deg2rad;
394 
395  // Rotate
396  rotate(phi_rad, theta_rad);
397 
398  // Return
399  return;
400 }
401 
402 
403 /***********************************************************************//**
404  * @brief Precess sky direction
405  *
406  * @param[in] from_epoch Epoch of the current coordinate.
407  * @param[in] to_epoch Epoch of the returned precessed coordinate.
408  *
409  * Precesses the sky direction from one epoch to another.
410  *
411  * The method is adapted from a set of fortran subroutines based on
412  * (a) pages 30-34 of the Explanatory Supplement to the AE,
413  * (b) Lieske, et al. (1977) A&A 58, 1-16, and
414  * (c) Lieske (1979) A&A 73, 282-284.
415  ***************************************************************************/
416 void GSkyDir::precess(const double& from_epoch, const double& to_epoch)
417 {
418  // Set constants
419  const double arcsec = gammalib::pi / 180.0 / 3600.0;
420 
421  // Continue only if epochs differ
422  if (from_epoch != to_epoch) {
423 
424  // t0, t below correspond to Lieske's big T and little T
425  double t0 = (from_epoch - 2000.0) / 100.0;
426  double t = (to_epoch - from_epoch) / 100.0;
427  double t02 = t0*t0;
428  double t2 = t*t;
429  double t3 = t2*t;
430 
431  // a,b,c below correspond to Lieske's zeta_A, z_A and theta_A
432  double a = ((2306.2181 + 1.39656 * t0 - 0.000139 * t02) * t +
433  (0.30188 - 0.000344 * t0) * t2 + 0.017998 * t3) * arcsec;
434  double b = ((2306.2181 + 1.39656 * t0 - 0.000139 * t02) * t +
435  (1.09468 + 0.000066 * t0) * t2 + 0.018203 * t3) * arcsec;
436  double c = ((2004.3109 - 0.85330 * t0 - 0.000217 * t02) * t +
437  (-0.42665 - 0.000217 * t0) * t2 - 0.041833 * t3) * arcsec;
438 
439  // Compute sines and cosines
440  double sina = std::sin(a);
441  double cosa = std::cos(a);
442  double sinb = std::sin(b);
443  double cosb = std::cos(b);
444  double sinc = std::sin(c);
445  double cosc = std::cos(c);
446 
447  // Setup precession rotation matrix
448  double xx = cosa*cosc*cosb - sina*sinb;
449  double yx = -sina*cosc*cosb - cosa*sinb;
450  double zx = -sinc*cosb;
451  double xy = cosa*cosc*sinb + sina*cosb;
452  double yy = -sina*cosc*sinb + cosa*cosb;
453  double zy = -sinc*sinb;
454  double xz = cosa*sinc;
455  double yz = -sina*sinc;
456  double zz = cosc;
457 
458  // Get vector of sky position
459  GVector vector = celvector();
460 
461  // Perform the rotation:
462  double x2 = xx*vector[0] + yx*vector[1] + zx*vector[2];
463  double y2 = xy*vector[0] + yy*vector[1] + zy*vector[2];
464  double z2 = xz*vector[0] + yz*vector[1] + zz*vector[2];
465 
466  // Setup rotated vector
467  GVector rotated_vector(x2, y2, z2);
468 
469  // Transform vector to sky direction
470  celvector(rotated_vector);
471 
472  } // endif: epochs differed
473 
474  // Return
475  return;
476 }
477 
478 
479 /***********************************************************************//**
480  * @brief Set sky direction to direction of Sun
481  *
482  * @param[in] time Time.
483  * @param[in] epoch Julian epoch.
484  *
485  * Sets the sky direction to the direction of the Sun at a given @p time.
486  *
487  * The equations were taken from the Astronomical Almanac. They calculate the
488  * apparent coordinates of the Sun to a precision of about 0.01 degrees
489  * (36 arcsec), for dates between 1950 and 2050.
490  *
491  * See https://en.wikipedia.org/wiki/Position_of_the_Sun
492  ***************************************************************************/
493 void GSkyDir::sun(const GTime& time, const double& epoch)
494 {
495  // Compute number of days since Greenwich noon, Terrestrial Time, on
496  // 1 January 2000
497  double n = time.jd() - 2451545.0;
498 
499  // Compute obliquity of the ecliptic in radians
500  double eps_rad = (23.439 - 0.0000004 * n) * gammalib::deg2rad;
501 
502  // Compute mean longitude of the Sun in degrees, corrected for the
503  // aberration of light. Put the mean longitude in the interval [0,360[
504  double mean_longitude = 280.460 + 0.9856474 * n;
505  while (mean_longitude < 0.0) {
506  mean_longitude += 360.0;
507  }
508  while (mean_longitude >= 360.0) {
509  mean_longitude -= 360.0;
510  }
511 
512  // Compute the mean anomaly of the Sun in radians
513  double mean_anomaly = (357.528 + 0.9856003 * n) * gammalib::deg2rad;
514 
515  // Compute the ecliptic longitude of the Sun in degrees
516  double ecliptic_longitude = mean_longitude +
517  1.915 * std::sin(mean_anomaly) +
518  0.020 * std::sin(2.0*mean_anomaly);
519 
520  // Compute sine and cosine
521  double ecliptic_longitude_rad = ecliptic_longitude * gammalib::deg2rad;
522  double sin_ecliptic_longitude = std::sin(ecliptic_longitude_rad);
523  double cos_ecliptic_longitude = std::cos(ecliptic_longitude_rad);
524 
525  // Compute Right Ascension and Declination of the Sun in degrees
526  double ra = std::atan2(std::cos(eps_rad) * sin_ecliptic_longitude,
527  cos_ecliptic_longitude) * gammalib::rad2deg;
528  double dec = std::asin(std::sin(eps_rad) * sin_ecliptic_longitude) *
530 
531  // Set position
532  radec_deg(ra, dec);
533 
534  // Precess sky position to requested epoch
535  precess(time.julian_epoch(), epoch);
536 
537  // Return
538  return;
539 }
540 
541 
542 /***********************************************************************//**
543  * @brief Set sky direction to direction of Moon
544  *
545  * @param[in] time Time.
546  * @param[in] epoch Julian epoch.
547  *
548  * Sets the sky direction to the direction of the Moon at a given @p time.
549  *
550  * The equations are derived from the Chapront ELP2000/82 Lunar Theory
551  * (Chapront-Touze' and Chapront, 1983, 124, 50), as described by Jean Meeus
552  * in Chapter 47 of ``Astronomical Algorithms'' (Willmann-Bell, Richmond),
553  * 2nd edition, 1998. Meeus quotes an approximate accuracy of 10 arcsec in
554  * longitude and 4 arcsec in latitude, but he does not give the time range
555  * for this accuracy.
556  *
557  * The method was compared to the results obtained using PyEphem for the
558  * period 2000 - 2050. Differences in Right Ascension and Declination are
559  * always smaller than 0.015 degrees, corresponding to about 1 arcmin.
560  * The differences increase over time, meaning that the largest distances
561  * are reached by the end of the period.
562  ***************************************************************************/
563 void GSkyDir::moon(const GTime& time, const double& epoch)
564 {
565  // Set longitude constants. Note that the cos_lng contants are commented
566  // out since they are actually not used. They could be useful in case
567  // that the distance to the Moon should be computed, therefore we leave
568  // them in comments in the code.
569  const int d_lng[] = { 0, 2, 2, 0, 0, 0, 2, 2, 2, 2,
570  0, 1, 0, 2, 0, 0, 4, 0, 4, 2,
571  2, 1, 1, 2, 2, 4, 2, 0, 2, 2,
572  1, 2, 0, 0, 2, 2, 2, 4, 0, 3,
573  2, 4, 0, 2, 2, 2, 4, 0, 4, 1,
574  2, 0, 1, 3, 4, 2, 0, 1, 2, 2};
575  const int m_lng[] = { 0, 0, 0, 0, 1, 0, 0,-1, 0,-1,
576  1, 0, 1, 0, 0, 0, 0, 0, 0, 1,
577  1, 0, 1,-1, 0, 0, 0, 1, 0,-1,
578  0,-2, 1, 2,-2, 0, 0,-1, 0, 0,
579  1,-1, 2, 2, 1,-1, 0, 0,-1, 0,
580  1, 0, 1, 0, 0,-1, 2, 1, 0,0};
581  const int mp_lng[] = { 1,-1, 0, 2, 0, 0,-2,-1, 1, 0,
582  -1, 0, 1, 0, 1, 1,-1, 3,-2,-1,
583  0,-1, 0, 1, 2, 0,-3,-2,-1,-2,
584  1, 0, 2, 0,-1, 1, 0,-1, 2,-1,
585  1,-2,-1,-1,-2, 0, 1, 4, 0,-2,
586  0, 2, 1,-2,-3, 2, 1,-1, 3,-1};
587  const int f_lng[] = { 0, 0, 0, 0, 0, 2, 0, 0, 0, 0,
588  0, 0, 0,-2, 2,-2, 0, 0, 0, 0,
589  0, 0, 0, 0, 0, 0, 0, 0, 2, 0,
590  0, 0, 0, 0, 0,-2, 2, 0, 2, 0,
591  0, 0, 0, 0, 0,-2, 0, 0, 0, 0,
592  -2,-2, 0, 0, 0, 0, 0, 0, 0,-2};
593  const int sin_lng[] = { 6288774, 1274027, 658314, 213618,-185116,
594  -114332, 58793, 57066, 53322, 45758,
595  -40923, -34720, -30383, 15327, -12528,
596  10980, 10675, 10034, 8548, -7888,
597  -6766, -5163, 4987, 4036, 3994,
598  3861, 3665, -2689, -2602, 2390,
599  -2348, 2236, -2120, -2069, 2048,
600  -1773, -1595, 1215, -1110, -892,
601  -810, 759, -713, -700, 691,
602  596, 549, 537, 520, -487,
603  -399, -381, 351, -340, 330,
604  327, -323, 299, 294, 0};
605  /*
606  const int cos_lng[] = {-20905355,-3699111,-2955968,-569925, 48888,
607  -3149, 246158, -152138,-170733,-204586,
608  -129620, 108743, 104755, 10321, 0,
609  79661, -34782, -23210, -21636, 24208,
610  30824, -8379, -16675, -12831, -10445,
611  -11650, 14403, -7003, 0, 10056,
612  6322, -9884, 5751, 0, -4950,
613  4130, 0, -3958, 0, 3258,
614  2616, -1897, -2117, 2354, 0,
615  0, -1423, -1117, -1571, -1739,
616  0, -4421, 0, 0, 0,
617  0, 1165, 0, 0, 8752};
618  */
619 
620  // Set latitude constants
621  const int d_lat[] = { 0, 0, 0, 2, 2, 2, 2, 0, 2, 0,
622  2, 2, 2, 2, 2, 2, 2, 0, 4, 0,
623  0, 0, 1, 0, 0, 0, 1, 0, 4, 4,
624  0, 4, 2, 2, 2, 2, 0, 2, 2, 2,
625  2, 4, 2, 2, 0, 2, 1, 1, 0, 2,
626  1, 2, 0, 4, 4, 1, 4, 1, 4, 2};
627  const int m_lat[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
628  -1, 0, 0, 1,-1,-1,-1, 1, 0, 1,
629  0, 1, 0, 1, 1, 1, 0, 0, 0, 0,
630  0, 0, 0, 0,-1, 0, 0, 0, 0, 1,
631  1, 0,-1,-2, 0, 1, 1, 1, 1, 1,
632  0,-1, 1, 0,-1, 0, 0, 0,-1,-2};
633  const int mp_lat[] = { 0, 1, 1, 0,-1,-1, 0, 2, 1, 2,
634  0,-2, 1, 0,-1, 0,-1,-1,-1, 0,
635  0,-1, 0, 1, 1, 0, 0, 3, 0,-1,
636  1,-2, 0, 2, 1,-2, 3, 2,-3,-1,
637  0, 0, 1, 0, 1, 1, 0, 0,-2,-1,
638  1,-2, 2,-2,-1, 1, 1,-1, 0, 0};
639  const int f_lat[] = { 1, 1,-1,-1, 1,-1, 1, 1,-1,-1,
640  -1,-1, 1,-1, 1, 1,-1,-1,-1, 1,
641  3, 1, 1, 1,-1,-1,-1, 1,-1, 1,
642  -3, 1,-3,-1,-1, 1,-1, 1,-1, 1,
643  1, 1, 1,-1, 3,-1,-1, 1,-1,-1,
644  1,-1, 1,-1,-1,-1,-1,-1,-1, 1};
645  const int sin_lat[] = {5128122, 280602, 277693, 173237, 55413,
646  46271, 32573, 17198, 9266, 8822,
647  8216, 4324, 4200, -3359, 2463,
648  2211, 2065, -1870, 1828, -1794,
649  -1749, -1565, -1491, -1475, -1410,
650  -1344, -1335, 1107, 1021, 833,
651  777, 671, 607, 596, 491,
652  -451, 439, 422, 421, -366,
653  -351, 331, 315, 302, -283,
654  -229, 223, 223, -220, -220,
655  -185, 181, -177, 176, 166,
656  -164, 132, -119, 115, 107};
657 
658  // Set constants for computation of nutation. The constants have been
659  // taken from https://idlastro.gsfc.nasa.gov/ftp/pro/astro/nutate.pro
660  // which is based on Chapter 22 of "Astronomical Algorithms" by Jean
661  // Meeus (1998, 2nd ed.) which is based on the 1980 IAU Theory of
662  // Nutation and includes all terms larger than 0.0003 arcsec.
663  const int nutate_d_lng[] = { 0,-2, 0, 0, 0, 0,-2, 0, 0,-2,
664  -2,-2, 0, 2, 0, 2, 0, 0,-2, 0,
665  2, 0, 0,-2, 0,-2, 0, 0, 2,-2,
666  0,-2, 0, 0, 2, 2, 0,-2, 0, 2,
667  2,-2,-2, 2, 2, 0,-2,-2, 0,-2,
668  -2, 0,-1,-2, 1, 0, 0,-1, 0, 0,
669  2, 0, 2};
670  const int nutate_m_lng[] = { 0, 0, 0, 0, 1, 0, 1, 0, 0,-1,
671  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
672  0, 0, 0, 0, 0, 0, 0, 2, 0, 2,
673  1, 0,-1, 0, 0, 0, 1, 1,-1, 0,
674  0, 0, 0, 0, 0,-1,-1, 0, 0, 0,
675  1, 0, 0, 1, 0, 0, 0,-1, 1,-1,
676  -1, 0,-1};
677  const int nutate_mp_lng[] = { 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
678  1, 0,-1, 0, 1,-1,-1, 1, 2,-2,
679  0, 2, 2, 1, 0, 0,-1, 0,-1, 0,
680  0, 1, 0, 2,-1, 1, 0, 1, 0, 0,
681  1, 2, 1,-2, 0, 1, 0, 0, 2, 2,
682  0, 1, 1, 0, 0, 1,-2, 1, 1, 1,
683  -1, 3, 0};
684  const int nutate_f_lng[] = { 0, 2, 2, 0, 0, 0, 2, 2, 2, 2,
685  0, 2, 2, 0, 0, 2, 0, 2, 0, 2,
686  2, 2, 0, 2, 2, 2, 2, 0, 0, 2,
687  0, 0, 0,-2, 2, 2, 2, 0, 2, 2,
688  0, 2, 2, 0, 0, 0, 2, 0, 2, 0,
689  2,-2, 0, 0, 0, 2, 2, 0, 0, 2,
690  2, 2, 2};
691  const int nutate_om_lng[] = { 1, 2, 2, 2, 0, 0, 2, 1, 2, 2,
692  0, 1, 2, 0, 1, 2, 1, 1, 0, 1,
693  2, 2, 0, 2, 0, 0, 1, 0, 1, 2,
694  1, 1, 1, 0, 1, 2, 2, 0, 2, 1,
695  0, 2, 1, 1, 1, 0, 1, 1, 1, 1,
696  1, 0, 0, 0, 0, 0, 2, 0, 0, 2,
697  2, 2, 2};
698  const int nutate_sin_lng[] = {-171996, -13187, -2274, 2062, 1426,
699  712, -517, -386, -301, 217,
700  -158, 129, 123, 63, 63,
701  -59, -58, -51, 48, 46,
702  -38, -31, 29, 29, 26,
703  -22, 21, 17, 16, -16,
704  -15, -13, -12, 11, -10,
705  -8, 7, -7, -7, -7,
706  6, 6, 6, -6, -6,
707  5, -5, -5, -5, 4,
708  4, 4, -4, -4, -4,
709  3, -3, -3, -3, -3,
710  -3, -3, -3};
711  const int nutate_cos_lng[] = { 92025, 5736, 977, -895, 54,
712  -7, 224, 200, 129, -95,
713  0, -70, -53, 0, -33,
714  26, 32, 27, 0, -24,
715  16, 13, 0, -12, 0,
716  0, -10, 0, -8, 7,
717  9, 7, 6, 0, 5,
718  3, -3, 0, 3, 3,
719  0, -3, -3, 3, 3,
720  0, 3, 3, 3, 0,
721  0, 0, 0, 0, 0,
722  0, 0, 0, 0, 0,
723  0, 0, 0};
724  const double nutate_sdelt[] = {-174.2, -1.6, -0.2, 0.2, -3.4,
725  0.1, 1.2, -0.4, 0.0, -0.5,
726  0.0, 0.1, 0.0, 0.0, 0.1,
727  0.0, -0.1, 0.0, 0.0, 0.0,
728  0.0, 0.0, 0.0, 0.0, 0.0,
729  0.0, 0.0, -0.1, 0.0, 0.1,
730  0.0, 0.0, 0.0, 0.0, 0.0,
731  0.0, 0.0, 0.0, 0.0, 0.0,
732  0.0, 0.0, 0.0, 0.0, 0.0,
733  0.0, 0.0, 0.0, 0.0, 0.0,
734  0.0, 0.0, 0.0, 0.0, 0.0,
735  0.0, 0.0, 0.0, 0.0, 0.0,
736  0.0, 0.0, 0.0};
737  const double nutate_cdelt[] = {8.9, -3.1, -0.5, 0.5, -0.1,
738  0.0, -0.6, 0.0, -0.1, 0.3,
739  0.0, 0.0, 0.0, 0.0, 0.0,
740  0.0, 0.0, 0.0, 0.0, 0.0,
741  0.0, 0.0, 0.0, 0.0, 0.0,
742  0.0, 0.0, 0.0, 0.0, 0.0,
743  0.0, 0.0, 0.0, 0.0, 0.0,
744  0.0, 0.0, 0.0, 0.0, 0.0,
745  0.0, 0.0, 0.0, 0.0, 0.0,
746  0.0, 0.0, 0.0, 0.0, 0.0,
747  0.0, 0.0, 0.0, 0.0, 0.0,
748  0.0, 0.0, 0.0, 0.0, 0.0,
749  0.0, 0.0, 0.0};
750 
751  // Compute Julian centuries from 1900.0
752  double t = (time.jd() - 2451545.0) / 36525.0;
753 
754  // Compute Moon's mean longitude referred to mean equinox of the date
755  // in degrees and radians
756  double Lprimed = 218.3164477 +
757  (481267.88123421 +
758  (-0.0015786 +
759  (1.0/538841.0 - 1.0/65194000.0 * t) * t) * t) * t;
760  double Lprime = Lprimed * gammalib::deg2rad;
761 
762  // Compute Moon's mean elongation in degrees
763  double D = (297.8501921 +
764  (445267.1114034 +
765  (-0.0018819 +
766  (1.0/545868.0 - 1.0/113065000.0 * t) * t) * t) * t) *
768 
769  // Compute Sun's mean anomaly in radians
770  double M = (357.5291092 +
771  (35999.0502909 +
772  (-0.0001536 + 1.0/24490000.0 * t) * t) * t) *
773  gammalib::deg2rad;
774 
775  // Compute Moon's mean anomaly in radians
776  double Mprime = (134.9633964 +
777  (477198.8675055 +
778  (0.0087414 +
779  (1.0/69699.0 - 1.0/14712000.0 * t) * t) * t) * t) *
781 
782  // Compute Moon's argument of latitude in radians
783  double F = (93.2720950 +
784  (483202.0175233 +
785  (-0.0036539 +
786  (-1.0/35260000 + 1.0/863310000 * t) * t) * t) * t) *
788 
789  // Compute longitude of the ascending node of the Moon's mean orbit on
790  // the ecliptic, measured from the mean equinox of the date in radians
791  double Omega = (125.04452 +
792  (-1934.136261 +
793  (0.0020708 + 1.0/4.5e5 * t) * t) * t) *
794  gammalib::deg2rad;
795 
796  // Compute actions of Venus and Jupiter in radians
797  double A1 = (119.75 + 131.849 * t) * gammalib::deg2rad;
798  double A2 = ( 53.09 + 479264.290 * t) * gammalib::deg2rad;
799  double A3 = (313.45 + 481266.484 * t) * gammalib::deg2rad;
800 
801  // Compute sum quantities
802  double Al = 3958.0 * std::sin(A1) +
803  1962.0 * std::sin(Lprime - F) +
804  318.0 * std::sin(A2);
805  double Ab = -2235.0 * std::sin(Lprime) +
806  382.0 * std::sin(A3) +
807  175.0 * std::sin(A1 - F) +
808  175.0 * std::sin(A1 + F) +
809  127.0 * std::sin(Lprime - Mprime) -
810  115.0 * std::sin(Lprime + Mprime);
811 
812  // Compute eccentricity of Earth's orbit around the Sun
813  double E = 1.0 - (0.002516 + 7.4e-6 * t) * t;
814  double E2 = E * E;
815 
816  // Sum the periodic terms
817  double sum_l = Al;
818  double sum_b = Ab;
819  //double sum_r = 0.0; // A_r = 0
820  for (int i = 0; i < 60; ++i) {
821 
822  // Compute sine and cosine of longitude and sine of latitude
823  double sinlng = sin_lng[i];
824  //double coslng = cos_lng[i];
825  double sinlat = sin_lat[i];
826  if (std::abs(m_lng[i]) == 1) {
827  sinlng *= E;
828  // coslng *= E;
829  }
830  else if (std::abs(m_lng[i]) == 2) {
831  sinlng *= E2;
832  // coslng *= E2;
833  }
834  if (std::abs(m_lat[i]) == 1) {
835  sinlat *= E;
836  }
837  else if (std::abs(m_lat[i]) == 2) {
838  sinlat *= E2;
839  }
840 
841  // Compute coefficient
842  double s_l = d_lng[i] * D + m_lng[i] * M + mp_lng[i] * Mprime +
843  f_lng[i] * F;
844  double s_b = d_lat[i] * D + m_lat[i] * M + mp_lat[i] * Mprime +
845  f_lat[i] * F;
846 
847  // Compute longitude
848  sum_l += sinlng * std::sin(s_l);
849  sum_b += sinlat * std::sin(s_b);
850  //sum_r += coslng * std::sin(s_l);
851 
852  } // endfor: looped over sum terms
853 
854  // Compute Moon's ecliptic coordinates in degrees and distance in km
855  double geolong = Lprimed + 1.0e-6 * sum_l;
856  double geolat = 1.0e-6 * sum_b;
857  //double distance = 385000.56 + 1.0e-3 * sum_r;
858 
859  // Compute nutation in arcsec, based on Chapter 22 of "Astronomical
860  // Algorithms" by Jean Meeus (1998, 2nd ed.) which is based on the 1980
861  // IAU Theory of Nutation and includes all terms larger than 0.0003 arcsec.
862  double nut_long = 0.0;
863  double nut_obliq = 0.0;
864  for (int i = 0; i < 63; ++i) {
865  double arg = nutate_d_lng[i] * D + nutate_m_lng[i] * M +
866  nutate_mp_lng[i] * Mprime + nutate_f_lng[i] * F +
867  nutate_om_lng[i] * Omega;
868  nut_long += (nutate_sdelt[i] * t + nutate_sin_lng[i]) * std::sin(arg);
869  nut_obliq += (nutate_cdelt[i] * t + nutate_cos_lng[i]) * std::cos(arg);
870  }
871  nut_long *= 0.0001;
872  nut_obliq *= 0.0001;
873 
874  // Add nutation in longitude
875  geolong += nut_long / 3600.0;
876 
877  // nutate, jd, nlong, elong
878  // geolong= geolong + nlong/3.6d3
879  // cirrange,geolong
880  double lambda = geolong * gammalib::deg2rad;
881  double beta = geolat * gammalib::deg2rad;
882 
883  // Find mean obliquity using J. Laskar's formula. For the formula, see
884  // http://www.neoprogrammics.com/obliquity_of_the_ecliptic/
885  double y = t / 100.0; // Time in years
886  double epsilon = (84381.448 + // In degrees
887  (4680.93 +
888  (-1.55 +
889  (1999.25 +
890  (-51.38 +
891  (-249.67 +
892  (-39.05 +
893  (7.12 +
894  (27.87 +
895  (5.79 + 2.45 * y)*y)*y)*y)*y)*y)*y)*y)*y)*y)/3600.0;
896  double eps = (epsilon + nut_obliq/3600.0) * gammalib::deg2rad;
897 
898  // Convert (lambda,beta) in radians to (RA,Dec) in degrees
899  double ra = std::atan2(std::sin(lambda) * std::cos(eps) -
900  std::tan(beta) * std::sin(eps), std::cos(lambda)) *
902  double dec = std::asin(std::sin(beta) * std::cos(eps) +
903  std::cos(beta) * std::sin(eps) * std::sin(lambda)) *
905 
906  // Set position
907  radec_deg(ra, dec);
908 
909  // Precess sky position to requested epoch
910  precess(time.julian_epoch(), epoch);
911 
912  // Return
913  return;
914 }
915 
916 
917 /***********************************************************************//**
918  * @brief Return sky direction as 3D vector in celestial coordinates
919  *
920  * @return Sky direction as 3D vector in celestial coordinates.
921  ***************************************************************************/
923 {
924  // If we have no equatorial coordinates then get them now
925  if (!m_has_radec && m_has_lb) {
926  gal2equ();
927  }
928 
929  // Compute 3D vector
930  double cosra = std::cos(m_ra);
931  double sinra = std::sin(m_ra);
932  #if defined(G_SINCOS_CACHE)
933  if (!m_has_radec_cache) {
934  m_has_radec_cache = true;
937  }
938  GVector vector(m_cos_dec*cosra, m_cos_dec*sinra, m_sin_dec);
939  #else
940  double cosdec = std::cos(m_dec);
941  double sindec = std::sin(m_dec);
942  GVector vector(cosdec*cosra, cosdec*sinra, sindec);
943  #endif
944 
945  // Return vector
946  return vector;
947 }
948 
949 
950 /***********************************************************************//**
951  * @brief Compute cosine of angular distance between sky directions
952  *
953  * @param[in] dir Sky direction to which cosine of distance is to be computed.
954  * @return Cosine of angular distance.
955  *
956  * Computes the cosine of the angular distance between two sky directions in
957  * radians.
958  ***************************************************************************/
959 double GSkyDir::cos_dist(const GSkyDir& dir) const
960 {
961  // Initialise cosine of distance
962  double cosdis;
963 
964  // Compute dependent on coordinate system availability. This speeds
965  // up things by avoiding unnecessary coordinate transformations.
966  if (m_has_lb) {
967  #if defined(G_SINCOS_CACHE)
968  if (!m_has_lb_cache) {
969  m_has_lb_cache = true;
970  m_sin_b = std::sin(m_b);
971  m_cos_b = std::cos(m_b);
972  }
973  #endif
974  if (dir.m_has_lb) {
975  #if defined(G_SINCOS_CACHE)
976  if (!dir.m_has_lb_cache) {
977  dir.m_has_lb_cache = true;
978  dir.m_sin_b = std::sin(dir.m_b);
979  dir.m_cos_b = std::cos(dir.m_b);
980  }
981  cosdis = m_sin_b * dir.m_sin_b +
982  m_cos_b * dir.m_cos_b *
983  std::cos(dir.m_l - m_l);
984  #else
985  cosdis = std::sin(m_b) * std::sin(dir.m_b) +
986  std::cos(m_b) * std::cos(dir.m_b) *
987  std::cos(dir.m_l - m_l);
988  #endif
989  }
990  else {
991  #if defined(G_SINCOS_CACHE)
992  cosdis = m_sin_b * std::sin(dir.b()) +
993  m_cos_b * std::cos(dir.b()) *
994  std::cos(dir.l() - m_l);
995  #else
996  cosdis = std::sin(m_b) * std::sin(dir.b()) +
997  std::cos(m_b) * std::cos(dir.b()) *
998  std::cos(dir.l() - m_l);
999  #endif
1000  }
1001  }
1002  else if (m_has_radec) {
1003  #if defined(G_SINCOS_CACHE)
1004  if (!m_has_radec_cache) {
1005  m_has_radec_cache = true;
1008  }
1009  #endif
1010  if (dir.m_has_radec) {
1011  #if defined(G_SINCOS_CACHE)
1012  if (!dir.m_has_radec_cache) {
1013  dir.m_has_radec_cache = true;
1014  dir.m_sin_dec = std::sin(dir.m_dec);
1015  dir.m_cos_dec = std::cos(dir.m_dec);
1016  }
1017  cosdis = m_sin_dec * dir.m_sin_dec +
1018  m_cos_dec * dir.m_cos_dec *
1019  std::cos(dir.m_ra - m_ra);
1020  #else
1021  cosdis = std::sin(m_dec) * std::sin(dir.m_dec) +
1022  std::cos(m_dec) * std::cos(dir.m_dec) *
1023  std::cos(dir.m_ra - m_ra);
1024  #endif
1025  }
1026  else {
1027  #if defined(G_SINCOS_CACHE)
1028  cosdis = m_sin_dec * sin(dir.dec()) +
1029  m_cos_dec * cos(dir.dec()) *
1030  std::cos(dir.ra() - m_ra);
1031  #else
1032  cosdis = sin(m_dec) * sin(dir.dec()) +
1033  cos(m_dec) * cos(dir.dec()) *
1034  std::cos(dir.ra() - m_ra);
1035  #endif
1036  }
1037  }
1038  else {
1039  cosdis = std::sin(dec()) * std::sin(dir.dec()) +
1040  std::cos(dec()) * std::cos(dir.dec()) *
1041  std::cos(dir.ra() - ra());
1042  }
1043 
1044  // Return cosine of distance
1045  return cosdis;
1046 }
1047 
1048 
1049 /***********************************************************************//**
1050  * @brief Compute position angle between sky directions in radians
1051  *
1052  * @param[in] dir Sky direction.
1053  * @param[in] coordsys Coordinate system ("CEL" or "GAL")
1054  * @return Position angle (rad).
1055  *
1056  * Computes the position angle of a specified sky direction with respect to
1057  * the sky direction of the instance in radians. If "CEL" is specified as
1058  * @p coordsys the position angle is counted counterclockwise from celestial
1059  * North, if "GAL" is specified the position angle is counted
1060  * counterclockwise from Galactic North.
1061  *
1062  * If @p coordsys is "CEL" the position angle is computed using
1063  *
1064  * \f[PA = \arctan \left(
1065  * \frac{\sin( \alpha - \alpha_0 )}
1066  * {\cos \delta_0 \tan \delta -
1067  * \sin \delta_0 \cos(\alpha - \alpha_0)} \right)\f]
1068  * where
1069  * * \f$(\alpha_0,\delta_0)\f$ are the celestial coordinates of the instance,
1070  * and
1071  * * \f$(\alpha,\delta)\f$ are the celestial coordinates of sky direction.
1072  *
1073  * If @p coordsys is "GAL" the position angle is computed using
1074  *
1075  * \f[PA = \arctan \left(
1076  * \frac{\sin( \l - \l_0 )}
1077  * {\cos \b_0 \tan \b - \sin \b_0 \cos(\l - \l_0)} \right)\f]
1078  * where
1079  * * \f$(\l_0,\l_0)\f$ are the Galactic coordinates of the instance,
1080  * and
1081  * * \f$(\l,\b)\f$ are the Galactic coordinates of sky direction.
1082  ***************************************************************************/
1083 double GSkyDir::posang(const GSkyDir& dir, const std::string& coordsys) const
1084 {
1085  // Initialise arguments of arctan
1086  double arg_1;
1087  double arg_2;
1088 
1089  // Case A: Computation in celestial coordinates
1090  if (coordsys == "CEL") {
1091 
1092  // If required then compute sin(dec) and cos(dec) cache values
1093  #if defined(G_SINCOS_CACHE)
1094  if (!m_has_radec_cache) {
1095  m_has_radec_cache = true;
1096  m_sin_dec = std::sin(dec());
1097  m_cos_dec = std::cos(dec());
1098  }
1099  #endif
1100 
1101  // Compute arguments of arctan
1102  arg_1 = std::sin(dir.ra() - ra());
1103  #if defined(G_SINCOS_CACHE)
1104  arg_2 = m_cos_dec * std::tan(dir.dec()) -
1105  m_sin_dec * std::cos(dir.ra() - ra());
1106  #else
1107  arg_2 = std::cos(dec()) * std::tan(dir.dec()) -
1108  std::sin(dec()) * std::cos(dir.ra() - ra());
1109  #endif
1110 
1111  } // endif: computation in celestial coordinates
1112 
1113  // Case B: Computation in Galactic coordinates
1114  else if (coordsys == "GAL") {
1115 
1116  // If required then compute sin(b) and cos(b) cache values
1117  #if defined(G_SINCOS_CACHE)
1118  if (!m_has_lb_cache) {
1119  m_has_lb_cache = true;
1120  m_sin_b = std::sin(b());
1121  m_cos_b = std::cos(b());
1122  }
1123  #endif
1124 
1125  // Compute arguments of arctan
1126  arg_1 = std::sin(dir.l() - l());
1127  #if defined(G_SINCOS_CACHE)
1128  arg_2 = m_cos_b * std::tan(dir.b()) -
1129  m_sin_b * std::cos(dir.l() - l());
1130  #else
1131  arg_2 = std::cos(b()) * std::tan(dir.b()) -
1132  std::sin(b()) * std::cos(dir.l() - l());
1133  #endif
1134 
1135  } // endelse: computation in Galactic coordinates
1136 
1137  // ... otherwise throw an exception
1138  else {
1139  std::string msg = "Invalid coordinate system \""+coordsys+"\" "
1140  "specified. Please specify either \"CEL\" or "
1141  "\"GAL\".";
1143  }
1144 
1145  // Compute position angle
1146  double pa = std::atan2(arg_1, arg_2);
1147 
1148  // Return position angle
1149  return pa;
1150 }
1151 
1152 
1153 /***********************************************************************//**
1154  * @brief Print sky direction information
1155  *
1156  * @param[in] chatter Chattiness.
1157  * @return String containing sky direction information.
1158  ***************************************************************************/
1159 std::string GSkyDir::print(const GChatter& chatter) const
1160 {
1161  // Initialise result string
1162  std::string result;
1163 
1164  // Continue only if chatter is not silent
1165  if (chatter != SILENT) {
1166 
1167  // Put coordinates in string
1168  if (m_has_lb) {
1169  result = "(l,b)=("+gammalib::str(m_l*gammalib::rad2deg) + "," +
1170  gammalib::str(m_b*gammalib::rad2deg)+")";
1171  }
1172  else if (m_has_radec) {
1173  result = "(RA,Dec)=("+gammalib::str(m_ra*gammalib::rad2deg) + "," +
1174  gammalib::str(m_dec*gammalib::rad2deg)+")";
1175  }
1176  else {
1177  result = "(RA,Dec)=(not initialised)";
1178  }
1179 
1180  } // endif: chatter was not silent
1181 
1182  // Return result
1183  return result;
1184 }
1185 
1186 
1187 /*==========================================================================
1188  = =
1189  = Private methods =
1190  = =
1191  ==========================================================================*/
1192 
1193 /***********************************************************************//**
1194  * @brief Initialise class members
1195  ***************************************************************************/
1197 {
1198  // Initialise members
1199  m_has_lb = false;
1200  m_has_radec = false;
1201  m_l = 0.0;
1202  m_b = 0.0;
1203  m_ra = 0.0;
1204  m_dec = 0.0;
1205 
1206  // Initialise sincos cache
1207  #if defined(G_SINCOS_CACHE)
1208  m_has_lb_cache = false;
1209  m_has_radec_cache = false;
1210  m_sin_b = 0.0;
1211  m_cos_b = 0.0;
1212  m_sin_dec = 0.0;
1213  m_cos_dec = 0.0;
1214  #endif
1215 
1216  // Return
1217  return;
1218 }
1219 
1220 
1221 /***********************************************************************//**
1222  * @brief Copy class members
1223  *
1224  * @param[in] dir Sky direction.
1225  ***************************************************************************/
1227 {
1228  // Copy members
1229  m_has_lb = dir.m_has_lb;
1230  m_has_radec = dir.m_has_radec;
1231  m_l = dir.m_l;
1232  m_b = dir.m_b;
1233  m_ra = dir.m_ra;
1234  m_dec = dir.m_dec;
1235 
1236  // Copy sincos cache
1237  #if defined(G_SINCOS_CACHE)
1240  m_sin_b = dir.m_sin_b;
1241  m_cos_b = dir.m_cos_b;
1242  m_sin_dec = dir.m_sin_dec;
1243  m_cos_dec = dir.m_cos_dec;
1244  #endif
1245 
1246  // Return
1247  return;
1248 }
1249 
1250 
1251 /***********************************************************************//**
1252  * @brief Delete class members
1253  ***************************************************************************/
1255 {
1256  // Return
1257  return;
1258 }
1259 
1260 
1261 /***********************************************************************//**
1262  * @brief Convert equatorial to galactic coordinates
1263  ***************************************************************************/
1264 void GSkyDir::equ2gal(void) const
1265 {
1266  // Get non-const pointers to data members. This allows to circumvent
1267  // the const correctness and allows treating GSkyDir access methods
1268  // as const
1269  double* l = (double*)&m_l;
1270  double* b = (double*)&m_b;
1271 
1272  // Convert from equatorial to galactic
1273  euler(0, m_ra, m_dec, l, b);
1274 
1275  // Signal that galactic coordinates are available
1276  m_has_lb = true;
1277 
1278  // Return
1279  return;
1280 }
1281 
1282 
1283 /***********************************************************************//**
1284  * @brief Convert galactic to equatorial coordinates
1285  ***************************************************************************/
1286 void GSkyDir::gal2equ(void) const
1287 {
1288  // Get non-const pointers to data members. This allows to circumvent
1289  // the const correctness and allows treating GSkyDir access methods
1290  // as const
1291  double* ra = (double*)&m_ra;
1292  double* dec = (double*)&m_dec;
1293 
1294  // Convert from galactic to equatorial
1295  euler(1, m_l, m_b, ra, dec);
1296 
1297  // Signal that equatorial coordinates are available
1298  m_has_radec = true;
1299 
1300  // Return
1301  return;
1302 }
1303 
1304 
1305 /***********************************************************************//**
1306  * @brief General coordinate transformation routine for J2000
1307  *
1308  * @param[in] type Conversion type (0=equ2gal, 1=gal2equ)
1309  * @param[in] xin Input longitude (RA or GLON) in radians.
1310  * @param[in] yin Input latitude (Dec or GLAT) in radians.
1311  * @param[out] xout Output longitude in radians.
1312  * @param[out] yout Output latitude in radians.
1313  ***************************************************************************/
1314 void GSkyDir::euler(const int& type, const double& xin, const double &yin,
1315  double* xout, double *yout) const
1316 {
1317  // Set transformation constants
1318  const double psi[] = {0.57477043300, 4.9368292465};
1319  const double stheta[] = {0.88998808748, -0.88998808748};
1320  const double ctheta[] = {0.45598377618, 0.45598377618};
1321  const double phi[] = {4.9368292465, 0.57477043300};
1322 
1323  // Perform transformation
1324  double a = xin - phi[type];
1325  double b = yin;
1326  double sb = std::sin(b);
1327  double cb = std::cos(b);
1328  double cbsa = cb * std::sin(a);
1329 
1330  //
1331  a = std::atan2(ctheta[type] * cbsa + stheta[type] * sb, cb * std::cos(a));
1332  b = -stheta[type] * cbsa + ctheta[type] * sb;
1333  if (b > 1.0)
1334  b = 1.0;
1335 
1336  //
1337  *yout = std::asin(b);
1338  *xout = gammalib::modulo((a+psi[type] + gammalib::fourpi), gammalib::twopi);
1339 
1340  // Return
1341  return;
1342 }
1343 
1344 
1345 /*==========================================================================
1346  = =
1347  = Friends =
1348  = =
1349  ==========================================================================*/
1350 
1351 /***********************************************************************//**
1352  * @brief Equality operator
1353  *
1354  * @param[in] a First sky direction.
1355  * @param[in] b Second sky direction.
1356  *
1357  * Compare two sky directions. If the coordinate is at the pole, the Right
1358  * Ascension or Longitude value is irrelevant.
1359  *
1360  * Comparisons are done dependent on the available coordinate system. This
1361  * speeds up things and avoids unnecessary coordinate transformations.
1362  ***************************************************************************/
1363 bool operator==(const GSkyDir &a, const GSkyDir &b)
1364 {
1365  // Initialise result
1366  bool equal = false;
1367 
1368  // Compute dependent on coordinate system availability. This speeds
1369  // up things by avoiding unnecessary coordinate transformations.
1370 
1371  // Check if both have equatorial coordinates
1372  if (a.m_has_radec && b.m_has_radec) {
1373  if (std::abs(a.m_dec) == 90.0) {
1374  equal = (a.m_dec == b.m_dec);
1375  }
1376  else {
1377  equal = (a.m_dec == b.m_dec && a.m_ra == b.m_ra);
1378  }
1379  }
1380  // ... check if both have Galactic coordinates
1381  else if (a.m_has_lb && b.m_has_lb) {
1382  if (std::abs(a.m_b) == 90.0) {
1383  equal = (a.m_b == b.m_b);
1384  }
1385  else {
1386  equal = (a.m_b == b.m_b && a.m_l == b.m_l);
1387  }
1388  }
1389  // ... otherwise the coordinate systems are different
1390  else {
1391  if (a.m_has_lb) {
1392  if (std::abs(a.m_b) == 90.0) {
1393  equal = (a.m_b == b.b());
1394  }
1395  else {
1396  equal = (a.m_b == b.b() && a.m_l == b.l());
1397  }
1398  }
1399  else if (a.m_has_radec) {
1400  if (std::abs(a.m_dec) == 90.0) {
1401  equal = (a.m_dec == b.dec());
1402  }
1403  else {
1404  equal = (a.m_dec == b.dec() && a.m_ra == b.ra());
1405  }
1406  }
1407  else {
1408  if (std::abs(b.dec()) == 90.0) {
1409  equal = (b.dec() == a.dec());
1410  }
1411  else {
1412  equal = (b.dec() == a.dec() && b.ra() == a.ra());
1413  }
1414  }
1415  }
1416  // Return equality
1417  return equal;
1418 }
1419 
1420 
1421 /***********************************************************************//**
1422  * @brief Non equality operator
1423  *
1424  * @param[in] a First sky direction.
1425  * @param[in] b Second sky direction.
1426  ***************************************************************************/
1427 bool operator!=(const GSkyDir &a, const GSkyDir &b)
1428 {
1429  // Return non equality
1430  return (!(a==b));
1431 }
double m_sin_dec
Definition: GSkyDir.hpp:134
double m_cos_b
Definition: GSkyDir.hpp:133
GSkyDir & operator=(const GSkyDir &dir)
Assignment operator.
Definition: GSkyDir.cpp:134
const double pi
Definition: GMath.hpp:35
Sky direction class interface definition.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1170
void lb(const double &l, const double &b)
Set galactic sky direction (radians)
Definition: GSkyDir.cpp:251
double m_sin_b
Definition: GSkyDir.hpp:132
void moon(const GTime &time, const double &epoch=2000.0)
Set sky direction to direction of Moon.
Definition: GSkyDir.cpp:563
const double & b(void) const
Return galactic latitude in radians.
Definition: GSkyDir.hpp:185
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1107
Generic matrix class definition.
void gal2equ(void) const
Convert galactic to equatorial coordinates.
Definition: GSkyDir.cpp:1286
double m_dec
Declination in radians.
Definition: GSkyDir.hpp:126
double modulo(const double &v1, const double &v2)
Returns the remainder of the division.
Definition: GMath.cpp:526
GSkyDir * clone(void) const
Clone sky direction.
Definition: GSkyDir.cpp:182
bool m_has_radec_cache
Definition: GSkyDir.hpp:131
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
bool m_has_lb_cache
Definition: GSkyDir.hpp:130
std::string print(const GChatter &chatter=NORMAL) const
Print sky direction information.
Definition: GSkyDir.cpp:1159
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
const double & ra(void) const
Return Right Ascension in radians.
Definition: GSkyDir.hpp:212
double jd(void) const
Return time in Julian Days (TT)
Definition: GTime.cpp:283
virtual ~GSkyDir(void)
Destructor.
Definition: GSkyDir.cpp:112
double m_b
Galactic latitude in radians.
Definition: GSkyDir.hpp:124
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
Definition: GSkyDir.cpp:197
void copy_members(const GSkyDir &dir)
Copy class members.
Definition: GSkyDir.cpp:1226
const double deg2rad
Definition: GMath.hpp:43
void sun(const GTime &time, const double &epoch=2000.0)
Set sky direction to direction of Sun.
Definition: GSkyDir.cpp:493
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1296
#define G_POSANG
Definition: GSkyDir.cpp:41
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:389
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
void init_members(void)
Initialise class members.
Definition: GSkyDir.cpp:1196
void rotate(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:343
bool m_has_lb
Has galactic coordinates.
Definition: GSkyDir.hpp:121
GChatter
Definition: GTypemaps.hpp:33
const double & l(void) const
Return galactic longitude in radians.
Definition: GSkyDir.hpp:158
const double & dec(void) const
Return Declination in radians.
Definition: GSkyDir.hpp:239
GSkyDir(void)
Constructor.
Definition: GSkyDir.cpp:60
void free_members(void)
Delete class members.
Definition: GSkyDir.cpp:1254
Vector class interface definition.
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1023
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1083
double cos_dist(const GSkyDir &dir) const
Compute cosine of angular distance between sky directions.
Definition: GSkyDir.cpp:959
GVector celvector(void) const
Return sky direction as 3D vector in celestial coordinates.
Definition: GSkyDir.cpp:922
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:164
void euler(const int &type, const double &xin, const double &yin, double *xout, double *yout) const
General coordinate transformation routine for J2000.
Definition: GSkyDir.cpp:1314
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:278
double m_l
Galactic longitude in radians.
Definition: GSkyDir.hpp:123
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1233
Exception handler interface definition.
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
const double fourpi
Definition: GMath.hpp:37
bool m_has_radec
Has equatorial coordinates.
Definition: GSkyDir.hpp:122
const double twopi
Definition: GMath.hpp:36
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
double m_cos_dec
Definition: GSkyDir.hpp:135
const double rad2deg
Definition: GMath.hpp:44
bool operator==(const GEnergy &a, const GEnergy &b)
Energy equality operator friend.
Definition: GEnergy.hpp:297
double m_ra
Right Ascension in radians.
Definition: GSkyDir.hpp:125
Sky direction class.
Definition: GSkyDir.hpp:62
double julian_epoch(void) const
Return Julian epoch in native reference (TT)
Definition: GTime.cpp:424
Time class interface definition.
void precess(const double &from_epoch, const double &to_epoch)
Precess sky direction.
Definition: GSkyDir.cpp:416
bool operator!=(const GEbounds &a, const GEbounds &b)
Energy boundaries inequality operator friend.
Definition: GEbounds.hpp:212
void equ2gal(void) const
Convert equatorial to galactic coordinates.
Definition: GSkyDir.cpp:1264
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415