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