Justine Devin / gammalib

Commit 8a3a3b13a0cb56f9f35e11e6ec3c7ae57713198b

Authored by Jürgen Knödlseder
1 parent 775699ab

Implement sin() and cos() pre computations cache. 

A pre computation cache has been implement for the sin() and the cos() of the declination or galactic latitiude. The values of the cache are used in distance computation, speeding up the repeated usage of the dist() method.
Inline Side-by-side
Showing 5 changed files with 411 additions and 172 deletions
ChangeLog
  View file @8a3a3b1
1   -2013-03-14 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2013-03-29 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.8.0 released
4 4 ========================
... ... @@ -9,6 +9,7 @@
9 9 Add Virtual Observatory support
10 10 Introduce GUrl classes for transparent URL handling
11 11 Make use of GUrl class in XML I/O
  12 + Optimize computations
12 13  
13 14 2012-12-18 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
14 15  
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -March, 18th, 2013
  5 +March, 29th, 2013
6 6  
7 7  
8 8 1. Introduction
... ... @@ -139,7 +139,9 @@ None
139 139  
140 140 12. Sky module
141 141 --------------
142   -None
  142 +Add a sin() and cos() pre computation cache to speed up repeated
  143 +distance computations using dist() or dist_deg(). Coordinate access
  144 +methods are now declared inline to provide also faster data access.
143 145  
144 146  
145 147 13. Support module
... ...
include/GSkyDir.hpp
  View file @8a3a3b1
1 1 /***************************************************************************
2   - * GSkyDir.hpp - Class that implements a sky direction *
  2 + * GSkyDir.hpp - Sky direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2012 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,6 +32,9 @@
32 32 #include "GBase.hpp"
33 33 #include "GVector.hpp"
34 34  
  35 +/* __ Compile options ____________________________________________________ */
  36 +#define G_SINCOS_CACHE
  37 +
35 38  
36 39 /***********************************************************************//**
37 40 * @class GSkyDir
... ... @@ -63,28 +66,28 @@ public:
63 66 GSkyDir& operator= (const GSkyDir& dir);
64 67  
65 68 // Methods
66   - void clear(void);
67   - GSkyDir* clone(void) const;
68   - void radec(const double& ra, const double& dec);
69   - void radec_deg(const double& ra, const double& dec);
70   - void lb(const double& l, const double& b);
71   - void lb_deg(const double& l, const double& b);
72   - void celvector(const GVector& vector);
73   - void rotate_deg(const double& phi, const double& theta);
74   - double l(void) const;
75   - double l_deg(void) const;
76   - double b(void) const;
77   - double b_deg(void) const;
78   - double ra(void) const;
79   - double ra_deg(void) const;
80   - double dec(void) const;
81   - double dec_deg(void) const;
82   - GVector celvector(void) const;
83   - double dist(const GSkyDir& dir) const;
84   - double dist_deg(const GSkyDir& dir) const;
85   - double posang(const GSkyDir& dir) const;
86   - double posang_deg(const GSkyDir& dir) const;
87   - std::string print(void) const;
  69 + void clear(void);
  70 + GSkyDir* clone(void) const;
  71 + void radec(const double& ra, const double& dec);
  72 + void radec_deg(const double& ra, const double& dec);
  73 + void lb(const double& l, const double& b);
  74 + void lb_deg(const double& l, const double& b);
  75 + void celvector(const GVector& vector);
  76 + void rotate_deg(const double& phi, const double& theta);
  77 + const double& l(void) const;
  78 + const double& b(void) const;
  79 + const double& ra(void) const;
  80 + const double& dec(void) const;
  81 + double l_deg(void) const;
  82 + double b_deg(void) const;
  83 + double ra_deg(void) const;
  84 + double dec_deg(void) const;
  85 + GVector celvector(void) const;
  86 + double dist(const GSkyDir& dir) const;
  87 + double dist_deg(const GSkyDir& dir) const;
  88 + double posang(const GSkyDir& dir) const;
  89 + double posang_deg(const GSkyDir& dir) const;
  90 + std::string print(void) const;
88 91  
89 92 private:
90 93 // Private methods
... ... @@ -103,6 +106,77 @@ private:
103 106 double m_b; //!< Galactic latitude in radians
104 107 double m_ra; //!< Right Ascension in radians
105 108 double m_dec; //!< Declination in radians
  109 +
  110 + // Sincos cache
  111 + #if defined(G_SINCOS_CACHE)
  112 + mutable bool m_has_lb_cache;
  113 + mutable bool m_has_radec_cache;
  114 + mutable double m_sin_b;
  115 + mutable double m_cos_b;
  116 + mutable double m_sin_dec;
  117 + mutable double m_cos_dec;
  118 + #endif
106 119 };
107 120  
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return galactic longitude in radians
  124 + *
  125 + * @return Galactic longitude in radians.
  126 + ***************************************************************************/
  127 +inline
  128 +const double& GSkyDir::l(void) const
  129 +{
  130 + if (!m_has_lb && m_has_radec) {
  131 + equ2gal();
  132 + }
  133 + return m_l;
  134 +}
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Return galactic latitude in radians
  139 + *
  140 + * @return Galactic latitude in radians.
  141 + ***************************************************************************/
  142 +inline
  143 +const double& GSkyDir::b(void) const
  144 +{
  145 + if (!m_has_lb && m_has_radec) {
  146 + equ2gal();
  147 + }
  148 + return m_b;
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Return Right Ascension in radians
  154 + *
  155 + * @return Right Ascension in radians.
  156 + ***************************************************************************/
  157 +inline
  158 +const double& GSkyDir::ra(void) const
  159 +{
  160 + if (!m_has_radec && m_has_lb) {
  161 + gal2equ();
  162 + }
  163 + return m_ra;
  164 +}
  165 +
  166 +
  167 +/***********************************************************************//**
  168 + * @brief Return Declination in radians
  169 + *
  170 + * @return Declination in radians.
  171 + ***************************************************************************/
  172 +inline
  173 +const double& GSkyDir::dec(void) const
  174 +{
  175 + if (!m_has_radec && m_has_lb) {
  176 + gal2equ();
  177 + }
  178 + return m_dec;
  179 +}
  180 +
  181 +
108 182 #endif /* GSKYDIR_HPP */
... ...
pyext/GSkyDir.i
  View file @8a3a3b1
1 1 /***************************************************************************
2   - * GSkyDir.i - Sky direction class *
  2 + * GSkyDir.i - Sky direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2012 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -44,27 +44,27 @@ public:
44 44 virtual ~GSkyDir(void);
45 45  
46 46 // Methods
47   - void clear(void);
48   - GSkyDir* clone(void) const;
49   - void radec(const double& ra, const double& dec);
50   - void radec_deg(const double& ra, const double& dec);
51   - void lb(const double& l, const double& b);
52   - void lb_deg(const double& l, const double& b);
53   - void celvector(const GVector& vector);
54   - void rotate_deg(const double& phi, const double& theta);
55   - double l(void) const;
56   - double l_deg(void) const;
57   - double b(void) const;
58   - double b_deg(void) const;
59   - double ra(void) const;
60   - double ra_deg(void) const;
61   - double dec(void) const;
62   - double dec_deg(void) const;
63   - GVector celvector(void) const;
64   - double dist(const GSkyDir& dir) const;
65   - double dist_deg(const GSkyDir& dir) const;
66   - double posang(const GSkyDir& dir) const;
67   - double posang_deg(const GSkyDir& dir) const;
  47 + void clear(void);
  48 + GSkyDir* clone(void) const;
  49 + void radec(const double& ra, const double& dec);
  50 + void radec_deg(const double& ra, const double& dec);
  51 + void lb(const double& l, const double& b);
  52 + void lb_deg(const double& l, const double& b);
  53 + void celvector(const GVector& vector);
  54 + void rotate_deg(const double& phi, const double& theta);
  55 + const double& l(void) const;
  56 + const double& b(void) const;
  57 + const double& ra(void) const;
  58 + const double& dec(void) const;
  59 + double l_deg(void) const;
  60 + double b_deg(void) const;
  61 + double ra_deg(void) const;
  62 + double dec_deg(void) const;
  63 + GVector celvector(void) const;
  64 + double dist(const GSkyDir& dir) const;
  65 + double dist_deg(const GSkyDir& dir) const;
  66 + double posang(const GSkyDir& dir) const;
  67 + double posang_deg(const GSkyDir& dir) const;
68 68 };
69 69  
70 70  
... ...
src/sky/GSkyDir.cpp
  View file @8a3a3b1
1 1 /***************************************************************************
2   - * GSkyDir.cpp - Class that implements a sky direction *
  2 + * GSkyDir.hpp - Sky direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2012 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -105,6 +105,7 @@ GSkyDir::~GSkyDir(void)
105 105 * @brief Assignment operator
106 106 *
107 107 * @param[in] dir Sky direction.
  108 + * @return Sky direction.
108 109 ***************************************************************************/
109 110 GSkyDir& GSkyDir::operator= (const GSkyDir& dir)
110 111 {
... ... @@ -150,11 +151,13 @@ void GSkyDir::clear(void)
150 151  
151 152  
152 153 /***********************************************************************//**
153   - * @brief Clone object
  154 + * @brief Clone sky direction
  155 + *
  156 + * @return Pointer to deep copy of sky direction.
154 157 ***************************************************************************/
155 158 GSkyDir* GSkyDir::clone(void) const
156 159 {
157   - // Clone this image
  160 + // Clone sky direction
158 161 return new GSkyDir(*this);
159 162 }
160 163  
... ... @@ -164,12 +167,18 @@ GSkyDir* GSkyDir::clone(void) const
164 167 *
165 168 * @param[in] ra Right Ascension in radians.
166 169 * @param[in] dec Declination in radians.
  170 + *
  171 + * Sets Right Ascension and Declination in radians.
167 172 ***************************************************************************/
168 173 void GSkyDir::radec(const double& ra, const double& dec)
169 174 {
170 175 // Set attributes
171 176 m_has_lb = false;
172 177 m_has_radec = true;
  178 + #if defined(G_SINCOS_CACHE)
  179 + m_has_lb_cache = false;
  180 + m_has_radec_cache = false;
  181 + #endif
173 182  
174 183 // Set direction
175 184 m_ra = ra;
... ... @@ -185,12 +194,18 @@ void GSkyDir::radec(const double&amp; ra, const double&amp; dec)
185 194 *
186 195 * @param[in] ra Right Ascension in degrees.
187 196 * @param[in] dec Declination in degrees.
  197 + *
  198 + * Sets Right Ascension and Declination in degrees.
188 199 ***************************************************************************/
189 200 void GSkyDir::radec_deg(const double& ra, const double& dec)
190 201 {
191 202 // Set attributes
192 203 m_has_lb = false;
193 204 m_has_radec = true;
  205 + #if defined(G_SINCOS_CACHE)
  206 + m_has_lb_cache = false;
  207 + m_has_radec_cache = false;
  208 + #endif
194 209  
195 210 // Set direction
196 211 m_ra = ra * deg2rad;
... ... @@ -206,12 +221,18 @@ void GSkyDir::radec_deg(const double&amp; ra, const double&amp; dec)
206 221 *
207 222 * @param[in] l Galactic longitude in radians.
208 223 * @param[in] b Galactic latitude in radians.
  224 + *
  225 + * Sets Galactic longitude and latitude in radians.
209 226 ***************************************************************************/
210 227 void GSkyDir::lb(const double& l, const double& b)
211 228 {
212 229 // Set attributes
213 230 m_has_lb = true;
214 231 m_has_radec = false;
  232 + #if defined(G_SINCOS_CACHE)
  233 + m_has_lb_cache = false;
  234 + m_has_radec_cache = false;
  235 + #endif
215 236  
216 237 // Set direction
217 238 m_l = l;
... ... @@ -227,12 +248,18 @@ void GSkyDir::lb(const double&amp; l, const double&amp; b)
227 248 *
228 249 * @param[in] l Galactic longitude in degrees.
229 250 * @param[in] b Galactic latitude in degrees.
  251 + *
  252 + * Sets Galactic longitude and latitude in degrees.
230 253 ***************************************************************************/
231 254 void GSkyDir::lb_deg(const double& l, const double& b)
232 255 {
233 256 // Set attributes
234 257 m_has_lb = true;
235 258 m_has_radec = false;
  259 + #if defined(G_SINCOS_CACHE)
  260 + m_has_lb_cache = false;
  261 + m_has_radec_cache = false;
  262 + #endif
236 263  
237 264 // Set direction
238 265 m_l = l * deg2rad;
... ... @@ -247,12 +274,27 @@ void GSkyDir::lb_deg(const double&amp; l, const double&amp; b)
247 274 * @brief Set sky direction from 3D vector in celestial coordinates
248 275 *
249 276 * @param[in] vector 3D vector.
  277 + *
  278 + * Convert a 3-dimensional vector in celestial coordinates into a sky
  279 + * direction. The transformation is given by
  280 + *
  281 + * \f[
  282 + * \alpha = \atan \left( \frac{x_1}{x_0} \right)
  283 + * \f]
  284 + *
  285 + * \f[
  286 + * \delta = \asin x_2
  287 + * \f]
250 288 ***************************************************************************/
251 289 void GSkyDir::celvector(const GVector& vector)
252 290 {
253 291 // Set attributes
254 292 m_has_lb = false;
255 293 m_has_radec = true;
  294 + #if defined(G_SINCOS_CACHE)
  295 + m_has_lb_cache = false;
  296 + m_has_radec_cache = false;
  297 + #endif
256 298  
257 299 // Convert vector into sky position
258 300 m_dec = std::asin(vector[2]);
... ... @@ -277,8 +319,9 @@ void GSkyDir::celvector(const GVector&amp; vector)
277 319 void GSkyDir::rotate_deg(const double& phi, const double& theta)
278 320 {
279 321 // If we have no equatorial coordinates then get them now
280   - if (!m_has_radec && m_has_lb)
  322 + if (!m_has_radec && m_has_lb) {
281 323 gal2equ();
  324 + }
282 325  
283 326 // Allocate Euler and rotation matrices
284 327 GMatrix ry;
... ... @@ -311,28 +354,19 @@ void GSkyDir::rotate_deg(const double&amp; phi, const double&amp; theta)
311 354 }
312 355  
313 356  
314   -/***********************************************************************//**
315   - * @brief Returns galactic longitude in radians
316   - ***************************************************************************/
317   -double GSkyDir::l(void) const
318   -{
319   - // If we have no galactic coordinates then get them now
320   - if (!m_has_lb && m_has_radec)
321   - equ2gal();
322   -
323   - // Return galactic longitude
324   - return m_l;
325   -}
326 357  
327 358  
328 359 /***********************************************************************//**
329   - * @brief Returns galactic longitude in degrees
  360 + * @brief Return galactic longitude in degrees
  361 + *
  362 + * @return Galactic longitude in degrees.
330 363 ***************************************************************************/
331 364 double GSkyDir::l_deg(void) const
332 365 {
333 366 // If we have no galactic coordinates then get them now
334   - if (!m_has_lb && m_has_radec)
  367 + if (!m_has_lb && m_has_radec) {
335 368 equ2gal();
  369 + }
336 370  
337 371 // Return galactic longitude
338 372 return m_l * rad2deg;
... ... @@ -340,27 +374,16 @@ double GSkyDir::l_deg(void) const
340 374  
341 375  
342 376 /***********************************************************************//**
343   - * @brief Returns galactic latitude in radians
344   - ***************************************************************************/
345   -double GSkyDir::b(void) const
346   -{
347   - // If we have no galactic coordinates then get them now
348   - if (!m_has_lb && m_has_radec)
349   - equ2gal();
350   -
351   - // Return galactic latitude
352   - return m_b;
353   -}
354   -
355   -
356   -/***********************************************************************//**
357 377 * @brief Returns galactic latitude in degrees
  378 + *
  379 + * @return Galactic latitude in degrees.
358 380 ***************************************************************************/
359 381 double GSkyDir::b_deg(void) const
360 382 {
361 383 // If we have no galactic coordinates then get them now
362   - if (!m_has_lb && m_has_radec)
  384 + if (!m_has_lb && m_has_radec) {
363 385 equ2gal();
  386 + }
364 387  
365 388 // Return galactic latitude
366 389 return m_b * rad2deg;
... ... @@ -368,27 +391,16 @@ double GSkyDir::b_deg(void) const
368 391  
369 392  
370 393 /***********************************************************************//**
371   - * @brief Returns Right Ascension in radians
372   - ***************************************************************************/
373   -double GSkyDir::ra(void) const
374   -{
375   - // If we have no equatorial coordinates then get them now
376   - if (!m_has_radec && m_has_lb)
377   - gal2equ();
378   -
379   - // Return Right Ascension
380   - return m_ra;
381   -}
382   -
383   -
384   -/***********************************************************************//**
385 394 * @brief Returns Right Ascension in degrees
  395 + *
  396 + * @return Right Ascension in degrees.
386 397 ***************************************************************************/
387 398 double GSkyDir::ra_deg(void) const
388 399 {
389 400 // If we have no equatorial coordinates then get them now
390   - if (!m_has_radec && m_has_lb)
  401 + if (!m_has_radec && m_has_lb) {
391 402 gal2equ();
  403 + }
392 404  
393 405 // Return Right Ascension
394 406 return m_ra * rad2deg;
... ... @@ -396,27 +408,16 @@ double GSkyDir::ra_deg(void) const
396 408  
397 409  
398 410 /***********************************************************************//**
399   - * @brief Returns Declination in radians
400   - ***************************************************************************/
401   -double GSkyDir::dec(void) const
402   -{
403   - // If we have no equatorial coordinates then get them now
404   - if (!m_has_radec && m_has_lb)
405   - gal2equ();
406   -
407   - // Return Declination
408   - return m_dec;
409   -}
410   -
411   -
412   -/***********************************************************************//**
413 411 * @brief Returns Declination in degrees
  412 + *
  413 + * @return Declination in degrees.
414 414 ***************************************************************************/
415 415 double GSkyDir::dec_deg(void) const
416 416 {
417 417 // If we have no equatorial coordinates then get them now
418   - if (!m_has_radec && m_has_lb)
  418 + if (!m_has_radec && m_has_lb) {
419 419 gal2equ();
  420 + }
420 421  
421 422 // Return Declination
422 423 return m_dec * rad2deg;
... ... @@ -424,20 +425,32 @@ double GSkyDir::dec_deg(void) const
424 425  
425 426  
426 427 /***********************************************************************//**
427   - * @brief Returns sky direction as 3D vector in celestial coordinates
  428 + * @brief Return sky direction as 3D vector in celestial coordinates
  429 + *
  430 + * @return Sky direction as 3D vector in celestial coordinates.
428 431 ***************************************************************************/
429 432 GVector GSkyDir::celvector(void) const
430 433 {
431 434 // If we have no equatorial coordinates then get them now
432   - if (!m_has_radec && m_has_lb)
  435 + if (!m_has_radec && m_has_lb) {
433 436 gal2equ();
  437 + }
434 438  
435 439 // Compute 3D vector
436 440 double cosra = std::cos(m_ra);
437 441 double sinra = std::sin(m_ra);
  442 + #if defined(G_SINCOS_CACHE)
  443 + if (!m_has_radec_cache) {
  444 + m_has_radec_cache = true;
  445 + m_sin_dec = std::sin(m_dec);
  446 + m_cos_dec = std::cos(m_dec);
  447 + }
  448 + GVector vector(m_cos_dec*cosra, m_cos_dec*sinra, m_sin_dec);
  449 + #else
438 450 double cosdec = std::cos(m_dec);
439 451 double sindec = std::sin(m_dec);
440 452 GVector vector(cosdec*cosra, cosdec*sinra, sindec);
  453 + #endif
441 454  
442 455 // Return vector
443 456 return vector;
... ... @@ -448,6 +461,9 @@ GVector GSkyDir::celvector(void) const
448 461 * @brief Compute angular distance between sky directions in radians
449 462 *
450 463 * @param[in] dir Sky direction to which distance is to be computed.
  464 + * @return Angular distance in radians.
  465 + *
  466 + * Computes the angular distance between two sky directions in radians.
451 467 ***************************************************************************/
452 468 double GSkyDir::dist(const GSkyDir& dir) const
453 469 {
... ... @@ -456,25 +472,82 @@ double GSkyDir::dist(const GSkyDir&amp; dir) const
456 472  
457 473 // Compute dependent on coordinate system availability. This speeds
458 474 // up things by avoiding unnecessary coordinate transformations.
459   - if (m_has_lb && dir.m_has_lb) {
460   - cosdis = std::sin(m_b)*std::sin(dir.m_b) +
461   - std::cos(m_b)*std::cos(dir.m_b) * std::cos(dir.m_l - m_l);
462   - }
463   - else if (m_has_radec && dir.m_has_radec) {
464   - cosdis = std::sin(m_dec)*std::sin(dir.m_dec) +
465   - std::cos(m_dec)*std::cos(dir.m_dec) * std::cos(dir.m_ra - m_ra);
466   - }
467   - else if (m_has_lb) {
468   - cosdis = std::sin(m_b)*std::sin(dir.b()) +
469   - std::cos(m_b)*std::cos(dir.b()) * std::cos(dir.l() - m_l);
  475 + if (m_has_lb) {
  476 + #if defined(G_SINCOS_CACHE)
  477 + if (!m_has_lb_cache) {
  478 + m_has_lb_cache = true;
  479 + m_sin_b = std::sin(m_b);
  480 + m_cos_b = std::cos(m_b);
  481 + }
  482 + #endif
  483 + if (dir.m_has_lb) {
  484 + #if defined(G_SINCOS_CACHE)
  485 + if (!dir.m_has_lb_cache) {
  486 + dir.m_has_lb_cache = true;
  487 + dir.m_sin_b = std::sin(dir.m_b);
  488 + dir.m_cos_b = std::cos(dir.m_b);
  489 + }
  490 + cosdis = m_sin_b * dir.m_sin_b +
  491 + m_cos_b * dir.m_cos_b *
  492 + std::cos(dir.m_l - m_l);
  493 + #else
  494 + cosdis = std::sin(m_b) * std::sin(dir.m_b) +
  495 + std::cos(m_b) * std::cos(dir.m_b) *
  496 + std::cos(dir.m_l - m_l);
  497 + #endif
  498 + }
  499 + else {
  500 + #if defined(G_SINCOS_CACHE)
  501 + cosdis = m_sin_b * std::sin(dir.b()) +
  502 + m_cos_b * std::cos(dir.b()) *
  503 + std::cos(dir.l() - m_l);
  504 + #else
  505 + cosdis = std::sin(m_b) * std::sin(dir.b()) +
  506 + std::cos(m_b) * std::cos(dir.b()) *
  507 + std::cos(dir.l() - m_l);
  508 + #endif
  509 + }
470 510 }
471 511 else if (m_has_radec) {
472   - cosdis = sin(m_dec)*sin(dir.dec()) +
473   - cos(m_dec)*cos(dir.dec()) * std::cos(dir.ra() - m_ra);
  512 + #if defined(G_SINCOS_CACHE)
  513 + if (!m_has_radec_cache) {
  514 + m_has_radec_cache = true;
  515 + m_sin_dec = std::sin(m_dec);
  516 + m_cos_dec = std::cos(m_dec);
  517 + }
  518 + #endif
  519 + if (dir.m_has_radec) {
  520 + #if defined(G_SINCOS_CACHE)
  521 + if (!dir.m_has_radec_cache) {
  522 + dir.m_has_radec_cache = true;
  523 + dir.m_sin_dec = std::sin(dir.m_dec);
  524 + dir.m_cos_dec = std::cos(dir.m_dec);
  525 + }
  526 + cosdis = m_sin_dec * dir.m_sin_dec +
  527 + m_cos_dec * dir.m_cos_dec *
  528 + std::cos(dir.m_ra - m_ra);
  529 + #else
  530 + cosdis = std::sin(m_dec) * std::sin(dir.m_dec) +
  531 + std::cos(m_dec) * std::cos(dir.m_dec) *
  532 + std::cos(dir.m_ra - m_ra);
  533 + #endif
  534 + }
  535 + else {
  536 + #if defined(G_SINCOS_CACHE)
  537 + cosdis = m_sin_dec * sin(dir.dec()) +
  538 + m_cos_dec * cos(dir.dec()) *
  539 + std::cos(dir.ra() - m_ra);
  540 + #else
  541 + cosdis = sin(m_dec) * sin(dir.dec()) +
  542 + cos(m_dec) * cos(dir.dec()) *
  543 + std::cos(dir.ra() - m_ra);
  544 + #endif
  545 + }
474 546 }
475 547 else {
476   - cosdis = std::sin(dec())*std::sin(dir.dec()) +
477   - std::cos(dec())*std::cos(dir.dec()) * std::cos(dir.ra() - ra());
  548 + cosdis = std::sin(dec()) * std::sin(dir.dec()) +
  549 + std::cos(dec()) * std::cos(dir.dec()) *
  550 + std::cos(dir.ra() - ra());
478 551 }
479 552  
480 553 // Compute distance (use argument save GTools function)
... ... @@ -489,6 +562,9 @@ double GSkyDir::dist(const GSkyDir&amp; dir) const
489 562 * @brief Compute angular distance between sky directions in degrees
490 563 *
491 564 * @param[in] dir Sky direction to which distance is to be computed.
  565 + * @return Angular distance in degrees.
  566 + *
  567 + * Computes the angular distance between two sky directions in degrees.
492 568 ***************************************************************************/
493 569 double GSkyDir::dist_deg(const GSkyDir& dir) const
494 570 {
... ... @@ -501,6 +577,7 @@ double GSkyDir::dist_deg(const GSkyDir&amp; dir) const
501 577 * @brief Compute position angle between sky directions in radians
502 578 *
503 579 * @param[in] dir Sky direction.
  580 + * @return Position angle in radians.
504 581 *
505 582 * Computes the position angle using
506 583 * \f[PA = \arctan \left(
... ... @@ -522,25 +599,68 @@ double GSkyDir::posang(const GSkyDir&amp; dir) const
522 599  
523 600 // Compute dependent on coordinate system availability. This speeds
524 601 // up things by avoiding unnecessary coordinate transformations.
525   - if (m_has_lb && dir.m_has_lb) {
526   - arg_1 = std::sin(dir.m_l - m_l);
527   - arg_2 = std::cos(m_b)*std::tan(dir.m_b) - std::sin(m_b)*std::cos(dir.m_l - m_l);
528   - }
529   - else if (m_has_radec && dir.m_has_radec) {
530   - arg_1 = std::sin(dir.m_ra - m_ra);
531   - arg_2 = std::cos(m_dec)*std::tan(dir.m_dec) - std::sin(m_dec)*std::cos(dir.m_ra - m_ra);
532   - }
533   - else if (m_has_lb) {
534   - arg_1 = std::sin(dir.l() - m_l);
535   - arg_2 = std::cos(m_b)*std::tan(dir.b()) - std::sin(m_b)*std::cos(dir.l() - m_l);
  602 + if (m_has_lb) {
  603 + #if defined(G_SINCOS_CACHE)
  604 + if (!m_has_lb_cache) {
  605 + m_has_lb_cache = true;
  606 + m_sin_b = std::sin(m_b);
  607 + m_cos_b = std::cos(m_b);
  608 + }
  609 + #endif
  610 + if (dir.m_has_lb) {
  611 + arg_1 = std::sin(dir.m_l - m_l);
  612 + #if defined(G_SINCOS_CACHE)
  613 + arg_2 = m_cos_b * std::tan(dir.m_b) -
  614 + m_sin_b * std::cos(dir.m_l - m_l);
  615 + #else
  616 + arg_2 = std::cos(m_b) * std::tan(dir.m_b) -
  617 + std::sin(m_b) * std::cos(dir.m_l - m_l);
  618 + #endif
  619 + }
  620 + else {
  621 + arg_1 = std::sin(dir.l() - m_l);
  622 + #if defined(G_SINCOS_CACHE)
  623 + arg_2 = m_cos_b * std::tan(dir.b()) -
  624 + m_sin_b * std::cos(dir.l() - m_l);
  625 + #else
  626 + arg_2 = std::cos(m_b) * std::tan(dir.b()) -
  627 + std::sin(m_b) * std::cos(dir.l() - m_l);
  628 + #endif
  629 + }
536 630 }
537 631 else if (m_has_radec) {
538   - arg_1 = std::sin(dir.ra() - m_ra);
539   - arg_2 = std::cos(m_dec)*std::tan(dir.dec()) - std::sin(m_dec)*std::cos(dir.ra() - m_ra);
  632 + #if defined(G_SINCOS_CACHE)
  633 + if (!m_has_radec_cache) {
  634 + m_has_radec_cache = true;
  635 + m_sin_dec = std::sin(m_dec);
  636 + m_cos_dec = std::cos(m_dec);
  637 + }
  638 + #endif
  639 + if (dir.m_has_radec) {
  640 + arg_1 = std::sin(dir.m_ra - m_ra);
  641 + #if defined(G_SINCOS_CACHE)
  642 + arg_2 = m_cos_dec * std::tan(dir.m_dec) -
  643 + m_sin_dec * std::cos(dir.m_ra - m_ra);
  644 + #else
  645 + arg_2 = std::cos(m_dec) * std::tan(dir.m_dec) -
  646 + std::sin(m_dec) * std::cos(dir.m_ra - m_ra);
  647 + #endif
  648 + }
  649 + else {
  650 + arg_1 = std::sin(dir.ra() - m_ra);
  651 + #if defined(G_SINCOS_CACHE)
  652 + arg_2 = m_cos_dec * std::tan(dir.dec()) -
  653 + m_sin_dec * std::cos(dir.ra() - m_ra);
  654 + #else
  655 + arg_2 = std::cos(m_dec) * std::tan(dir.dec()) -
  656 + std::sin(m_dec) * std::cos(dir.ra() - m_ra);
  657 + #endif
  658 + }
540 659 }
541 660 else {
542 661 arg_1 = std::sin(dir.ra() - ra());
543   - arg_2 = std::cos(dec())*std::tan(dir.dec()) - std::sin(dec())*std::cos(dir.ra() - ra());
  662 + arg_2 = std::cos(dec())*std::tan(dir.dec()) -
  663 + std::sin(dec())*std::cos(dir.ra() - ra());
544 664 }
545 665  
546 666 // Compute position angle
... ... @@ -555,8 +675,11 @@ double GSkyDir::posang(const GSkyDir&amp; dir) const
555 675 * @brief Compute position angle between sky directions in degrees
556 676 *
557 677 * @param[in] dir Sky direction.
  678 + * @return Position angle in degrees.
  679 + *
  680 + * Returns the position angle of a sky direction in degrees.
558 681 *
559   - * See GSkyDir::posang for more information about the computation of the
  682 + * See the posang() method for more information about the computation of the
560 683 * position angle.
561 684 ***************************************************************************/
562 685 double GSkyDir::posang_deg(const GSkyDir& dir) const
... ... @@ -567,7 +690,9 @@ double GSkyDir::posang_deg(const GSkyDir&amp; dir) const
567 690  
568 691  
569 692 /***********************************************************************//**
570   - * @brief Print vector information
  693 + * @brief Print sky direction information
  694 + *
  695 + * @return String containing sky direction information.
571 696 ***************************************************************************/
572 697 std::string GSkyDir::print(void) const
573 698 {
... ... @@ -575,12 +700,15 @@ std::string GSkyDir::print(void) const
575 700 std::string result;
576 701  
577 702 // Put coordinates in string
578   - if (m_has_lb)
  703 + if (m_has_lb) {
579 704 result = "(l,b)=("+str(m_l*rad2deg)+","+str(m_b*rad2deg)+")";
580   - else if (m_has_radec)
  705 + }
  706 + else if (m_has_radec) {
581 707 result = "(RA,Dec)=("+str(m_ra*rad2deg)+","+str(m_dec*rad2deg)+")";
582   - else
  708 + }
  709 + else {
583 710 result = "(RA,Dec)=(not initialised)";
  711 + }
584 712  
585 713 // Return result
586 714 return result;
... ... @@ -606,6 +734,16 @@ void GSkyDir::init_members(void)
606 734 m_ra = 0.0;
607 735 m_dec = 0.0;
608 736  
  737 + // Initialise sincos cache
  738 + #if defined(G_SINCOS_CACHE)
  739 + m_has_lb_cache = false;
  740 + m_has_radec_cache = false;
  741 + m_sin_b = 0.0;
  742 + m_cos_b = 0.0;
  743 + m_sin_dec = 0.0;
  744 + m_cos_dec = 0.0;
  745 + #endif
  746 +
609 747 // Return
610 748 return;
611 749 }
... ... @@ -626,6 +764,16 @@ void GSkyDir::copy_members(const GSkyDir&amp; dir)
626 764 m_ra = dir.m_ra;
627 765 m_dec = dir.m_dec;
628 766  
  767 + // Copy sincos cache
  768 + #if defined(G_SINCOS_CACHE)
  769 + m_has_lb_cache = dir.m_has_lb_cache;
  770 + m_has_radec_cache = dir.m_has_radec_cache;
  771 + m_sin_b = dir.m_sin_b;
  772 + m_cos_b = dir.m_cos_b;
  773 + m_sin_dec = dir.m_sin_dec;
  774 + m_cos_dec = dir.m_cos_dec;
  775 + #endif
  776 +
629 777 // Return
630 778 return;
631 779 }
... ... @@ -744,35 +892,49 @@ bool operator==(const GSkyDir &amp;a, const GSkyDir &amp;b)
744 892  
745 893 // Compute dependent on coordinate system availability. This speeds
746 894 // up things by avoiding unnecessary coordinate transformations.
747   - if (a.m_has_lb && b.m_has_lb) {
748   - if (std::abs(a.m_b) == 90.0)
749   - equal = (a.m_b == b.m_b);
750   - else
751   - equal = (a.m_b == b.m_b && a.m_l == b.m_l);
752   - }
753   - else if (a.m_has_radec && b.m_has_radec) {
754   - if (std::abs(a.m_dec) == 90.0)
755   - equal = (a.m_dec == b.m_dec);
756   - else
757   - equal = (a.m_dec == b.m_dec && a.m_ra == b.m_ra);
758   - }
759   - else if (a.m_has_lb) {
760   - if (std::abs(a.m_b) == 90.0)
761   - equal = (a.m_b == b.b());
762   - else
763   - equal = (a.m_b == b.b() && a.m_l == b.l());
  895 + if (a.m_has_lb) {
  896 + if (b.m_has_lb) {
  897 + if (std::abs(a.m_b) == 90.0) {
  898 + equal = (a.m_b == b.m_b);
  899 + }
  900 + else {
  901 + equal = (a.m_b == b.m_b && a.m_l == b.m_l);
  902 + }
  903 + }
  904 + else {
  905 + if (std::abs(a.m_b) == 90.0) {
  906 + equal = (a.m_b == b.b());
  907 + }
  908 + else {
  909 + equal = (a.m_b == b.b() && a.m_l == b.l());
  910 + }
  911 + }
764 912 }
765 913 else if (a.m_has_radec) {
766   - if (std::abs(a.m_dec) == 90.0)
  914 + if (b.m_has_radec) {
  915 + if (std::abs(a.m_dec) == 90.0) {
  916 + equal = (a.m_dec == b.m_dec);
  917 + }
  918 + else {
  919 + equal = (a.m_dec == b.m_dec && a.m_ra == b.m_ra);
  920 + }
  921 + }
  922 + else {
  923 + if (std::abs(a.m_dec) == 90.0) {
767 924 equal = (a.m_dec == b.dec());
768   - else
769   - equal = (a.m_dec == b.dec() && a.m_ra == b.ra());
  925 + }
  926 + else {
  927 + equal = (a.m_dec == b.dec() && a.m_ra == b.ra());
  928 + }
  929 + }
770 930 }
771 931 else {
772   - if (std::abs(b.dec()) == 90.0)
  932 + if (std::abs(b.dec()) == 90.0) {
773 933 equal = (b.dec() == a.dec());
774   - else
  934 + }
  935 + else {
775 936 equal = (b.dec() == a.dec() && b.ra() == a.ra());
  937 + }
776 938 }
777 939  
778 940 // Return equality
... ...