GWcs::wcs_set() method now sets radesys and equinox members (#2866)

Jürgen Knödlseder
1 parent 1ea91214
Showing 5 changed files with 102 additions and 6 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
include/GWcs.hpp
src/sky/GWcs.cpp
-2019-04-03
+2019-04-08
  
         * Version 1.6.0 released
         ========================
  
+        GWcs::wcs_set() method now sets radesys and equinox members (#2866)
         Add GCTAPointing::is_valid() method and remove table methods
         Set CTA IRF integration iterations to 6 to cope with H.E.S.S. PSF (#2860)
         Add kludge to handle incomplete Fermi-LAT source maps (#2850)
 New Features and Important Changes in GammaLib 1.6.0
  
-3 April 2019
+8 April 2019
  
  
 1. Introduction
@@ -251,6 +251,10 @@ None
  
 12. Sky module
 --------------
+The GWcs::wcs_set() method now sets m_radesys and m_equinox members through the
+new protected GWcs::wcs_set_radesys() method. By default, sky maps now have
+a RADECSYS=ICRS keyword (#2866).
+
 An OpenMP locking mechanism was added to GWcs::dir2pix that improves the speed of
 the code parallelization by avoiding mutual lockups of the threads (#2770).
  
@@ -19,6 +19,8 @@ In particular, this release provides:
 Bug fixes
 ---------
  
+* [`2866 <https://cta-redmine.irap.omp.eu/issues/2866>`_] -
+  ``GWcs::wcs_set()`` method now sets radesys and equinox members
 * [`2860 <https://cta-redmine.irap.omp.eu/issues/2860>`_] -
   Set CTA IRF integration iterations to 6 to cope with H.E.S.S. PSF
 * [`2838 <https://cta-redmine.irap.omp.eu/issues/2838>`_] -
@@ -62,7 +62,7 @@ public:
     virtual ~GWcs(void);
  
     // Operators
-    virtual GWcs& operator= (const GWcs& wcs);
+    virtual GWcs& operator=(const GWcs& wcs);
  
     // Pure virtual base class methods
     virtual void        clear(void) = 0;
@@ -104,11 +104,12 @@ protected:
                              const double& crpix1, const double& crpix2,
                              const double& cdelt1, const double& cdelt2);
     virtual bool compare(const GSkyProjection& proj) const;
-    bool         undefined(const double& value) const { return (value==UNDEFINED); }
+    bool         undefined(const double& value) const;
  
     // Methods adapted from wcslib::wcs.c 
     void        wcs_ini(int naxis);
     void        wcs_set(void) const;
+    void        wcs_set_radesys(void) const;
     void        wcs_set_ctype(void) const;
     void        wcs_p2s(int ncoord, int nelem, const double* pixcrd, double* imgcrd,
                         double* phi, double* theta, double* world, int* stat) const;
@@ -177,8 +178,8 @@ protected:
     mutable double                   m_latpole; //!< LATPOLEa keyvalue
     double                           m_restfrq; //!< RESTFRQa keyvalue
     double                           m_restwav; //!< RESTWAVa keyvalue
-    std::string                      m_radesys; //!< RADESYS keyvalue
-    double                           m_equinox; //!< EQUINOX keyvalue
+    mutable std::string              m_radesys; //!< RADESYS keyvalue
+    mutable double                   m_equinox; //!< EQUINOX keyvalue
     std::vector<double>              m_cd;      //!< CDi_ja linear transformation matrix
     std::vector<double>              m_crota;   //!< CROTAia keyvalues for each coord axis
     int                              m_lng;     //!< Longitude axis
@@ -252,6 +253,19 @@ int GWcs::size(void) const
  
  
 /***********************************************************************//**
+ * @brief Check if value is undefined
+ * 
+ * @param[in] value Value to check
+ * @return True if @p value is undefined, false otherwise
+ ***************************************************************************/
+inline
+bool GWcs::undefined(const double& value) const
+{
+    return (value==UNDEFINED);
+}
+
+
+/***********************************************************************//**
  * @brief Initializes an OpenMP lock with a specific name
  *
  * @param[in] lock_id Identifier of the lock that should be initialized
@@ -299,6 +299,16 @@ void GWcs::read(const GFitsHDU&amp; hdu)
         // Set standard parameters
         set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
  
+        // Optionally read "RADESYS" and "EQUINOX" keywords (needs to come
+        // before wcs_set() since the method will eventually complement
+        // missing keywords)
+        if (hdu.has_card("RADESYS")) {
+            m_radesys = hdu.string("RADESYS");
+        }
+        if (hdu.has_card("EQUINOX")) {
+            m_equinox = hdu.real("EQUINOX");
+        }
+
         // Setup WCS derived parameters
         wcs_set();
  
@@ -1282,6 +1292,9 @@ void GWcs::wcs_set(void) const
     // Initialize the linear transformation
     lin_set();
  
+    // Set defaults for radesys and equinox
+    wcs_set_radesys();
+
     // Signal that WCS is set
     m_wcsset = true;
  
@@ -1291,6 +1304,68 @@ void GWcs::wcs_set(void) const
  
  
 /***********************************************************************//**
+ * @brief Set radesys and equinox members
+ *
+ * Set default values for radesys and equinox for equatorial, ecliptic or
+ * helioecliptic coordinates. For other coordinates the radesys and equinox
+ * members are set to empty/undefined values.
+ *
+ * This method has been inspired by code from wcs.c::wcsset.
+ ***************************************************************************/
+void GWcs::wcs_set_radesys(void) const
+{
+    // Set defaults for radesys and equinox for equatorial, ecliptic or
+    // helioecliptic coordinates
+    if ((m_coordsys == 0) || (m_coordsys == 2) || (m_coordsys == 3)) {
+
+        // If radesys is undefined then set it dependend on the equinox.
+        // If equinox is undefined then use "ICRS", otherwise use "FK4" if
+        // equinox is smaller than 1984 and "FK5" for later equinoxes
+        if (m_radesys.empty()) {
+            if (undefined(m_equinox)) {
+                m_radesys = "ICRS";
+            }
+            else if (m_equinox < 1984.0) {
+                m_radesys = "FK4";
+            }
+            else {
+                m_radesys = "FK5";
+            }
+        }
+
+        // ... otherwise, if radesys is "ICRS" or "GAPPT" then set the equinox
+        // to undefined since the equinox is not applicable for these coordinate
+        // systems
+        else if ((m_radesys == "ICRS") || (m_radesys == "GAPPT")) {
+            m_equinox = UNDEFINED;
+        }
+
+        // ... otherwise, if equinox is undefined and radesys is "FK5" then
+        // set equinox to 2000.0 and if RADESYS is "FK4" or "FK4-NO-E" then
+        // set equinox to 1950.0
+        else if (undefined(m_equinox)) {
+            if (m_radesys == "FK5") {
+                m_equinox = 2000.0;
+            }
+            else if ((m_radesys == "FK4") || (m_radesys == "FK4-NO-E")) {
+                m_equinox = 1950.0;
+            }
+        }
+
+    } // endif: coordinate system requires RADESYS
+
+    // ... otherwise reset radesys and set equinox to undefined
+    else {
+        m_radesys.clear();
+        m_equinox = UNDEFINED;
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Set CTYPEa keywords
  *
  * @exception GException::wcs_invalid