From 964b56a964c0a5e28bf3d3f8cbd1cfcdf5b4c77a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 13 Apr 2020 17:52:32 +0200 Subject: [PATCH] Fix ambiguity problem in GSkyDir::posang() and posang_deg() (#3199) --- ChangeLog | 3 ++- NEWS | 9 ++++++++- doc/source/admin/release_history/1.7.rst | 4 +++- include/GSkyDir.hpp | 24 ++++++++++++++++-------- inst/com/include/GCOMOad.hpp | 9 ++++----- inst/cta/src/GCTAOnOffObservation.cpp | 4 ++-- inst/cta/src/GCTAPointing.cpp | 4 ++-- inst/cta/src/GCTAResponseCube.cpp | 4 ++-- inst/cta/src/GCTAResponseIrf.cpp | 10 +++++----- inst/spi/include/GSPIResponse.hpp | 4 ++-- inst/spi/src/GSPIResponse.cpp | 6 +----- pyext/GSkyDir.i | 10 ++++++---- src/model/GModelSpatialElliptical.cpp | 4 ++-- src/sky/GSkyDir.cpp | 147 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------- test/test_GSky.cpp | 34 +++++++++++++++++++++++++++++++++- 15 files changed, 165 insertions(+), 111 deletions(-) diff --git a/ChangeLog b/ChangeLog index e252f80..2e58a53 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2020-02-17 +2020-04-13 * Version 1.7.0 released ======================== + Fix ambiguity problem in GSkyDir::posang() and posang_deg() (#3199) Avoid rounding errors in GSkyMap operators (#3126) Add GObservations::npred() methods (#3062) Exclude upper edge boundary in GSkyMap::contains() method (#3012) diff --git a/NEWS b/NEWS index 57a3917..3703a71 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.7.0 -11 April 2020 +13 April 2020 1. Introduction @@ -56,6 +56,8 @@ The following methods have been renamed: The arguments for the following methods have been changed: - GTime::utc(int&) - Add optional precision value +- GSkyDir::posang(GSkyDir&, std::string&) - Add optional coordinate system string +- GSkyDir::posang_dir(GSkyDir&, std::string&) - Add optional coordinate system string The return value of the following methods has been changed: - none @@ -142,6 +144,11 @@ None 12. Sky module -------------- +The position angle is now computed by GSkyDir::posang() and GSkyDir::posang_deg() +for the coordinate system that is specified by the coordsys argument. By default +this argument is set to "CEL", hence by default the methods return the position +angle with respect to celestial coordinates (#3199). + The GSkyMap operators used the sky map value access operator that bi-linearly interpolates sky map values even for maps with identical definition. This could lead to rounding errors. For maps with identical definition the operators now diff --git a/doc/source/admin/release_history/1.7.rst b/doc/source/admin/release_history/1.7.rst index 9e6b7fd..491b72c 100644 --- a/doc/source/admin/release_history/1.7.rst +++ b/doc/source/admin/release_history/1.7.rst @@ -12,12 +12,14 @@ In particular, this release provides: * Support for INTEGRAL/SPI data analysis * Response caching for all instruments -* Support for extended models for all instruments +* ... Bug fixes --------- +* [`3199 `_] - + Fix ambiguity problem in ``GSkyDir::posang()`` and ``GSkyDir::posang_deg()`` * [`3126 `_] - Avoid rounding errors in ``GSkyMap`` operators. * [`3012 `_] - diff --git a/include/GSkyDir.hpp b/include/GSkyDir.hpp index c35481a..a539e21 100644 --- a/include/GSkyDir.hpp +++ b/include/GSkyDir.hpp @@ -87,8 +87,8 @@ public: void rotate(const double& phi, const double& theta); void rotate_deg(const double& phi, const double& theta); void precess(const double& from_epoch, const double& to_epoch); - void sun(const GTime& time, const double& epoch = 2000.0); - void moon(const GTime& time, const double& epoch = 2000.0); + void sun(const GTime& time, const double& epoch = 2000.0); + void moon(const GTime& time, const double& epoch = 2000.0); const double& l(void) const; const double& b(void) const; const double& ra(void) const; @@ -101,8 +101,10 @@ public: double cos_dist(const GSkyDir& dir) const; double dist(const GSkyDir& dir) const; double dist_deg(const GSkyDir& dir) const; - double posang(const GSkyDir& dir) const; - double posang_deg(const GSkyDir& dir) const; + double posang(const GSkyDir& dir, + const std::string& coordsys = "CEL") const; + double posang_deg(const GSkyDir& dir, + const std::string& coordsys = "CEL") const; std::string print(const GChatter& chatter = NORMAL) const; private: @@ -289,17 +291,23 @@ double GSkyDir::dist_deg(const GSkyDir& dir) const * @brief Compute position angle between sky directions in degrees * * @param[in] dir Sky direction. - * @return Position angle in degrees. + * @param[in] coordsys Coordinate system ("CEL" or "GAL") + * @return Position angle (deg). * - * Returns the position angle of a sky direction in degrees. + * Computes the position angle of a specified sky direction with respect to + * the sky direction of the instance in degrees. If "CEL" is specified as + * @p coordsys the position angle is counted counterclockwise from celestial + * North, if "GAL" is specified the position angle is counted + * counterclockwise from Galactic North. * * See the posang() method for more information about the computation of the * position angle. ***************************************************************************/ inline -double GSkyDir::posang_deg(const GSkyDir& dir) const +double GSkyDir::posang_deg(const GSkyDir& dir, + const std::string& coordsys) const { - return (posang(dir) * gammalib::rad2deg); + return (posang(dir, coordsys) * gammalib::rad2deg); } #endif /* GSKYDIR_HPP */ diff --git a/inst/com/include/GCOMOad.hpp b/inst/com/include/GCOMOad.hpp index 0fc47c9..91b0597 100644 --- a/inst/com/include/GCOMOad.hpp +++ b/inst/com/include/GCOMOad.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMOad.hpp - COMPTEL Orbit Aspect Data class * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Juergen Knodlseder * + * copyright (C) 2017-2020 by Juergen Knodlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -409,12 +409,11 @@ double GCOMOad::theta(const GSkyDir& sky) const inline double GCOMOad::phi(const GSkyDir& sky) const { - // If position angle has not be initialised the do it now + // If position angle has not be initialised then do it now if (m_posang > 1.0e20) { - m_posang = m_zaxis.posang_deg(m_xaxis); + m_posang = m_zaxis.posang_deg(m_xaxis); // Celestial system } - //return (m_zaxis.posang_deg(m_xaxis) - m_zaxis.posang_deg(sky)); - return (m_posang - m_zaxis.posang_deg(sky)); + return (m_posang - m_zaxis.posang_deg(sky)); // Celestial system } #endif /* GCOMOAD_HPP */ diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp index c123898..3a1a637 100644 --- a/inst/cta/src/GCTAOnOffObservation.cpp +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -1582,7 +1582,7 @@ void GCTAOnOffObservation::compute_arf_cut(const GCTAObservation& obs, // Compute position of pixel centre in instrument coordinates double theta = obsdir.dist(pixdir); - double phi = obsdir.posang(pixdir); + double phi = obsdir.posang(pixdir); // Celestial system // Add up effective area m_arf[i] += rsp.aeff(theta, phi, @@ -2068,7 +2068,7 @@ void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs, // Compute position of pixel centre in instrument coordinates double theta = obsdir.dist(pixdir); - double phi = obsdir.posang(pixdir); + double phi = obsdir.posang(pixdir); // Celestial system // Loop over true energy for (int itrue = 0; itrue < ntrue; ++itrue) { diff --git a/inst/cta/src/GCTAPointing.cpp b/inst/cta/src/GCTAPointing.cpp index 8c04541..569166c 100644 --- a/inst/cta/src/GCTAPointing.cpp +++ b/inst/cta/src/GCTAPointing.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAPointing.cpp - CTA pointing class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2019 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -262,7 +262,7 @@ GCTAInstDir GCTAPointing::instdir(const GSkyDir& skydir) const } #else double theta = m_dir.dist(skydir); - double phi = m_dir.posang(skydir); + double phi = m_dir.posang(skydir); // Celestial system double detx(0.0); double dety(0.0); if (theta > 0.0 ) { diff --git a/inst/cta/src/GCTAResponseCube.cpp b/inst/cta/src/GCTAResponseCube.cpp index b59c871..34b83fb 100644 --- a/inst/cta/src/GCTAResponseCube.cpp +++ b/inst/cta/src/GCTAResponseCube.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAResponseCube.cpp - CTA cube analysis response function class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2018 by Juergen Knoedlseder * + * copyright (C) 2014-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -1407,7 +1407,7 @@ double GCTAResponseCube::irf_elliptical(const GEvent& event, // Compute angle between model centre and measured photon direction and // position angle (radians) double rho_obs = model->dir().dist(obsDir); - double posangle_obs = model->dir().posang(obsDir); + double posangle_obs = model->dir().posang(obsDir); // Celestial // Get livetime (in seconds) double livetime = exposure().livetime(); diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index 0db513b..f533e27 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAResponseIrf.cpp - CTA instrument response function class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2019 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -2519,12 +2519,12 @@ double GCTAResponseIrf::irf_elliptical(const GEvent& event, // centre and position angle of observed photon direction seen from the // model centre [radians] double rho_obs = centre.dist(obsDir); - double posangle_obs = centre.posang(obsDir); + double posangle_obs = centre.posang(obsDir); // Celestial system // Determine angular distance between model centre and pointing direction // [radians] double rho_pnt = centre.dist(pnt.dir()); - double posangle_pnt = centre.posang(pnt.dir()); + double posangle_pnt = centre.posang(pnt.dir()); // Celestial system // Compute azimuth angle of pointing in model coordinate system [radians] double omega_pnt = posangle_pnt - posangle_obs; @@ -2957,7 +2957,7 @@ double GCTAResponseIrf::nroi_radial(const GModelSky& model, // Compute position angle of ROI centre with respect to model // centre (radians) - double omega0 = centre.posang(roi.centre().dir()); + double omega0 = centre.posang(roi.centre().dir()); // Celestial system // Setup integration kernel cta_nroi_radial_kern_rho integrand(this, @@ -3162,7 +3162,7 @@ double GCTAResponseIrf::nroi_elliptical(const GModelSky& model, // Compute position angle of ROI centre with respect to model // centre (radians) - double posangle_roi = centre.posang(roi.centre().dir()); + double posangle_roi = centre.posang(roi.centre().dir()); // Celestial system // Setup integration kernel cta_nroi_elliptical_kern_rho integrand(this, diff --git a/inst/spi/include/GSPIResponse.hpp b/inst/spi/include/GSPIResponse.hpp index 4c97d87..782e338 100644 --- a/inst/spi/include/GSPIResponse.hpp +++ b/inst/spi/include/GSPIResponse.hpp @@ -130,7 +130,7 @@ private: // Private cache std::vector m_spix; //!< SPI pointing direction - std::vector m_posang; //!< Position angle of Y axis (radians) + std::vector m_posang; //!< Position angle of Y axis (CEL, radians) mutable bool m_has_wcs; //!< Has WCS information mutable double m_wcs_xmin; //!< Minimum X value (radians) mutable double m_wcs_ymin; //!< Minimum Y value (radians) @@ -241,7 +241,7 @@ double GSPIResponse::zenith(const int& ipt, const GSkyDir& dir) const inline double GSPIResponse::azimuth(const int& ipt, const GSkyDir& dir) const { - double azimuth = m_posang[ipt] - m_spix[ipt].posang(dir); + double azimuth = m_posang[ipt] - m_spix[ipt].posang(dir); // Celestial system if (azimuth < 0.0) { azimuth += gammalib::twopi; } diff --git a/inst/spi/src/GSPIResponse.cpp b/inst/spi/src/GSPIResponse.cpp index 94deb13..e3146db 100644 --- a/inst/spi/src/GSPIResponse.cpp +++ b/inst/spi/src/GSPIResponse.cpp @@ -727,10 +727,6 @@ void GSPIResponse::init_members(void) m_wcs_ypix_max = 0.0; m_max_zenith = 180.0 * gammalib::deg2rad; - // Set energy scale for response cache to MeV - //m_irf_cache.energy_scale(GEnergy(1.0, "MeV")); - //m_nroi_cache.energy_scale(GEnergy(1.0, "MeV")); - // Return return; } @@ -1287,7 +1283,7 @@ void GSPIResponse::set_cache(const GSPIEventCube* cube) // Loop over all pointings for (int ipt = 0; ipt < npt; ++ipt) { - // Compute position angle + // Compute position angle is celestial coordinates double posang = cube->spi_x(ipt).posang(cube->spi_z(ipt)) + gammalib::pihalf; // Store angles diff --git a/pyext/GSkyDir.i b/pyext/GSkyDir.i index 878837a..42cf195 100644 --- a/pyext/GSkyDir.i +++ b/pyext/GSkyDir.i @@ -55,8 +55,8 @@ public: void rotate(const double& phi, const double& theta); void rotate_deg(const double& phi, const double& theta); void precess(const double& from_epoch, const double& to_epoch); - void sun(const GTime& time, const double& epoch = gammalib::epoch2000); - void moon(const GTime& time, const double& epoch = gammalib::epoch2000); + void sun(const GTime& time, const double& epoch = 2000.0); + void moon(const GTime& time, const double& epoch = 2000.0); const double& l(void) const; const double& b(void) const; const double& ra(void) const; @@ -69,8 +69,10 @@ public: double cos_dist(const GSkyDir& dir) const; double dist(const GSkyDir& dir) const; double dist_deg(const GSkyDir& dir) const; - double posang(const GSkyDir& dir) const; - double posang_deg(const GSkyDir& dir) const; + double posang(const GSkyDir& dir, + const std::string& coordsys = "CEL") const; + double posang_deg(const GSkyDir& dir, + const std::string& coordsys = "CEL") const; }; diff --git a/src/model/GModelSpatialElliptical.cpp b/src/model/GModelSpatialElliptical.cpp index 6d5f1df..259b420 100644 --- a/src/model/GModelSpatialElliptical.cpp +++ b/src/model/GModelSpatialElliptical.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialElliptical.cpp - Abstract elliptical spatial model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2018 by Juergen Knoedlseder * + * copyright (C) 2013-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -176,7 +176,7 @@ double GModelSpatialElliptical::eval(const GPhoton& photon, // Compute distance from source and position angle (in radians) const GSkyDir& srcDir = photon.dir(); double theta = dir().dist(srcDir); - double posang = dir().posang(srcDir); + double posang = dir().posang(srcDir); // Celestial system // Evaluate model double value = eval(theta, posang, photon.energy(), photon.time(), gradients); diff --git a/src/sky/GSkyDir.cpp b/src/sky/GSkyDir.cpp index 38ebaad..527a4bb 100644 --- a/src/sky/GSkyDir.cpp +++ b/src/sky/GSkyDir.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyDir.cpp - Sky direction class * * ----------------------------------------------------------------------- * - * copyright (C) 2008-2019 by Juergen Knoedlseder * + * copyright (C) 2008-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -38,6 +38,7 @@ #include "GVector.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_POSANG "GSkyDir::posang(GSkyDir&, std::string&)" /* __ Macros _____________________________________________________________ */ @@ -507,7 +508,7 @@ void GSkyDir::sun(const GTime& time, const double& epoch) while (mean_longitude >= 360.0) { mean_longitude -= 360.0; } - + // Compute the mean anomaly of the Sun in radians double mean_anomaly = (357.528 + 0.9856003 * n) * gammalib::deg2rad; @@ -872,7 +873,7 @@ void GSkyDir::moon(const GTime& time, const double& epoch) // Add nutation in longitude geolong += nut_long / 3600.0; - + // nutate, jd, nlong, elong // geolong= geolong + nlong/3.6d3 // cirrange,geolong @@ -1049,90 +1050,96 @@ double GSkyDir::cos_dist(const GSkyDir& dir) const * @brief Compute position angle between sky directions in radians * * @param[in] dir Sky direction. - * @return Position angle in radians. + * @param[in] coordsys Coordinate system ("CEL" or "GAL") + * @return Position angle (rad). + * + * Computes the position angle of a specified sky direction with respect to + * the sky direction of the instance in radians. If "CEL" is specified as + * @p coordsys the position angle is counted counterclockwise from celestial + * North, if "GAL" is specified the position angle is counted + * counterclockwise from Galactic North. + * + * If @p coordsys is "CEL" the position angle is computed using * - * Computes the position angle using - * \f[PA = \arctan \left( - * \frac{\sin( \alpha_1 - \alpha_0 )} - * {\cos \delta_0 \tan \delta_1 - - * \sin \delta_0 \cos(\alpha_1 - \alpha_0)} \right)\f] + * \f[PA = \arctan \left( + * \frac{\sin( \alpha - \alpha_0 )} + * {\cos \delta_0 \tan \delta - + * \sin \delta_0 \cos(\alpha - \alpha_0)} \right)\f] * where - * \f$(\alpha_0,\delta_0)\f$ are the coordinates of the reference point, and - * \f$(\alpha_1,\delta_1)\f$ are the coordinates of sky direction for which - * the position angle is to be computed. + * * \f$(\alpha_0,\delta_0)\f$ are the celestial coordinates of the instance, + * and + * * \f$(\alpha,\delta)\f$ are the celestial coordinates of sky direction. * - * The position angle is counted counterclockwise from north. + * If @p coordsys is "GAL" the position angle is computed using + * + * \f[PA = \arctan \left( + * \frac{\sin( \l - \l_0 )} + * {\cos \b_0 \tan \b - \sin \b_0 \cos(\l - \l_0)} \right)\f] + * where + * * \f$(\l_0,\l_0)\f$ are the Galactic coordinates of the instance, + * and + * * \f$(\l,\b)\f$ are the Galactic coordinates of sky direction. ***************************************************************************/ -double GSkyDir::posang(const GSkyDir& dir) const +double GSkyDir::posang(const GSkyDir& dir, const std::string& coordsys) const { // Initialise arguments of arctan double arg_1; double arg_2; - // Compute dependent on coordinate system availability. This speeds - // up things by avoiding unnecessary coordinate transformations. - if (m_has_lb) { + // Case A: Computation in celestial coordinates + if (coordsys == "CEL") { + + // If required then compute sin(dec) and cos(dec) cache values + #if defined(G_SINCOS_CACHE) + if (!m_has_radec_cache) { + m_has_radec_cache = true; + m_sin_dec = std::sin(dec()); + m_cos_dec = std::cos(dec()); + } + #endif + + // Compute arguments of arctan + arg_1 = std::sin(dir.ra() - ra()); + #if defined(G_SINCOS_CACHE) + arg_2 = m_cos_dec * std::tan(dir.dec()) - + m_sin_dec * std::cos(dir.ra() - ra()); + #else + arg_2 = std::cos(dec()) * std::tan(dir.dec()) - + std::sin(dec()) * std::cos(dir.ra() - ra()); + #endif + + } // endif: computation in celestial coordinates + + // Case B: Computation in Galactic coordinates + else if (coordsys == "GAL") { + + // If required then compute sin(b) and cos(b) cache values #if defined(G_SINCOS_CACHE) if (!m_has_lb_cache) { m_has_lb_cache = true; - m_sin_b = std::sin(m_b); - m_cos_b = std::cos(m_b); + m_sin_b = std::sin(b()); + m_cos_b = std::cos(b()); } #endif - if (dir.m_has_lb) { - arg_1 = std::sin(dir.m_l - m_l); - #if defined(G_SINCOS_CACHE) - arg_2 = m_cos_b * std::tan(dir.m_b) - - m_sin_b * std::cos(dir.m_l - m_l); - #else - arg_2 = std::cos(m_b) * std::tan(dir.m_b) - - std::sin(m_b) * std::cos(dir.m_l - m_l); - #endif - } - else { - arg_1 = std::sin(dir.l() - m_l); - #if defined(G_SINCOS_CACHE) - arg_2 = m_cos_b * std::tan(dir.b()) - - m_sin_b * std::cos(dir.l() - m_l); - #else - arg_2 = std::cos(m_b) * std::tan(dir.b()) - - std::sin(m_b) * std::cos(dir.l() - m_l); - #endif - } - } - else if (m_has_radec) { + + // Compute arguments of arctan + arg_1 = std::sin(dir.l() - l()); #if defined(G_SINCOS_CACHE) - if (!m_has_radec_cache) { - m_has_radec_cache = true; - m_sin_dec = std::sin(m_dec); - m_cos_dec = std::cos(m_dec); - } + arg_2 = m_cos_b * std::tan(dir.b()) - + m_sin_b * std::cos(dir.l() - l()); + #else + arg_2 = std::cos(b()) * std::tan(dir.b()) - + std::sin(b()) * std::cos(dir.l() - l()); #endif - if (dir.m_has_radec) { - arg_1 = std::sin(dir.m_ra - m_ra); - #if defined(G_SINCOS_CACHE) - arg_2 = m_cos_dec * std::tan(dir.m_dec) - - m_sin_dec * std::cos(dir.m_ra - m_ra); - #else - arg_2 = std::cos(m_dec) * std::tan(dir.m_dec) - - std::sin(m_dec) * std::cos(dir.m_ra - m_ra); - #endif - } - else { - arg_1 = std::sin(dir.ra() - m_ra); - #if defined(G_SINCOS_CACHE) - arg_2 = m_cos_dec * std::tan(dir.dec()) - - m_sin_dec * std::cos(dir.ra() - m_ra); - #else - arg_2 = std::cos(m_dec) * std::tan(dir.dec()) - - std::sin(m_dec) * std::cos(dir.ra() - m_ra); - #endif - } - } + + } // endelse: computation in Galactic coordinates + + // ... otherwise throw an exception else { - arg_1 = std::sin(dir.ra() - ra()); - arg_2 = std::cos(dec())*std::tan(dir.dec()) - - std::sin(dec())*std::cos(dir.ra() - ra()); + std::string msg = "Invalid coordinate system \""+coordsys+"\" " + "specified. Please specify either \"CEL\" or " + "\"GAL\"."; + throw GException::invalid_argument(G_POSANG, msg); } // Compute position angle diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index dbdfc77..ef30abd 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -315,7 +315,7 @@ void TestGSky::test_GSkyDir(void) test_value(dir2.l_deg(), 184.55745, "Longitude of Crab"); test_value(dir2.b_deg(), -5.78436, "Latitude of Crab"); - // Test precess method + // Test precess() method dir2.precess(2000.0, 1950.0); test_value(dir2.ra_deg(), 82.88086799, "Right Ascension of Crab (epoch 1950)"); test_value(dir2.dec_deg(), 21.98181049, "Declination of Crab (epoch 1950)"); @@ -323,6 +323,38 @@ void TestGSky::test_GSkyDir(void) test_value(dir2.ra_deg(), 83.633083, "Right Ascension of Crab"); test_value(dir2.dec_deg(), 22.0145, "Declination of Crab"); + // Test posang() and posang_deg() methods + dir1.radec_deg(0.0, 0.0); + dir2.radec_deg(0.0, 10.0); + test_value(dir1.posang(dir2), 0.0, "1. Position angle posang()"); + test_value(dir1.posang(dir2, "CEL"), 0.0, "2. Position angle posang(CEL)"); + test_value(dir1.posang(dir2, "GAL"), 0.4098000, "3. Position angle posang(GAL)"); + test_value(dir1.posang_deg(dir2), 0.0, "4. Position angle posang_deg()"); + test_value(dir1.posang_deg(dir2, "CEL"), 0.0, "5. Position angle posang_deg(CEL)"); + test_value(dir1.posang_deg(dir2, "GAL"), 23.47981, "6. Position angle posang_deg(GAL)"); + dir2.radec_deg(90.0, 0.0); + test_value(dir1.posang(dir2), 1.5707963, "7. Position angle posang()"); + test_value(dir1.posang(dir2, "CEL"), 1.5707963, "8. Position angle posang(CEL)"); + test_value(dir1.posang(dir2, "GAL"), 1.9805963, "9. Position angle posang(GAL)"); + test_value(dir1.posang_deg(dir2), 90.0, "10. Position angle posang_deg()"); + test_value(dir1.posang_deg(dir2, "CEL"), 90.0, "11. Position angle posang_deg(CEL)"); + test_value(dir1.posang_deg(dir2, "GAL"), 113.47981, "12. Position angle posang_deg(GAL)"); + dir1.lb_deg(0.0, 0.0); + dir2.lb_deg(0.0, 10.0); + test_value(dir1.posang(dir2), -1.0227397, "13. Position angle posang()"); + test_value(dir1.posang(dir2, "CEL"), -1.0227397, "14. Position angle posang(CEL)"); + test_value(dir1.posang(dir2, "GAL"), 0.0, "15. Position angle posang(GAL)"); + test_value(dir1.posang_deg(dir2), -58.5986663, "16. Position angle posang_deg()"); + test_value(dir1.posang_deg(dir2, "CEL"), -58.5986663, "17. Position angle posang_deg(CEL)"); + test_value(dir1.posang_deg(dir2, "GAL"), 0.0, "18. Position angle posang_deg(GAL)"); + dir2.lb_deg(90.0, 0.0); + test_value(dir1.posang(dir2), 0.5480567, "19. Position angle posang()"); + test_value(dir1.posang(dir2, "CEL"), 0.5480567, "20. Position angle posang(CEL)"); + test_value(dir1.posang(dir2, "GAL"), 1.5707963, "21. Position angle posang(GAL)"); + test_value(dir1.posang_deg(dir2), 31.4013337, "22. Position angle posang_deg()"); + test_value(dir1.posang_deg(dir2, "CEL"), 31.4013337, "23. Position angle posang_deg(CEL)"); + test_value(dir1.posang_deg(dir2, "GAL"), 90.0, "24. Position angle posang_deg(GAL)"); + // Return return; } -- libgit2 0.21.4