Pierrick Martin / gammalib

Commit 964b56a964c0a5e28bf3d3f8cbd1cfcdf5b4c77a

Authored by Jürgen Knödlseder
1 parent b0fb9d94

Fix ambiguity problem in GSkyDir::posang() and posang_deg() (#3199)

Inline Side-by-side
Showing 15 changed files with 165 additions and 111 deletions
ChangeLog
  View file @964b56a
1   -2020-02-17
  1 +2020-04-13
2 2  
3 3 * Version 1.7.0 released
4 4 ========================
5 5  
  6 + Fix ambiguity problem in GSkyDir::posang() and posang_deg() (#3199)
6 7 Avoid rounding errors in GSkyMap operators (#3126)
7 8 Add GObservations::npred() methods (#3062)
8 9 Exclude upper edge boundary in GSkyMap::contains() method (#3012)
... ...
1 1 New Features and Important Changes in GammaLib 1.7.0
2 2  
3   -11 April 2020
  3 +13 April 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -56,6 +56,8 @@ The following methods have been renamed:
56 56  
57 57 The arguments for the following methods have been changed:
58 58 - GTime::utc(int&) - Add optional precision value
  59 +- GSkyDir::posang(GSkyDir&, std::string&) - Add optional coordinate system string
  60 +- GSkyDir::posang_dir(GSkyDir&, std::string&) - Add optional coordinate system string
59 61  
60 62 The return value of the following methods has been changed:
61 63 - none
... ... @@ -142,6 +144,11 @@ None
142 144  
143 145 12. Sky module
144 146 --------------
  147 +The position angle is now computed by GSkyDir::posang() and GSkyDir::posang_deg()
  148 +for the coordinate system that is specified by the coordsys argument. By default
  149 +this argument is set to "CEL", hence by default the methods return the position
  150 +angle with respect to celestial coordinates (#3199).
  151 +
145 152 The GSkyMap operators used the sky map value access operator that bi-linearly
146 153 interpolates sky map values even for maps with identical definition. This could
147 154 lead to rounding errors. For maps with identical definition the operators now
... ...
doc/source/admin/release_history/1.7.rst
  View file @964b56a
... ... @@ -12,12 +12,14 @@ In particular, this release provides:
12 12  
13 13 * Support for INTEGRAL/SPI data analysis
14 14 * Response caching for all instruments
15   -* Support for extended models for all instruments
  15 +* ...
16 16  
17 17  
18 18 Bug fixes
19 19 ---------
20 20  
  21 +* [`3199 <https://cta-redmine.irap.omp.eu/issues/3199>`_] -
  22 + Fix ambiguity problem in ``GSkyDir::posang()`` and ``GSkyDir::posang_deg()``
21 23 * [`3126 <https://cta-redmine.irap.omp.eu/issues/3126>`_] -
22 24 Avoid rounding errors in ``GSkyMap`` operators.
23 25 * [`3012 <https://cta-redmine.irap.omp.eu/issues/3012>`_] -
... ...
include/GSkyDir.hpp
  View file @964b56a
... ... @@ -87,8 +87,8 @@ public:
87 87 void rotate(const double& phi, const double& theta);
88 88 void rotate_deg(const double& phi, const double& theta);
89 89 void precess(const double& from_epoch, const double& to_epoch);
90   - void sun(const GTime& time, const double& epoch = 2000.0);
91   - void moon(const GTime& time, const double& epoch = 2000.0);
  90 + void sun(const GTime& time, const double& epoch = 2000.0);
  91 + void moon(const GTime& time, const double& epoch = 2000.0);
92 92 const double& l(void) const;
93 93 const double& b(void) const;
94 94 const double& ra(void) const;
... ... @@ -101,8 +101,10 @@ public:
101 101 double cos_dist(const GSkyDir& dir) const;
102 102 double dist(const GSkyDir& dir) const;
103 103 double dist_deg(const GSkyDir& dir) const;
104   - double posang(const GSkyDir& dir) const;
105   - double posang_deg(const GSkyDir& dir) const;
  104 + double posang(const GSkyDir& dir,
  105 + const std::string& coordsys = "CEL") const;
  106 + double posang_deg(const GSkyDir& dir,
  107 + const std::string& coordsys = "CEL") const;
106 108 std::string print(const GChatter& chatter = NORMAL) const;
107 109  
108 110 private:
... ... @@ -289,17 +291,23 @@ double GSkyDir::dist_deg(const GSkyDir&amp; dir) const
289 291 * @brief Compute position angle between sky directions in degrees
290 292 *
291 293 * @param[in] dir Sky direction.
292   - * @return Position angle in degrees.
  294 + * @param[in] coordsys Coordinate system ("CEL" or "GAL")
  295 + * @return Position angle (deg).
293 296 *
294   - * Returns the position angle of a sky direction in degrees.
  297 + * Computes the position angle of a specified sky direction with respect to
  298 + * the sky direction of the instance in degrees. If "CEL" is specified as
  299 + * @p coordsys the position angle is counted counterclockwise from celestial
  300 + * North, if "GAL" is specified the position angle is counted
  301 + * counterclockwise from Galactic North.
295 302 *
296 303 * See the posang() method for more information about the computation of the
297 304 * position angle.
298 305 ***************************************************************************/
299 306 inline
300   -double GSkyDir::posang_deg(const GSkyDir& dir) const
  307 +double GSkyDir::posang_deg(const GSkyDir& dir,
  308 + const std::string& coordsys) const
301 309 {
302   - return (posang(dir) * gammalib::rad2deg);
  310 + return (posang(dir, coordsys) * gammalib::rad2deg);
303 311 }
304 312  
305 313 #endif /* GSKYDIR_HPP */
... ...
inst/com/include/GCOMOad.hpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GCOMOad.hpp - COMPTEL Orbit Aspect Data class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017 by Juergen Knodlseder *
  4 + * copyright (C) 2017-2020 by Juergen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -409,12 +409,11 @@ double GCOMOad::theta(const GSkyDir&amp; sky) const
409 409 inline
410 410 double GCOMOad::phi(const GSkyDir& sky) const
411 411 {
412   - // If position angle has not be initialised the do it now
  412 + // If position angle has not be initialised then do it now
413 413 if (m_posang > 1.0e20) {
414   - m_posang = m_zaxis.posang_deg(m_xaxis);
  414 + m_posang = m_zaxis.posang_deg(m_xaxis); // Celestial system
415 415 }
416   - //return (m_zaxis.posang_deg(m_xaxis) - m_zaxis.posang_deg(sky));
417   - return (m_posang - m_zaxis.posang_deg(sky));
  416 + return (m_posang - m_zaxis.posang_deg(sky)); // Celestial system
418 417 }
419 418  
420 419 #endif /* GCOMOAD_HPP */
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @964b56a
... ... @@ -1582,7 +1582,7 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation&amp; obs,
1582 1582  
1583 1583 // Compute position of pixel centre in instrument coordinates
1584 1584 double theta = obsdir.dist(pixdir);
1585   - double phi = obsdir.posang(pixdir);
  1585 + double phi = obsdir.posang(pixdir); // Celestial system
1586 1586  
1587 1587 // Add up effective area
1588 1588 m_arf[i] += rsp.aeff(theta, phi,
... ... @@ -2068,7 +2068,7 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation&amp; obs,
2068 2068  
2069 2069 // Compute position of pixel centre in instrument coordinates
2070 2070 double theta = obsdir.dist(pixdir);
2071   - double phi = obsdir.posang(pixdir);
  2071 + double phi = obsdir.posang(pixdir); // Celestial system
2072 2072  
2073 2073 // Loop over true energy
2074 2074 for (int itrue = 0; itrue < ntrue; ++itrue) {
... ...
inst/cta/src/GCTAPointing.cpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GCTAPointing.cpp - CTA pointing class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -262,7 +262,7 @@ GCTAInstDir GCTAPointing::instdir(const GSkyDir&amp; skydir) const
262 262 }
263 263 #else
264 264 double theta = m_dir.dist(skydir);
265   - double phi = m_dir.posang(skydir);
  265 + double phi = m_dir.posang(skydir); // Celestial system
266 266 double detx(0.0);
267 267 double dety(0.0);
268 268 if (theta > 0.0 ) {
... ...
inst/cta/src/GCTAResponseCube.cpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GCTAResponseCube.cpp - CTA cube analysis response function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1407,7 +1407,7 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp; event,
1407 1407 // Compute angle between model centre and measured photon direction and
1408 1408 // position angle (radians)
1409 1409 double rho_obs = model->dir().dist(obsDir);
1410   - double posangle_obs = model->dir().posang(obsDir);
  1410 + double posangle_obs = model->dir().posang(obsDir); // Celestial
1411 1411  
1412 1412 // Get livetime (in seconds)
1413 1413 double livetime = exposure().livetime();
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GCTAResponseIrf.cpp - CTA instrument response function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -2519,12 +2519,12 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp; event,
2519 2519 // centre and position angle of observed photon direction seen from the
2520 2520 // model centre [radians]
2521 2521 double rho_obs = centre.dist(obsDir);
2522   - double posangle_obs = centre.posang(obsDir);
  2522 + double posangle_obs = centre.posang(obsDir); // Celestial system
2523 2523  
2524 2524 // Determine angular distance between model centre and pointing direction
2525 2525 // [radians]
2526 2526 double rho_pnt = centre.dist(pnt.dir());
2527   - double posangle_pnt = centre.posang(pnt.dir());
  2527 + double posangle_pnt = centre.posang(pnt.dir()); // Celestial system
2528 2528  
2529 2529 // Compute azimuth angle of pointing in model coordinate system [radians]
2530 2530 double omega_pnt = posangle_pnt - posangle_obs;
... ... @@ -2957,7 +2957,7 @@ double GCTAResponseIrf::nroi_radial(const GModelSky&amp; model,
2957 2957  
2958 2958 // Compute position angle of ROI centre with respect to model
2959 2959 // centre (radians)
2960   - double omega0 = centre.posang(roi.centre().dir());
  2960 + double omega0 = centre.posang(roi.centre().dir()); // Celestial system
2961 2961  
2962 2962 // Setup integration kernel
2963 2963 cta_nroi_radial_kern_rho integrand(this,
... ... @@ -3162,7 +3162,7 @@ double GCTAResponseIrf::nroi_elliptical(const GModelSky&amp; model,
3162 3162  
3163 3163 // Compute position angle of ROI centre with respect to model
3164 3164 // centre (radians)
3165   - double posangle_roi = centre.posang(roi.centre().dir());
  3165 + double posangle_roi = centre.posang(roi.centre().dir()); // Celestial system
3166 3166  
3167 3167 // Setup integration kernel
3168 3168 cta_nroi_elliptical_kern_rho integrand(this,
... ...
inst/spi/include/GSPIResponse.hpp
  View file @964b56a
... ... @@ -130,7 +130,7 @@ private:
130 130  
131 131 // Private cache
132 132 std::vector<GSkyDir> m_spix; //!< SPI pointing direction
133   - std::vector<double> m_posang; //!< Position angle of Y axis (radians)
  133 + std::vector<double> m_posang; //!< Position angle of Y axis (CEL, radians)
134 134 mutable bool m_has_wcs; //!< Has WCS information
135 135 mutable double m_wcs_xmin; //!< Minimum X value (radians)
136 136 mutable double m_wcs_ymin; //!< Minimum Y value (radians)
... ... @@ -241,7 +241,7 @@ double GSPIResponse::zenith(const int&amp; ipt, const GSkyDir&amp; dir) const
241 241 inline
242 242 double GSPIResponse::azimuth(const int& ipt, const GSkyDir& dir) const
243 243 {
244   - double azimuth = m_posang[ipt] - m_spix[ipt].posang(dir);
  244 + double azimuth = m_posang[ipt] - m_spix[ipt].posang(dir); // Celestial system
245 245 if (azimuth < 0.0) {
246 246 azimuth += gammalib::twopi;
247 247 }
... ...
inst/spi/src/GSPIResponse.cpp
  View file @964b56a
... ... @@ -727,10 +727,6 @@ void GSPIResponse::init_members(void)
727 727 m_wcs_ypix_max = 0.0;
728 728 m_max_zenith = 180.0 * gammalib::deg2rad;
729 729  
730   - // Set energy scale for response cache to MeV
731   - //m_irf_cache.energy_scale(GEnergy(1.0, "MeV"));
732   - //m_nroi_cache.energy_scale(GEnergy(1.0, "MeV"));
733   -
734 730 // Return
735 731 return;
736 732 }
... ... @@ -1287,7 +1283,7 @@ void GSPIResponse::set_cache(const GSPIEventCube* cube)
1287 1283 // Loop over all pointings
1288 1284 for (int ipt = 0; ipt < npt; ++ipt) {
1289 1285  
1290   - // Compute position angle
  1286 + // Compute position angle is celestial coordinates
1291 1287 double posang = cube->spi_x(ipt).posang(cube->spi_z(ipt)) + gammalib::pihalf;
1292 1288  
1293 1289 // Store angles
... ...
pyext/GSkyDir.i
  View file @964b56a
... ... @@ -55,8 +55,8 @@ public:
55 55 void rotate(const double& phi, const double& theta);
56 56 void rotate_deg(const double& phi, const double& theta);
57 57 void precess(const double& from_epoch, const double& to_epoch);
58   - void sun(const GTime& time, const double& epoch = gammalib::epoch2000);
59   - void moon(const GTime& time, const double& epoch = gammalib::epoch2000);
  58 + void sun(const GTime& time, const double& epoch = 2000.0);
  59 + void moon(const GTime& time, const double& epoch = 2000.0);
60 60 const double& l(void) const;
61 61 const double& b(void) const;
62 62 const double& ra(void) const;
... ... @@ -69,8 +69,10 @@ public:
69 69 double cos_dist(const GSkyDir& dir) const;
70 70 double dist(const GSkyDir& dir) const;
71 71 double dist_deg(const GSkyDir& dir) const;
72   - double posang(const GSkyDir& dir) const;
73   - double posang_deg(const GSkyDir& dir) const;
  72 + double posang(const GSkyDir& dir,
  73 + const std::string& coordsys = "CEL") const;
  74 + double posang_deg(const GSkyDir& dir,
  75 + const std::string& coordsys = "CEL") const;
74 76 };
75 77  
76 78  
... ...
src/model/GModelSpatialElliptical.cpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GModelSpatialElliptical.cpp - Abstract elliptical spatial model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -176,7 +176,7 @@ double GModelSpatialElliptical::eval(const GPhoton&amp; photon,
176 176 // Compute distance from source and position angle (in radians)
177 177 const GSkyDir& srcDir = photon.dir();
178 178 double theta = dir().dist(srcDir);
179   - double posang = dir().posang(srcDir);
  179 + double posang = dir().posang(srcDir); // Celestial system
180 180  
181 181 // Evaluate model
182 182 double value = eval(theta, posang, photon.energy(), photon.time(), gradients);
... ...
src/sky/GSkyDir.cpp
  View file @964b56a
1 1 /***************************************************************************
2 2 * GSkyDir.cpp - Sky direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -38,6 +38,7 @@
38 38 #include "GVector.hpp"
39 39  
40 40 /* __ Method name definitions ____________________________________________ */
  41 +#define G_POSANG "GSkyDir::posang(GSkyDir&, std::string&)"
41 42  
42 43 /* __ Macros _____________________________________________________________ */
43 44  
... ... @@ -507,7 +508,7 @@ void GSkyDir::sun(const GTime&amp; time, const double&amp; epoch)
507 508 while (mean_longitude >= 360.0) {
508 509 mean_longitude -= 360.0;
509 510 }
510   -
  511 +
511 512 // Compute the mean anomaly of the Sun in radians
512 513 double mean_anomaly = (357.528 + 0.9856003 * n) * gammalib::deg2rad;
513 514  
... ... @@ -872,7 +873,7 @@ void GSkyDir::moon(const GTime&amp; time, const double&amp; epoch)
872 873  
873 874 // Add nutation in longitude
874 875 geolong += nut_long / 3600.0;
875   -
  876 +
876 877 // nutate, jd, nlong, elong
877 878 // geolong= geolong + nlong/3.6d3
878 879 // cirrange,geolong
... ... @@ -1049,90 +1050,96 @@ double GSkyDir::cos_dist(const GSkyDir&amp; dir) const
1049 1050 * @brief Compute position angle between sky directions in radians
1050 1051 *
1051 1052 * @param[in] dir Sky direction.
1052   - * @return Position angle in radians.
  1053 + * @param[in] coordsys Coordinate system ("CEL" or "GAL")
  1054 + * @return Position angle (rad).
  1055 + *
  1056 + * Computes the position angle of a specified sky direction with respect to
  1057 + * the sky direction of the instance in radians. If "CEL" is specified as
  1058 + * @p coordsys the position angle is counted counterclockwise from celestial
  1059 + * North, if "GAL" is specified the position angle is counted
  1060 + * counterclockwise from Galactic North.
  1061 + *
  1062 + * If @p coordsys is "CEL" the position angle is computed using
1053 1063 *
1054   - * Computes the position angle using
1055   - * \f[PA = \arctan \left(
1056   - * \frac{\sin( \alpha_1 - \alpha_0 )}
1057   - * {\cos \delta_0 \tan \delta_1 -
1058   - * \sin \delta_0 \cos(\alpha_1 - \alpha_0)} \right)\f]
  1064 + * \f[PA = \arctan \left(
  1065 + * \frac{\sin( \alpha - \alpha_0 )}
  1066 + * {\cos \delta_0 \tan \delta -
  1067 + * \sin \delta_0 \cos(\alpha - \alpha_0)} \right)\f]
1059 1068 * where
1060   - * \f$(\alpha_0,\delta_0)\f$ are the coordinates of the reference point, and
1061   - * \f$(\alpha_1,\delta_1)\f$ are the coordinates of sky direction for which
1062   - * the position angle is to be computed.
  1069 + * * \f$(\alpha_0,\delta_0)\f$ are the celestial coordinates of the instance,
  1070 + * and
  1071 + * * \f$(\alpha,\delta)\f$ are the celestial coordinates of sky direction.
1063 1072 *
1064   - * The position angle is counted counterclockwise from north.
  1073 + * If @p coordsys is "GAL" the position angle is computed using
  1074 + *
  1075 + * \f[PA = \arctan \left(
  1076 + * \frac{\sin( \l - \l_0 )}
  1077 + * {\cos \b_0 \tan \b - \sin \b_0 \cos(\l - \l_0)} \right)\f]
  1078 + * where
  1079 + * * \f$(\l_0,\l_0)\f$ are the Galactic coordinates of the instance,
  1080 + * and
  1081 + * * \f$(\l,\b)\f$ are the Galactic coordinates of sky direction.
1065 1082 ***************************************************************************/
1066   -double GSkyDir::posang(const GSkyDir& dir) const
  1083 +double GSkyDir::posang(const GSkyDir& dir, const std::string& coordsys) const
1067 1084 {
1068 1085 // Initialise arguments of arctan
1069 1086 double arg_1;
1070 1087 double arg_2;
1071 1088  
1072   - // Compute dependent on coordinate system availability. This speeds
1073   - // up things by avoiding unnecessary coordinate transformations.
1074   - if (m_has_lb) {
  1089 + // Case A: Computation in celestial coordinates
  1090 + if (coordsys == "CEL") {
  1091 +
  1092 + // If required then compute sin(dec) and cos(dec) cache values
  1093 + #if defined(G_SINCOS_CACHE)
  1094 + if (!m_has_radec_cache) {
  1095 + m_has_radec_cache = true;
  1096 + m_sin_dec = std::sin(dec());
  1097 + m_cos_dec = std::cos(dec());
  1098 + }
  1099 + #endif
  1100 +
  1101 + // Compute arguments of arctan
  1102 + arg_1 = std::sin(dir.ra() - ra());
  1103 + #if defined(G_SINCOS_CACHE)
  1104 + arg_2 = m_cos_dec * std::tan(dir.dec()) -
  1105 + m_sin_dec * std::cos(dir.ra() - ra());
  1106 + #else
  1107 + arg_2 = std::cos(dec()) * std::tan(dir.dec()) -
  1108 + std::sin(dec()) * std::cos(dir.ra() - ra());
  1109 + #endif
  1110 +
  1111 + } // endif: computation in celestial coordinates
  1112 +
  1113 + // Case B: Computation in Galactic coordinates
  1114 + else if (coordsys == "GAL") {
  1115 +
  1116 + // If required then compute sin(b) and cos(b) cache values
1075 1117 #if defined(G_SINCOS_CACHE)
1076 1118 if (!m_has_lb_cache) {
1077 1119 m_has_lb_cache = true;
1078   - m_sin_b = std::sin(m_b);
1079   - m_cos_b = std::cos(m_b);
  1120 + m_sin_b = std::sin(b());
  1121 + m_cos_b = std::cos(b());
1080 1122 }
1081 1123 #endif
1082   - if (dir.m_has_lb) {
1083   - arg_1 = std::sin(dir.m_l - m_l);
1084   - #if defined(G_SINCOS_CACHE)
1085   - arg_2 = m_cos_b * std::tan(dir.m_b) -
1086   - m_sin_b * std::cos(dir.m_l - m_l);
1087   - #else
1088   - arg_2 = std::cos(m_b) * std::tan(dir.m_b) -
1089   - std::sin(m_b) * std::cos(dir.m_l - m_l);
1090   - #endif
1091   - }
1092   - else {
1093   - arg_1 = std::sin(dir.l() - m_l);
1094   - #if defined(G_SINCOS_CACHE)
1095   - arg_2 = m_cos_b * std::tan(dir.b()) -
1096   - m_sin_b * std::cos(dir.l() - m_l);
1097   - #else
1098   - arg_2 = std::cos(m_b) * std::tan(dir.b()) -
1099   - std::sin(m_b) * std::cos(dir.l() - m_l);
1100   - #endif
1101   - }
1102   - }
1103   - else if (m_has_radec) {
  1124 +
  1125 + // Compute arguments of arctan
  1126 + arg_1 = std::sin(dir.l() - l());
1104 1127 #if defined(G_SINCOS_CACHE)
1105   - if (!m_has_radec_cache) {
1106   - m_has_radec_cache = true;
1107   - m_sin_dec = std::sin(m_dec);
1108   - m_cos_dec = std::cos(m_dec);
1109   - }
  1128 + arg_2 = m_cos_b * std::tan(dir.b()) -
  1129 + m_sin_b * std::cos(dir.l() - l());
  1130 + #else
  1131 + arg_2 = std::cos(b()) * std::tan(dir.b()) -
  1132 + std::sin(b()) * std::cos(dir.l() - l());
1110 1133 #endif
1111   - if (dir.m_has_radec) {
1112   - arg_1 = std::sin(dir.m_ra - m_ra);
1113   - #if defined(G_SINCOS_CACHE)
1114   - arg_2 = m_cos_dec * std::tan(dir.m_dec) -
1115   - m_sin_dec * std::cos(dir.m_ra - m_ra);
1116   - #else
1117   - arg_2 = std::cos(m_dec) * std::tan(dir.m_dec) -
1118   - std::sin(m_dec) * std::cos(dir.m_ra - m_ra);
1119   - #endif
1120   - }
1121   - else {
1122   - arg_1 = std::sin(dir.ra() - m_ra);
1123   - #if defined(G_SINCOS_CACHE)
1124   - arg_2 = m_cos_dec * std::tan(dir.dec()) -
1125   - m_sin_dec * std::cos(dir.ra() - m_ra);
1126   - #else
1127   - arg_2 = std::cos(m_dec) * std::tan(dir.dec()) -
1128   - std::sin(m_dec) * std::cos(dir.ra() - m_ra);
1129   - #endif
1130   - }
1131   - }
  1134 +
  1135 + } // endelse: computation in Galactic coordinates
  1136 +
  1137 + // ... otherwise throw an exception
1132 1138 else {
1133   - arg_1 = std::sin(dir.ra() - ra());
1134   - arg_2 = std::cos(dec())*std::tan(dir.dec()) -
1135   - std::sin(dec())*std::cos(dir.ra() - ra());
  1139 + std::string msg = "Invalid coordinate system \""+coordsys+"\" "
  1140 + "specified. Please specify either \"CEL\" or "
  1141 + "\"GAL\".";
  1142 + throw GException::invalid_argument(G_POSANG, msg);
1136 1143 }
1137 1144  
1138 1145 // Compute position angle
... ...
test/test_GSky.cpp
  View file @964b56a
... ... @@ -315,7 +315,7 @@ void TestGSky::test_GSkyDir(void)
315 315 test_value(dir2.l_deg(), 184.55745, "Longitude of Crab");
316 316 test_value(dir2.b_deg(), -5.78436, "Latitude of Crab");
317 317  
318   - // Test precess method
  318 + // Test precess() method
319 319 dir2.precess(2000.0, 1950.0);
320 320 test_value(dir2.ra_deg(), 82.88086799, "Right Ascension of Crab (epoch 1950)");
321 321 test_value(dir2.dec_deg(), 21.98181049, "Declination of Crab (epoch 1950)");
... ... @@ -323,6 +323,38 @@ void TestGSky::test_GSkyDir(void)
323 323 test_value(dir2.ra_deg(), 83.633083, "Right Ascension of Crab");
324 324 test_value(dir2.dec_deg(), 22.0145, "Declination of Crab");
325 325  
  326 + // Test posang() and posang_deg() methods
  327 + dir1.radec_deg(0.0, 0.0);
  328 + dir2.radec_deg(0.0, 10.0);
  329 + test_value(dir1.posang(dir2), 0.0, "1. Position angle posang()");
  330 + test_value(dir1.posang(dir2, "CEL"), 0.0, "2. Position angle posang(CEL)");
  331 + test_value(dir1.posang(dir2, "GAL"), 0.4098000, "3. Position angle posang(GAL)");
  332 + test_value(dir1.posang_deg(dir2), 0.0, "4. Position angle posang_deg()");
  333 + test_value(dir1.posang_deg(dir2, "CEL"), 0.0, "5. Position angle posang_deg(CEL)");
  334 + test_value(dir1.posang_deg(dir2, "GAL"), 23.47981, "6. Position angle posang_deg(GAL)");
  335 + dir2.radec_deg(90.0, 0.0);
  336 + test_value(dir1.posang(dir2), 1.5707963, "7. Position angle posang()");
  337 + test_value(dir1.posang(dir2, "CEL"), 1.5707963, "8. Position angle posang(CEL)");
  338 + test_value(dir1.posang(dir2, "GAL"), 1.9805963, "9. Position angle posang(GAL)");
  339 + test_value(dir1.posang_deg(dir2), 90.0, "10. Position angle posang_deg()");
  340 + test_value(dir1.posang_deg(dir2, "CEL"), 90.0, "11. Position angle posang_deg(CEL)");
  341 + test_value(dir1.posang_deg(dir2, "GAL"), 113.47981, "12. Position angle posang_deg(GAL)");
  342 + dir1.lb_deg(0.0, 0.0);
  343 + dir2.lb_deg(0.0, 10.0);
  344 + test_value(dir1.posang(dir2), -1.0227397, "13. Position angle posang()");
  345 + test_value(dir1.posang(dir2, "CEL"), -1.0227397, "14. Position angle posang(CEL)");
  346 + test_value(dir1.posang(dir2, "GAL"), 0.0, "15. Position angle posang(GAL)");
  347 + test_value(dir1.posang_deg(dir2), -58.5986663, "16. Position angle posang_deg()");
  348 + test_value(dir1.posang_deg(dir2, "CEL"), -58.5986663, "17. Position angle posang_deg(CEL)");
  349 + test_value(dir1.posang_deg(dir2, "GAL"), 0.0, "18. Position angle posang_deg(GAL)");
  350 + dir2.lb_deg(90.0, 0.0);
  351 + test_value(dir1.posang(dir2), 0.5480567, "19. Position angle posang()");
  352 + test_value(dir1.posang(dir2, "CEL"), 0.5480567, "20. Position angle posang(CEL)");
  353 + test_value(dir1.posang(dir2, "GAL"), 1.5707963, "21. Position angle posang(GAL)");
  354 + test_value(dir1.posang_deg(dir2), 31.4013337, "22. Position angle posang_deg()");
  355 + test_value(dir1.posang_deg(dir2, "CEL"), 31.4013337, "23. Position angle posang_deg(CEL)");
  356 + test_value(dir1.posang_deg(dir2, "GAL"), 90.0, "24. Position angle posang_deg(GAL)");
  357 +
326 358 // Return
327 359 return;
328 360 }
... ...