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