Eirik Hatlen / gammalib

Commit e47ea36954eda1ba469477b4662a1864b8f620f9

Authored by Jürgen Knödlseder
1 parent e9d77e48

Add ARC WCS projection (#2541)

Inline Side-by-side
Showing 18 changed files with 909 additions and 16 deletions
ChangeLog
  View file @e47ea36
1   -2018-07-03
  1 +2018-07-30
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Add ARC WCS projection (#2541)
6 7 Add GLog::buffer() method.
7 8 Renamed GPhoton::mcid() to GPhoton::mc_id()
8 9 Remove at() methods from Python interface (#2515)
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -5 July 2018
  3 +30 July 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -13,6 +13,7 @@ made since the last release of GammaLib.
13 13 ---------------------------------
14 14 The following classes have been added:
15 15 - GModelSpectralExponential
  16 +- GWcsARC
16 17  
17 18 The following classes have been removed:
18 19 - none
... ... @@ -205,6 +206,10 @@ None
205 206  
206 207 12. Sky module
207 208 --------------
  209 +The class GWcsARC was added to handle ARC WCS projections (#2541). In addition
  210 +the private GWcs::prj_bchk() method was added to implement boundary checking
  211 +similar to the checking done in WCS.
  212 +
208 213 A precomputation cache was added to the GWcs::dir2pix() and GWcs::pix2dir()
209 214 methods so that results are not recomputed for the same arguments. This should
210 215 provide some speed-up as it avoids unnecessary coordinate transformations.
... ...
doc/source/admin/release_history/1.6.rst
  View file @e47ea36
... ... @@ -34,6 +34,8 @@ Bug fixes
34 34 Improvements
35 35 ------------
36 36  
  37 +* [`2541 <https://cta-redmine.irap.omp.eu/issues/2541>`_] -
  38 + Add ``ARC`` WCS projection.
37 39 * Add ``GLog::buffer()`` method
38 40 * Renamed ``GPhoton::mcid()`` to ``GPhoton::mc_id()``
39 41 * [`2513 <https://cta-redmine.irap.omp.eu/issues/2513>`_] -
... ...
include/GWcs.hpp
  View file @e47ea36
... ... @@ -146,6 +146,13 @@ protected:
146 146 // Methods adapted from wcslib::prj.c
147 147 void prj_ini(void) const;
148 148 void prj_off(const double& phi0, const double& theta0) const;
  149 + int prj_bchk(const double& tol,
  150 + const int& nphi,
  151 + const int& ntheta,
  152 + const int& spt,
  153 + double* phi,
  154 + double* theta,
  155 + int* stat) const;
149 156 virtual void prj_set(void) const = 0;
150 157 virtual void prj_x2s(int nx, int ny, int sxy, int spt,
151 158 const double* x, const double* y,
... ...
include/GWcsARC.hpp 0 → 100644
  View file @e47ea36
  1 +/***************************************************************************
  2 + * GWcsARC.hpp - Zenithal/azimuthal equidistant (ARC) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GWcsARC.hpp
  23 + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GWCSARC_HPP
  28 +#define GWCSARC_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GWcs.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @class GWcsARC
  37 + *
  38 + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition
  39 + *
  40 + * This class implements the Zenithal/azimuthal equidistant projection for
  41 + * the World Coordinate System.
  42 + ***************************************************************************/
  43 +class GWcsARC : public GWcs {
  44 +
  45 +public:
  46 + // Constructors and destructors
  47 + GWcsARC(void);
  48 + GWcsARC(const std::string& coords,
  49 + const double& crval1, const double& crval2,
  50 + const double& crpix1, const double& crpix2,
  51 + const double& cdelt1, const double& cdelt2);
  52 + GWcsARC(const GWcsARC& wcs);
  53 + virtual ~GWcsARC(void);
  54 +
  55 + // Operators
  56 + GWcsARC& operator=(const GWcsARC& wcs);
  57 +
  58 + // Implemented pure virtual base class methods
  59 + virtual void clear(void);
  60 + virtual GWcsARC* clone(void) const;
  61 + virtual std::string classname(void) const;
  62 + virtual std::string code(void) const;
  63 + virtual std::string name(void) const;
  64 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  65 +
  66 +private:
  67 + // Private methods
  68 + void init_members(void);
  69 + void copy_members(const GWcsARC& wcs);
  70 + void free_members(void);
  71 + void prj_set(void) const;
  72 + void prj_x2s(int nx, int ny, int sxy, int spt,
  73 + const double* x, const double* y,
  74 + double* phi, double* theta, int* stat) const;
  75 + void prj_s2x(int nphi, int ntheta, int spt, int sxy,
  76 + const double* phi, const double* theta,
  77 + double* x, double* y, int* stat) const;
  78 +};
  79 +
  80 +
  81 +/***********************************************************************//**
  82 + * @brief Return class name
  83 + *
  84 + * @return String containing the class name ("GWcsARC").
  85 + ***************************************************************************/
  86 +inline
  87 +std::string GWcsARC::classname(void) const
  88 +{
  89 + return ("GWcsARC");
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return projection code
  95 + *
  96 + * @return Projection code.
  97 + *
  98 + * Returns the projection code "ARC".
  99 + ***************************************************************************/
  100 +inline
  101 +std::string GWcsARC::code(void) const
  102 +{
  103 + return "ARC";
  104 +}
  105 +
  106 +
  107 +/***********************************************************************//**
  108 + * @brief Return projection name
  109 + *
  110 + * @return Projection name.
  111 + *
  112 + * Returns the projection name.
  113 + ***************************************************************************/
  114 +inline
  115 +std::string GWcsARC::name(void) const
  116 +{
  117 + return "Zenithal/azimuthal equidistant";
  118 +}
  119 +
  120 +#endif /* GWCSARC_HPP */
... ...
include/GammaLib.hpp
  View file @e47ea36
1 1 /***************************************************************************
2 2 * GammaLib.hpp - Gamma-Ray Astronomy Library Header file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2006-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -155,6 +155,7 @@
155 155 #include "GWcsRegistry.hpp"
156 156 #include "GWcs.hpp"
157 157 #include "GWcsAIT.hpp"
  158 +#include "GWcsARC.hpp"
158 159 #include "GWcsAZP.hpp"
159 160 #include "GWcsCAR.hpp"
160 161 #include "GWcsGLS.hpp"
... ...
include/Makefile.am
  View file @e47ea36
1 1 #############################################################################
2 2 # Makefile.am for include module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2017 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2018 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -117,6 +117,7 @@ pkginclude_HEADERS = GammaLib.hpp \
117 117 GWcsRegistry.hpp \
118 118 GWcs.hpp \
119 119 GWcsAIT.hpp \
  120 + GWcsARC.hpp \
120 121 GWcsAZP.hpp \
121 122 GWcsCAR.hpp \
122 123 GWcsGLS.hpp \
... ...
pyext/GWcsARC.i 0 → 100644
  View file @e47ea36
  1 +/***************************************************************************
  2 + * GWcsARC.i - Zenithal/azimuthal equidistant (ARC) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GWcsARC.i
  23 + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GWcsARC.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GWcsARC
  34 + *
  35 + * @brief Zenithal/azimuthal equidistant (ARC) projection class
  36 + ***************************************************************************/
  37 +class GWcsARC : public GWcs {
  38 +public:
  39 + // Constructors and destructors
  40 + GWcsARC(void);
  41 + GWcsARC(const std::string& coords,
  42 + const double& crval1, const double& crval2,
  43 + const double& crpix1, const double& crpix2,
  44 + const double& cdelt1, const double& cdelt2);
  45 + GWcsARC(const GWcsARC& wcs);
  46 + virtual ~GWcsARC(void);
  47 +
  48 + // Implemented pure virtual methods
  49 + virtual void clear(void);
  50 + virtual GWcsARC* clone(void) const;
  51 + virtual std::string classname(void) const;
  52 + virtual std::string code(void) const;
  53 + virtual std::string name(void) const;
  54 +};
  55 +
  56 +
  57 +/***********************************************************************//**
  58 + * @brief GWcsARC class extension
  59 + ***************************************************************************/
  60 +%extend GWcsARC {
  61 + GWcsARC copy() {
  62 + return (*self);
  63 + }
  64 +};
... ...
pyext/gammalib/sky.i
  View file @e47ea36
1 1 /***************************************************************************
2 2 * sky.i - Sky module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2017 by Jurgen Knodlseder *
  4 + * copyright (C) 2011-2018 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -122,6 +122,7 @@ namespace std {
122 122 %include "GWcsRegistry.i"
123 123 %include "GWcs.i"
124 124 %include "GWcsAIT.i"
  125 +%include "GWcsARC.i"
125 126 %include "GWcsAZP.i"
126 127 %include "GWcsCAR.i"
127 128 %include "GWcsGLS.i"
... ...
src/sky/GWcs.cpp
  View file @e47ea36
... ... @@ -3090,3 +3090,109 @@ void GWcs::prj_off(const double&amp; phi0, const double&amp; theta0) const
3090 3090 // Return
3091 3091 return;
3092 3092 }
  3093 +
  3094 +
  3095 +/***********************************************************************//**
  3096 + * @brief Performs bounds checking on native spherical coordinates
  3097 + *
  3098 + * @param[in] tol Tolerance for the bounds check [deg]
  3099 + * @param[in] nphi Number of bins in phi.
  3100 + * @param[in] ntheta Number of bins in theta.
  3101 + * @param[in] spt Step size.
  3102 + * @param[in,out] phi Pointer to phi array [deg]
  3103 + * @param[in,out] theta Pointer to theta array [deg]
  3104 + * @param[in,out] stat Pointer to status array.
  3105 + *
  3106 + * Performs bounds checking on native spherical coordinates. As returned by
  3107 + * the deprojection (x2s) routines, native longitude is expected to lie in
  3108 + * the closed interval [-180,180], with latitude in [-90,90].
  3109 + *
  3110 + * A tolerance may be specified to provide a small allowance for numerical
  3111 + * imprecision. Values that lie outside the allowed range by not more than
  3112 + * the specified tolerance will be adjusted back into range.
  3113 + *
  3114 + * Code adapted from prj.c::prjbchk().
  3115 + ***************************************************************************/
  3116 +int GWcs::prj_bchk(const double& tol,
  3117 + const int& nphi,
  3118 + const int& ntheta,
  3119 + const int& spt,
  3120 + double* phi,
  3121 + double* theta,
  3122 + int* stat) const
  3123 +{
  3124 + // Initialize result
  3125 + int status = 0;
  3126 +
  3127 + // Continue only if bounds checking is requested
  3128 + if (m_bounds) {
  3129 +
  3130 + // Initialise array pointers
  3131 + register double* phip = phi;
  3132 + register double* thetap = theta;
  3133 + register int* statp = stat;
  3134 +
  3135 + // Loop over all theta bins
  3136 + for (int itheta = 0; itheta < ntheta; ++itheta) {
  3137 +
  3138 + // Loop over all phi bins
  3139 + for (int iphi = 0; iphi < nphi; ++iphi, phip += spt, thetap += spt, statp++) {
  3140 +
  3141 + // Skip values already marked as illegal
  3142 + if (*statp == 0) {
  3143 +
  3144 + // Check lower phi boundary
  3145 + if (*phip < -180.0) {
  3146 + if (*phip < -180.0-tol) {
  3147 + *statp = 1;
  3148 + status = 1;
  3149 + }
  3150 + else {
  3151 + *phip = -180.0;
  3152 + }
  3153 + }
  3154 +
  3155 + // Check upper phi boundary
  3156 + else if (180.0 < *phip) {
  3157 + if (180.0+tol < *phip) {
  3158 + *statp = 1;
  3159 + status = 1;
  3160 + }
  3161 + else {
  3162 + *phip = 180.0;
  3163 + }
  3164 + }
  3165 +
  3166 + // Check lower theta boundary
  3167 + if (*thetap < -90.0) {
  3168 + if (*thetap < -90.0-tol) {
  3169 + *statp = 1;
  3170 + status = 1;
  3171 + }
  3172 + else {
  3173 + *thetap = -90.0;
  3174 + }
  3175 + }
  3176 +
  3177 + // Check upper theta boundary
  3178 + else if (90.0 < *thetap) {
  3179 + if (90.0+tol < *thetap) {
  3180 + *statp = 1;
  3181 + status = 1;
  3182 + }
  3183 + else {
  3184 + *thetap = 90.0;
  3185 + }
  3186 + }
  3187 +
  3188 + } // endif: values were marked as valid
  3189 +
  3190 + } // endfor: looped over phi
  3191 +
  3192 + } // endfor: looped over theta
  3193 +
  3194 + } // endif: bounds checking was requested
  3195 +
  3196 + // Return status
  3197 + return status;
  3198 +}
... ...
src/sky/GWcsAIT.cpp
  View file @e47ea36
1 1 /***************************************************************************
2   - * GWcsAIT.cpp - Aitoff (AIT) projection class *
  2 + * GWcsAIT.cpp - Zenithal/azimuthal equidistant (AIT) projection class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2013-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -306,21 +306,21 @@ void GWcsAIT::prj_set(void) const
306 306  
307 307 // Initialise projection parameters
308 308 m_w.assign(4, 0.0);
309   -
  309 +
310 310 // Set undefined parameters
311 311 if (m_r0 == 0.0) {
312 312 m_r0 = gammalib::rad2deg;
313 313 }
314   -
  314 +
315 315 // Precompute
316 316 m_w[0] = 2.0 * m_r0 * m_r0;
317 317 m_w[1] = 1.0 / (2.0 * m_w[0]);
318 318 m_w[2] = m_w[1] / 4.0;
319 319 m_w[3] = 1.0 / (2.0 * m_r0);
320   -
  320 +
321 321 // Compute fiducial offset
322 322 prj_off(0.0, 0.0);
323   -
  323 +
324 324 // Return
325 325 return;
326 326 }
... ... @@ -381,7 +381,7 @@ void GWcsAIT::prj_x2s(int nx, int ny, int sxy, int spt,
381 381 // Initialise status code and statistics
382 382 int status = 0;
383 383 int n_invalid = 0;
384   -
  384 +
385 385 // Do x dependence
386 386 const double* xp = x;
387 387 int rowoff = 0;
... ... @@ -480,6 +480,11 @@ void GWcsAIT::prj_x2s(int nx, int ny, int sxy, int spt,
480 480 throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
481 481 }
482 482  
  483 + // Do bounds checking on the native coordinates
  484 + // if (prj->bounds&4 && prjbchk(1.0e-13, nx, my, spt, phi, theta, stat)) {
  485 + // if (!status) status = PRJERR_BAD_PIX_SET("aitx2s");
  486 + // }
  487 +
483 488 // Return
484 489 return;
485 490 }
... ... @@ -537,7 +542,7 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy,
537 542 // Initialise status code and statistics
538 543 int status = 0;
539 544 int n_invalid = 0;
540   -
  545 +
541 546 // Do phi dependence
542 547 const double* phip = phi;
543 548 int rowoff = 0;
... ... @@ -563,7 +568,7 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy,
563 568 double* yp = y;
564 569 int* statp = stat;
565 570 for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
566   -
  571 +
567 572 // Compute sin(theta) and cos(theta)
568 573 double sinthe;
569 574 double costhe;
... ... @@ -578,12 +583,12 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy,
578 583 } // endfor: phi
579 584  
580 585 } // endfor: theta
581   -
  586 +
582 587 // Handle status code
583 588 if (status == 4) {
584 589 throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
585 590 }
586   -
  591 +
587 592 // Return
588 593 return;
589 594 }
... ...
src/sky/GWcsARC.cpp 0 → 100644
  View file @e47ea36
  1 +/***************************************************************************
  2 + * GWcsARC.cpp - Zenithal/azimuthal equidistant (ARC) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GWcsARC.cpp
  23 + * @brief Zenithal/azimuthal equidistant (ARC) projection class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GMath.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GWcsARC.hpp"
  35 +#include "GWcsRegistry.hpp"
  36 +
  37 +/* __ Method name definitions ____________________________________________ */
  38 +#define G_PRJ_X2S "GWcsARC::prj_x2s(int, int, int, int, double*, double*,"\
  39 + " double*, double*, int*)"
  40 +#define G_PRJ_S2X "GWcsARC::prj_s2x(int, int, int, int, double*, double*,"\
  41 + " double*, double*, int*)"
  42 +
  43 +/* __ Macros _____________________________________________________________ */
  44 +
  45 +/* __ Coding definitions _________________________________________________ */
  46 +
  47 +/* __ Debug definitions __________________________________________________ */
  48 +
  49 +/* __ Local prototypes ___________________________________________________ */
  50 +
  51 +/* __ Constants __________________________________________________________ */
  52 +
  53 +/* __ Globals ____________________________________________________________ */
  54 +const GWcsARC g_wcs_arc_seed;
  55 +const GWcsRegistry g_wcs_arc_registry(&g_wcs_arc_seed);
  56 +
  57 +
  58 +/*==========================================================================
  59 + = =
  60 + = Constructors/destructors =
  61 + = =
  62 + ==========================================================================*/
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Void constructor
  66 + ***************************************************************************/
  67 +GWcsARC::GWcsARC(void) : GWcs()
  68 +{
  69 + // Initialise class members
  70 + init_members();
  71 +
  72 + // Return
  73 + return;
  74 +}
  75 +
  76 +
  77 +/***********************************************************************//**
  78 + * @brief Constructor
  79 + *
  80 + * @param[in] coords Coordinate system.
  81 + * @param[in] crval1 X value of reference pixel.
  82 + * @param[in] crval2 Y value of reference pixel.
  83 + * @param[in] crpix1 X index of reference pixel (starting from 1).
  84 + * @param[in] crpix2 Y index of reference pixel (starting from 1).
  85 + * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
  86 + * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
  87 + ***************************************************************************/
  88 +GWcsARC::GWcsARC(const std::string& coords,
  89 + const double& crval1, const double& crval2,
  90 + const double& crpix1, const double& crpix2,
  91 + const double& cdelt1, const double& cdelt2) :
  92 + GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
  93 +
  94 +{
  95 + // Initialise class members
  96 + init_members();
  97 +
  98 + // Setup WCS derived parameters
  99 + wcs_set();
  100 +
  101 + // Return
  102 + return;
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Copy constructor
  108 + *
  109 + * @param[in] wcs World Coordinate System.
  110 + ***************************************************************************/
  111 +GWcsARC::GWcsARC(const GWcsARC& wcs) : GWcs(wcs)
  112 +{
  113 + // Initialise class members for clean destruction
  114 + init_members();
  115 +
  116 + // Copy members
  117 + copy_members(wcs);
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Destructor
  126 + ***************************************************************************/
  127 +GWcsARC::~GWcsARC(void)
  128 +{
  129 + // Free members
  130 + free_members();
  131 +
  132 + // Return
  133 + return;
  134 +}
  135 +
  136 +
  137 +/*==========================================================================
  138 + = =
  139 + = Operators =
  140 + = =
  141 + ==========================================================================*/
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Assignment operator
  145 + *
  146 + * @param[in] wcs World Coordinate System.
  147 + * @return World Coordinate System.
  148 + ***************************************************************************/
  149 +GWcsARC& GWcsARC::operator=(const GWcsARC& wcs)
  150 +{
  151 + // Execute only if object is not identical
  152 + if (this != &wcs) {
  153 +
  154 + // Copy base class members
  155 + this->GWcs::operator=(wcs);
  156 +
  157 + // Free members
  158 + free_members();
  159 +
  160 + // Initialise private members for clean destruction
  161 + init_members();
  162 +
  163 + // Copy members
  164 + copy_members(wcs);
  165 +
  166 + } // endif: object was not identical
  167 +
  168 + // Return this object
  169 + return *this;
  170 +}
  171 +
  172 +
  173 +/*==========================================================================
  174 + = =
  175 + = Public methods =
  176 + = =
  177 + ==========================================================================*/
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Clear instance
  181 + *
  182 + * This method properly resets the object to an initial state.
  183 + ***************************************************************************/
  184 +void GWcsARC::clear(void)
  185 +{
  186 + // Free class members (base and derived classes, derived class first)
  187 + free_members();
  188 + this->GWcs::free_members();
  189 + this->GSkyProjection::free_members();
  190 +
  191 + // Initialise members
  192 + this->GSkyProjection::init_members();
  193 + this->GWcs::init_members();
  194 + init_members();
  195 +
  196 + // Return
  197 + return;
  198 +}
  199 +
  200 +
  201 +/***********************************************************************//**
  202 + * @brief Clone instance
  203 + *
  204 + * @return Pointer to deep copy of Zenithal/azimuthal equidistant projection.
  205 + ***************************************************************************/
  206 +GWcsARC* GWcsARC::clone(void) const
  207 +{
  208 + return new GWcsARC(*this);
  209 +}
  210 +
  211 +
  212 +
  213 +/***********************************************************************//**
  214 + * @brief Print WCS information
  215 + *
  216 + * @param[in] chatter Chattiness.
  217 + * @return String containing WCS information.
  218 + ***************************************************************************/
  219 +std::string GWcsARC::print(const GChatter& chatter) const
  220 +{
  221 + // Initialise result string
  222 + std::string result;
  223 +
  224 + // Continue only if chatter is not silent
  225 + if (chatter != SILENT) {
  226 +
  227 + // Append header
  228 + result.append("=== GWcsARC ===");
  229 +
  230 + // Append information
  231 + result.append(wcs_print(chatter));
  232 +
  233 + } // endif: chatter was not silent
  234 +
  235 + // Return result
  236 + return result;
  237 +}
  238 +
  239 +
  240 +/*==========================================================================
  241 + = =
  242 + = Private methods =
  243 + = =
  244 + ==========================================================================*/
  245 +
  246 +/***********************************************************************//**
  247 + * @brief Initialise class members
  248 + ***************************************************************************/
  249 +void GWcsARC::init_members(void)
  250 +{
  251 + // Return
  252 + return;
  253 +}
  254 +
  255 +
  256 +/***********************************************************************//**
  257 + * @brief Copy class members
  258 + *
  259 + * @param[in] wcs World Coordinate System.
  260 + ***************************************************************************/
  261 +void GWcsARC::copy_members(const GWcsARC& wcs)
  262 +{
  263 + // Return
  264 + return;
  265 +}
  266 +
  267 +
  268 +/***********************************************************************//**
  269 + * @brief Delete class members
  270 + ***************************************************************************/
  271 +void GWcsARC::free_members(void)
  272 +{
  273 + // Return
  274 + return;
  275 +}
  276 +
  277 +
  278 +/***********************************************************************//**
  279 + * @brief Setup of projection
  280 + *
  281 + * This method sets up the projection information. The method has been
  282 + * adapted from the wcslib function prj.c::arcset.
  283 + *
  284 + * Given and/or returned:
  285 + * m_r0 Reset to 180/pi if 0.
  286 + * m_phi0 Reset to 0.0 if undefined.
  287 + * m_theta0 Reset to 90.0 if undefined.
  288 + *
  289 + * Returned:
  290 + * m_x0 Fiducial offset in x.
  291 + * m_y0 Fiducial offset in y.
  292 + * m_w[0] r0*(pi/180)
  293 + * m_w[1] (180/pi)/r0
  294 + ***************************************************************************/
  295 +void GWcsARC::prj_set(void) const
  296 +{
  297 + // Signal that projection has been set (needs to be done before calling
  298 + // the prj_off() method to avoid an endless loop)
  299 + m_prjset = true;
  300 +
  301 + // Initialise projection parameters
  302 + m_w.assign(2, 0.0);
  303 +
  304 + // Set undefined parameters
  305 + if (m_r0 == 0.0) {
  306 + m_r0 = gammalib::rad2deg;
  307 + m_w[0] = 1.0;
  308 + m_w[1] = 1.0;
  309 + }
  310 + else {
  311 + m_w[0] = m_r0 * gammalib::deg2rad;
  312 + m_w[1] = 1.0/m_w[0];
  313 + }
  314 +
  315 + // Compute fiducial offset
  316 + prj_off(0.0, 0.0);
  317 +
  318 + // Return
  319 + return;
  320 +}
  321 +
  322 +
  323 +/***********************************************************************//**
  324 + * @brief Pixel-to-spherical deprojection
  325 + *
  326 + * @param[in] nx X vector length.
  327 + * @param[in] ny Y vector length (0=no replication).
  328 + * @param[in] sxy Input vector step.
  329 + * @param[in] spt Output vector step.
  330 + * @param[in] x Vector of projected x coordinates.
  331 + * @param[in] y Vector of projected y coordinates.
  332 + * @param[out] phi Longitude of the projected point in native spherical
  333 + * coordinates [deg].
  334 + * @param[out] theta Latitude of the projected point in native spherical
  335 + * coordinates [deg].
  336 + * @param[out] stat Status return value for each vector element (always 0)
  337 + *
  338 + * @exception GException::wcs_invalid_x_y
  339 + * One or more of the (x,y) coordinates were invalid, as indicated
  340 + * by the stat vector.
  341 + *
  342 + * Deproject pixel (x,y) coordinates in the plane of projection to native
  343 + * spherical coordinates (phi,theta).
  344 + *
  345 + * This method has been adapted from the wcslib function prj.c::aitx2s().
  346 + * The interface follows very closely that of wcslib. In contrast to the
  347 + * wcslib routine, however, the method assumes that the projection has been
  348 + * setup previously (as this will be done by the constructor).
  349 + ***************************************************************************/
  350 +void GWcsARC::prj_x2s(int nx, int ny, int sxy, int spt,
  351 + const double* x, const double* y,
  352 + double* phi, double* theta, int* stat) const
  353 +{
  354 + // Initialize projection if required
  355 + if (!m_prjset) {
  356 + prj_set();
  357 + }
  358 +
  359 + // Set value replication length mx,my
  360 + int mx;
  361 + int my;
  362 + if (ny > 0) {
  363 + mx = nx;
  364 + my = ny;
  365 + }
  366 + else {
  367 + mx = 1;
  368 + my = 1;
  369 + ny = nx;
  370 + }
  371 +
  372 + // Do x dependence
  373 + const double* xp = x;
  374 + int rowoff = 0;
  375 + int rowlen = nx * spt;
  376 + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
  377 + double xj = *xp + m_x0;
  378 + double* phip = phi + rowoff;
  379 + for (int iy = 0; iy < my; ++iy, phip += rowlen) {
  380 + *phip = xj;
  381 + }
  382 + }
  383 +
  384 + // Do y dependence
  385 + const double* yp = y;
  386 + double* phip = phi;
  387 + double* thetap = theta;
  388 + int* statp = stat;
  389 + for (int iy = 0; iy < ny; ++iy, yp += sxy) {
  390 + double yj = *yp + m_y0;
  391 + double yj2 = yj * yj;
  392 + for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
  393 + double xj = *phip;
  394 + double r = std::sqrt(xj*xj + yj2);
  395 + if (r == 0.0) {
  396 + *phip = 0.0;
  397 + *thetap = 90.0;
  398 + }
  399 + else {
  400 + *phip = gammalib::atan2d(xj, -yj);
  401 + *thetap = 90.0 - r * m_w[1];
  402 + }
  403 + *(statp++) = 0;
  404 + }
  405 + }
  406 +
  407 + // Do boundary checking
  408 + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat);
  409 +
  410 + // Handle status code
  411 + if (status != 0) {
  412 + std::string msg = "One or more of the (x, y) coordinates were invalid "
  413 + "for the "+name()+" ("+code()+") projection.";
  414 + throw GException::invalid_value(G_PRJ_X2S, msg);
  415 + }
  416 +
  417 + // Return
  418 + return;
  419 +}
  420 +
  421 +
  422 +/***********************************************************************//**
  423 + * @brief Generic spherical-to-pixel projection
  424 + *
  425 + * @param[in] nphi Longitude vector length.
  426 + * @param[in] ntheta Latitude vector length (0=no replication).
  427 + * @param[in] spt Input vector step.
  428 + * @param[in] sxy Output vector step.
  429 + * @param[in] phi Longitude vector of the projected point in native spherical
  430 + * coordinates [deg].
  431 + * @param[in] theta Latitude vector of the projected point in native spherical
  432 + * coordinates [deg].
  433 + * @param[out] x Vector of projected x coordinates.
  434 + * @param[out] y Vector of projected y coordinates.
  435 + * @param[out] stat Status return value for each vector element (always 0)
  436 + *
  437 + * @exception GException::wcs_invalid_phi_theta
  438 + * One or more of the (phi,theta) coordinates were invalid, as
  439 + * indicated by the stat vector.
  440 + *
  441 + * Project native spherical coordinates (phi,theta) to pixel (x,y)
  442 + * coordinates in the plane of projection.
  443 + *
  444 + * This method has been adapted from the wcslib function prj.c::aits2x().
  445 + * The interface follows very closely that of wcslib. In contrast to the
  446 + * wcslib routine, however, the method assumes that the projection has been
  447 + * setup previously (as this will be done by the constructor).
  448 + ***************************************************************************/
  449 +void GWcsARC::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  450 + const double* phi, const double* theta,
  451 + double* x, double* y, int* stat) const
  452 +{
  453 + // Initialize projection if required
  454 + if (!m_prjset) {
  455 + prj_set();
  456 + }
  457 +
  458 + // Set value replication length mphi,mtheta
  459 + int mphi;
  460 + int mtheta;
  461 + if (ntheta > 0) {
  462 + mphi = nphi;
  463 + mtheta = ntheta;
  464 + }
  465 + else {
  466 + mphi = 1;
  467 + mtheta = 1;
  468 + ntheta = nphi;
  469 + }
  470 +
  471 + // Initialise status code and statistics
  472 + //int status = 0;
  473 + //int n_invalid = 0;
  474 +
  475 + // Do phi dependence
  476 + const double* phip = phi;
  477 + int rowoff = 0;
  478 + int rowlen = nphi * sxy;
  479 + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
  480 + double sinphi;
  481 + double cosphi;
  482 + gammalib::sincosd(*phip, &sinphi, &cosphi);
  483 + double* xp = x + rowoff;
  484 + double* yp = y + rowoff;
  485 + for (int itheta = 0; itheta < mtheta; ++itheta) {
  486 + *xp = sinphi;
  487 + *yp = cosphi;
  488 + xp += rowlen;
  489 + yp += rowlen;
  490 + }
  491 + }
  492 +
  493 + // Do theta dependence
  494 + const double* thetap = theta;
  495 + double* xp = x;
  496 + double* yp = y;
  497 + int* statp = stat;
  498 + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
  499 +
  500 + // Compute r
  501 + double r = m_w[0] * (90.0 - *thetap);
  502 +
  503 + // Do phi dependence
  504 + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
  505 + *xp = r * (*xp) - m_x0;
  506 + *yp = -r * (*yp) - m_y0;
  507 + *(statp++) = 0;
  508 + } // endfor: phi
  509 +
  510 + } // endfor: theta
  511 +
  512 + // Handle status code
  513 + /*
  514 + if (status == 4) {
  515 + throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
  516 + }
  517 + */
  518 +
  519 + // Return
  520 + return;
  521 +}
... ...
src/sky/GWcsCAR.cpp
  View file @e47ea36
... ... @@ -34,6 +34,8 @@
34 34 #include "GWcsRegistry.hpp"
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
  37 +#define G_PRJ_X2S "GWcsCAR::prj_x2s(int, int, int, int, double*, double*,"\
  38 + " double*, double*, int*)"
37 39  
38 40 /* __ Macros _____________________________________________________________ */
39 41  
... ... @@ -383,6 +385,16 @@ void GWcsCAR::prj_x2s(int nx, int ny, int sxy, int spt,
383 385 }
384 386 }
385 387  
  388 + // Do boundary checking
  389 + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat);
  390 +
  391 + // Handle status code
  392 + if (status != 0) {
  393 + std::string msg = "One or more of the (x, y) coordinates were invalid "
  394 + "for the "+name()+" ("+code()+") projection.";
  395 + throw GException::invalid_value(G_PRJ_X2S, msg);
  396 + }
  397 +
386 398 // Return
387 399 return;
388 400 }
... ...
src/sky/GWcsMER.cpp
  View file @e47ea36
... ... @@ -34,6 +34,8 @@
34 34 #include "GWcsRegistry.hpp"
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
  37 +#define G_PRJ_X2S "GWcsMER::prj_x2s(int, int, int, int, double*, double*,"\
  38 + " double*, double*, int*)"
37 39  
38 40 /* __ Macros _____________________________________________________________ */
39 41  
... ... @@ -383,6 +385,16 @@ void GWcsMER::prj_x2s(int nx, int ny, int sxy, int spt,
383 385 }
384 386 }
385 387  
  388 + // Do boundary checking
  389 + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat);
  390 +
  391 + // Handle status code
  392 + if (status != 0) {
  393 + std::string msg = "One or more of the (x, y) coordinates were invalid "
  394 + "for the "+name()+" ("+code()+") projection.";
  395 + throw GException::invalid_value(G_PRJ_X2S, msg);
  396 + }
  397 +
386 398 // Return
387 399 return;
388 400 }
... ...
src/sky/GWcsMOL.cpp
  View file @e47ea36
... ... @@ -529,6 +529,16 @@ void GWcsMOL::prj_x2s(int nx, int ny, int sxy, int spt,
529 529 throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
530 530 }
531 531  
  532 + // Do boundary checking
  533 + status = prj_bchk(1.0e-11, nx, my, spt, phi, theta, stat);
  534 +
  535 + // Handle status code
  536 + if (status != 0) {
  537 + std::string msg = "One or more of the (x, y) coordinates were invalid "
  538 + "for the "+name()+" ("+code()+") projection.";
  539 + throw GException::invalid_value(G_PRJ_X2S, msg);
  540 + }
  541 +
532 542 // Return
533 543 return;
534 544 }
... ...
src/sky/GWcsSFL.cpp
  View file @e47ea36
... ... @@ -428,6 +428,16 @@ void GWcsSFL::prj_x2s(int nx, int ny, int sxy, int spt,
428 428 throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
429 429 }
430 430  
  431 + // Do boundary checking
  432 + status = prj_bchk(1.0e-12, nx, my, spt, phi, theta, stat);
  433 +
  434 + // Handle status code
  435 + if (status != 0) {
  436 + std::string msg = "One or more of the (x, y) coordinates were invalid "
  437 + "for the "+name()+" ("+code()+") projection.";
  438 + throw GException::invalid_value(G_PRJ_X2S, msg);
  439 + }
  440 +
431 441 // Return
432 442 return;
433 443 }
... ...
src/sky/Makefile.am
  View file @e47ea36
1 1 #############################################################################
2 2 # Makefile.am for sky module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2017 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2018 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -38,6 +38,7 @@ sources = GSkyDir.cpp \
38 38 GWcsRegistry.cpp \
39 39 GWcs.cpp \
40 40 GWcsAIT.cpp \
  41 + GWcsARC.cpp \
41 42 GWcsAZP.cpp \
42 43 GWcsCAR.cpp \
43 44 GWcsGLS.cpp \
... ...
test/test_GSky.py
  View file @e47ea36
... ... @@ -358,6 +358,17 @@ class Test(gammalib.GPythonTestSuite):
358 358 # Return
359 359 return
360 360  
  361 + # Test ARC projection
  362 + def _test_skymap_arc(self):
  363 + """
  364 + Test ARC projection.
  365 + """
  366 + # Execute generic test
  367 + self._test_skymap_proj('ARC')
  368 +
  369 + # Return
  370 + return
  371 +
361 372 # Test AZP projection
362 373 def _test_skymap_azp(self):
363 374 """
... ... @@ -565,6 +576,7 @@ class Test(gammalib.GPythonTestSuite):
565 576 test_support._pickeling(self, gammalib.GSkyRegionMap())
566 577 test_support._pickeling(self, gammalib.GSkyRegions())
567 578 test_support._pickeling(self, gammalib.GWcsAIT())
  579 + test_support._pickeling(self, gammalib.GWcsARC())
568 580 test_support._pickeling(self, gammalib.GWcsAZP())
569 581 test_support._pickeling(self, gammalib.GWcsCAR())
570 582 test_support._pickeling(self, gammalib.GWcsGLS())
... ... @@ -596,6 +608,7 @@ class Test(gammalib.GPythonTestSuite):
596 608 test_support._pickeling(self, gammalib.GSkyRegionMap(map))
597 609 test_support._pickeling(self, gammalib.GSkyRegions(regions))
598 610 test_support._pickeling(self, gammalib.GWcsAIT('GAL',1.0,2.0,1.0,1.0,0.1,0.1))
  611 + test_support._pickeling(self, gammalib.GWcsARC('GAL',1.0,2.0,1.0,1.0,0.1,0.1))
599 612 test_support._pickeling(self, gammalib.GWcsAZP('GAL',1.0,2.0,1.0,1.0,0.1,0.1))
600 613 test_support._pickeling(self, gammalib.GWcsCAR('GAL',1.0,2.0,1.0,1.0,0.1,0.1))
601 614 test_support._pickeling(self, gammalib.GWcsGLS('GAL',1.0,2.0,1.0,1.0,0.1,0.1))
... ... @@ -622,6 +635,7 @@ class Test(gammalib.GPythonTestSuite):
622 635 self.append(self._test_operators, 'Test sky map operators')
623 636 self.append(self._test_skymap_healpix, 'Test HEALPix map')
624 637 self.append(self._test_skymap_ait, 'Test AIT projection map')
  638 + self.append(self._test_skymap_arc, 'Test ARC projection map')
625 639 self.append(self._test_skymap_azp, 'Test AZP projection map')
626 640 self.append(self._test_skymap_car, 'Test CAR projection map')
627 641 self.append(self._test_skymap_gls, 'Test GLS projection map')
... ...