Luigi Tibaldo / gammalib

Commit f048480f8dcf9547d45e9d5af3ea779e7d50a448

Authored by Jürgen Knödlseder
1 parent 800ade35

Add Mollweide's WCS projection

Inline Side-by-side
Showing 11 changed files with 880 additions and 7 deletions
ChangeLog
  View file @f048480
1   -2015-02-03 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2015-02-04 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.0.0 released
4 4 ========================
... ... @@ -11,6 +11,7 @@
11 11 Rename CTA cube analysis classes (#1269)
12 12 Do not dump unqueried parameters (#1349)
13 13 Add GBilinear class
  14 + Add WCS projections (MOL)
14 15  
15 16  
16 17 2015-01-07 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.0.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -3 February 2015
  5 +4 February 2015
6 6  
7 7  
8 8 1. Introduction
... ... @@ -78,7 +78,8 @@ None
78 78  
79 79 12. Sky module
80 80 --------------
81   -None
  81 +The following WCS projections have been added:
  82 +* GWcsMOL (Mollweide's projection)
82 83  
83 84  
84 85 13. Support module
... ...
include/GWcsMOL.hpp 0 → 100644
  View file @f048480
  1 +/***************************************************************************
  2 + * GWcsMOL.hpp - Mollweide's projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 GWcsMOL.hpp
  23 + * @brief Mollweide's projection class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GWCSMOL_HPP
  28 +#define GWCSMOL_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GWcs.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @class GWcsMOL
  37 + *
  38 + * @brief Mollweide's projection class definition
  39 + *
  40 + * This class implements the Mollweide projection for the World Coordinate
  41 + * System.
  42 + ***************************************************************************/
  43 +class GWcsMOL : public GWcs {
  44 +
  45 +public:
  46 + // Constructors and destructors
  47 + GWcsMOL(void);
  48 + GWcsMOL(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 + GWcsMOL(const GWcsMOL& wcs);
  53 + virtual ~GWcsMOL(void);
  54 +
  55 + // Operators
  56 + GWcsMOL& operator=(const GWcsMOL& wcs);
  57 +
  58 + // Implemented pure virtual base class methods
  59 + virtual void clear(void);
  60 + virtual GWcsMOL* 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 GWcsMOL& 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 ("GWcsMOL").
  85 + ***************************************************************************/
  86 +inline
  87 +std::string GWcsMOL::classname(void) const
  88 +{
  89 + return ("GWcsMOL");
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return projection code
  95 + *
  96 + * @return Projection code.
  97 + *
  98 + * Returns the projection code "MOL".
  99 + ***************************************************************************/
  100 +inline
  101 +std::string GWcsMOL::code(void) const
  102 +{
  103 + return "MOL";
  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 GWcsMOL::name(void) const
  116 +{
  117 + return "Mollweide";
  118 +}
  119 +
  120 +#endif /* GWCSMOL_HPP */
... ...
include/GammaLib.hpp
  View file @f048480
... ... @@ -151,6 +151,7 @@
151 151 #include "GWcsAZP.hpp"
152 152 #include "GWcsCAR.hpp"
153 153 #include "GWcsMER.hpp"
  154 +#include "GWcsMOL.hpp"
154 155 #include "GWcsTAN.hpp"
155 156 #include "GWcsSTG.hpp"
156 157  
... ...
include/Makefile.am
  View file @f048480
... ... @@ -93,6 +93,7 @@ pkginclude_HEADERS = GammaLib.hpp \
93 93 GWcsCAR.hpp \
94 94 GWcsTAN.hpp \
95 95 GWcsMER.hpp \
  96 + GWcsMOL.hpp \
96 97 GWcsSTG.hpp \
97 98 GEnergy.hpp \
98 99 GEnergies.hpp \
... ...
pyext/GWcsMOL.i 0 → 100644
  View file @f048480
  1 +/***************************************************************************
  2 + * GWcsMOL.i - Mollweide's projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 GWcsMOL.i
  23 + * @brief Mollweide's projection class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GWcsMOL.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GWcsMOL
  34 + *
  35 + * @brief Mollweide's projection class definition
  36 + ***************************************************************************/
  37 +class GWcsMOL : public GWcs {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GWcsMOL(void);
  42 + GWcsMOL(const std::string& coords,
  43 + const double& crval1, const double& crval2,
  44 + const double& crpix1, const double& crpix2,
  45 + const double& cdelt1, const double& cdelt2);
  46 + GWcsMOL(const GWcsMOL& wcs);
  47 + virtual ~GWcsMOL(void);
  48 +
  49 + // Implemented pure virtual methods
  50 + virtual void clear(void);
  51 + virtual GWcsMOL* clone(void) const;
  52 + virtual std::string classname(void) const;
  53 + virtual std::string code(void) const;
  54 + virtual std::string name(void) const;
  55 +};
  56 +
  57 +
  58 +/***********************************************************************//**
  59 + * @brief GWcsMOL class extension
  60 + ***************************************************************************/
  61 +%extend GWcsMOL {
  62 + GWcsMOL copy() {
  63 + return (*self);
  64 + }
  65 +};
... ...
pyext/gammalib/sky.i
  View file @f048480
1 1 /***************************************************************************
2 2 * sky.i - Sky module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Jurgen Knodlseder *
  4 + * copyright (C) 2011-2015 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -113,5 +113,6 @@
113 113 %include "GWcsAZP.i"
114 114 %include "GWcsCAR.i"
115 115 %include "GWcsMER.i"
  116 +%include "GWcsMOL.i"
116 117 %include "GWcsTAN.i"
117 118 %include "GWcsSTG.i"
... ...
src/sky/GWcsMOL.cpp 0 → 100644
  View file @f048480
  1 +/***************************************************************************
  2 + * GWcsMOL.cpp - Mollweide's projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 GWcsMOL.cpp
  23 + * @brief Mollweide's 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 "GWcsMOL.hpp"
  35 +#include "GWcsRegistry.hpp"
  36 +
  37 +/* __ Method name definitions ____________________________________________ */
  38 +#define G_PRJ_X2S "GWcsMOL::prj_x2s(int, int, int, int, double*, double*,"\
  39 + " double*, double*, int*)"
  40 +#define G_PRJ_S2X "GWcsMOL::prj_s2x(int, int, int, int, double*, double*,"\
  41 + " double*, double*, int*)"
  42 +
  43 +/* __ Macros _____________________________________________________________ */
  44 +#define copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
  45 +
  46 +/* __ Coding definitions _________________________________________________ */
  47 +
  48 +/* __ Debug definitions __________________________________________________ */
  49 +//#define G_DEBUG_PRJ //!< Debug GWcsMOL::prj_x2s
  50 +
  51 +/* __ Local prototypes ___________________________________________________ */
  52 +
  53 +/* __ Constants __________________________________________________________ */
  54 +
  55 +/* __ Globals ____________________________________________________________ */
  56 +const GWcsMOL g_wcs_mol_seed;
  57 +const GWcsRegistry g_wcs_mol_registry(&g_wcs_mol_seed);
  58 +
  59 +
  60 +/*==========================================================================
  61 + = =
  62 + = Constructors/destructors =
  63 + = =
  64 + ==========================================================================*/
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Void constructor
  68 + ***************************************************************************/
  69 +GWcsMOL::GWcsMOL(void) : GWcs()
  70 +{
  71 + // Initialise class members
  72 + init_members();
  73 +
  74 + // Return
  75 + return;
  76 +}
  77 +
  78 +
  79 +/***********************************************************************//**
  80 + * @brief Constructor
  81 + *
  82 + * @param[in] coords Coordinate system.
  83 + * @param[in] crval1 X value of reference pixel.
  84 + * @param[in] crval2 Y value of reference pixel.
  85 + * @param[in] crpix1 X index of reference pixel (starting from 1).
  86 + * @param[in] crpix2 Y index of reference pixel (starting from 1).
  87 + * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
  88 + * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
  89 + ***************************************************************************/
  90 +GWcsMOL::GWcsMOL(const std::string& coords,
  91 + const double& crval1, const double& crval2,
  92 + const double& crpix1, const double& crpix2,
  93 + const double& cdelt1, const double& cdelt2) :
  94 + GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
  95 +
  96 +{
  97 + // Initialise class members
  98 + init_members();
  99 +
  100 + // Return
  101 + return;
  102 +}
  103 +
  104 +
  105 +/***********************************************************************//**
  106 + * @brief Copy constructor
  107 + *
  108 + * @param[in] wcs World Coordinate System.
  109 + ***************************************************************************/
  110 +GWcsMOL::GWcsMOL(const GWcsMOL& wcs) : GWcs(wcs)
  111 +{
  112 + // Initialise class members for clean destruction
  113 + init_members();
  114 +
  115 + // Copy members
  116 + copy_members(wcs);
  117 +
  118 + // Return
  119 + return;
  120 +}
  121 +
  122 +
  123 +/***********************************************************************//**
  124 + * @brief Destructor
  125 + ***************************************************************************/
  126 +GWcsMOL::~GWcsMOL(void)
  127 +{
  128 + // Free members
  129 + free_members();
  130 +
  131 + // Return
  132 + return;
  133 +}
  134 +
  135 +
  136 +/*==========================================================================
  137 + = =
  138 + = Operators =
  139 + = =
  140 + ==========================================================================*/
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Assignment operator
  144 + *
  145 + * @param[in] wcs World Coordinate System.
  146 + * @return World Coordinate System.
  147 + ***************************************************************************/
  148 +GWcsMOL& GWcsMOL::operator=(const GWcsMOL& wcs)
  149 +{
  150 + // Execute only if object is not identical
  151 + if (this != &wcs) {
  152 +
  153 + // Copy base class members
  154 + this->GWcs::operator=(wcs);
  155 +
  156 + // Free members
  157 + free_members();
  158 +
  159 + // Initialise private members for clean destruction
  160 + init_members();
  161 +
  162 + // Copy members
  163 + copy_members(wcs);
  164 +
  165 + } // endif: object was not identical
  166 +
  167 + // Return this object
  168 + return *this;
  169 +}
  170 +
  171 +
  172 +/*==========================================================================
  173 + = =
  174 + = Public methods =
  175 + = =
  176 + ==========================================================================*/
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Clear instance
  180 + *
  181 + * This method properly resets the object to an initial state.
  182 + ***************************************************************************/
  183 +void GWcsMOL::clear(void)
  184 +{
  185 + // Free class members (base and derived classes, derived class first)
  186 + free_members();
  187 + this->GWcs::free_members();
  188 + this->GSkyProjection::free_members();
  189 +
  190 + // Initialise members
  191 + this->GSkyProjection::init_members();
  192 + this->GWcs::init_members();
  193 + init_members();
  194 +
  195 + // Return
  196 + return;
  197 +}
  198 +
  199 +
  200 +/***********************************************************************//**
  201 + * @brief Clone instance
  202 + *
  203 + * @return Pointer to deep copy of Mollweide projection.
  204 + ***************************************************************************/
  205 +GWcsMOL* GWcsMOL::clone(void) const
  206 +{
  207 + return new GWcsMOL(*this);
  208 +}
  209 +
  210 +
  211 +
  212 +/***********************************************************************//**
  213 + * @brief Print WCS information
  214 + *
  215 + * @param[in] chatter Chattiness (defaults to NORMAL).
  216 + * @return String containing WCS information.
  217 + ***************************************************************************/
  218 +std::string GWcsMOL::print(const GChatter& chatter) const
  219 +{
  220 + // Initialise result string
  221 + std::string result;
  222 +
  223 + // Continue only if chatter is not silent
  224 + if (chatter != SILENT) {
  225 +
  226 + // Append header
  227 + result.append("=== GWcsMOL ===");
  228 +
  229 + // Append information
  230 + result.append(wcs_print(chatter));
  231 +
  232 + } // endif: chatter was not silent
  233 +
  234 + // Return result
  235 + return result;
  236 +}
  237 +
  238 +
  239 +/*==========================================================================
  240 + = =
  241 + = Private methods =
  242 + = =
  243 + ==========================================================================*/
  244 +
  245 +/***********************************************************************//**
  246 + * @brief Initialise class members
  247 + *
  248 + * This method sets up the World Coordinate System by calling wcs_set().
  249 + ***************************************************************************/
  250 +void GWcsMOL::init_members(void)
  251 +{
  252 + // Setup World Coordinate System
  253 + wcs_set();
  254 +
  255 + // Return
  256 + return;
  257 +}
  258 +
  259 +
  260 +/***********************************************************************//**
  261 + * @brief Copy class members
  262 + *
  263 + * @param[in] wcs World Coordinate System.
  264 + ***************************************************************************/
  265 +void GWcsMOL::copy_members(const GWcsMOL& wcs)
  266 +{
  267 + // Return
  268 + return;
  269 +}
  270 +
  271 +
  272 +/***********************************************************************//**
  273 + * @brief Delete class members
  274 + ***************************************************************************/
  275 +void GWcsMOL::free_members(void)
  276 +{
  277 + // Return
  278 + return;
  279 +}
  280 +
  281 +
  282 +/***********************************************************************//**
  283 + * @brief Setup of projection
  284 + *
  285 + * This method sets up the projection information. The method has been
  286 + * adapted from the wcslib function prj.c::molset (version 4.20).
  287 + *
  288 + * Given and/or returned:
  289 + * m_r0 Reset to 180/pi if 0.
  290 + * m_phi0 Reset to 0.0 if undefined.
  291 + * m_theta0 Reset to 0.0 if undefined.
  292 + *
  293 + * Returned:
  294 + * m_x0 Fiducial offset in x.
  295 + * m_y0 Fiducial offset in y.
  296 + * m_w[0] sqrt(2)*r0
  297 + * m_w[1] sqrt(2)*r0/90
  298 + * m_w[2] 1/(sqrt(2)*r0)
  299 + * m_w[3] 90/r0
  300 + * m_w[4] 2/pi
  301 + ***************************************************************************/
  302 +void GWcsMOL::prj_set(void) const
  303 +{
  304 + // Initialise projection parameters
  305 + m_w.assign(5, 0.0);
  306 +
  307 + // Set undefined parameters
  308 + if (m_r0 == 0.0) {
  309 + m_r0 = gammalib::rad2deg;
  310 + }
  311 +
  312 + // Precompute
  313 + m_w[0] = gammalib::sqrt_two * m_r0;
  314 + m_w[1] = m_w[0]/90.0;
  315 + m_w[2] = 1.0/m_w[0];
  316 + m_w[3] = 90.0/m_r0;
  317 + m_w[4] = 2.0/gammalib::pi;
  318 +
  319 + // Compute fiducial offset
  320 + prj_off(0.0, 0.0);
  321 +
  322 + // Signal that projection has been set
  323 + m_prjset = true;
  324 +
  325 + // Return
  326 + return;
  327 +}
  328 +
  329 +
  330 +/***********************************************************************//**
  331 + * @brief Pixel-to-spherical deprojection
  332 + *
  333 + * @param[in] nx X vector length.
  334 + * @param[in] ny Y vector length (0=no replication).
  335 + * @param[in] sxy Input vector step.
  336 + * @param[in] spt Output vector step.
  337 + * @param[in] x Vector of projected x coordinates.
  338 + * @param[in] y Vector of projected y coordinates.
  339 + * @param[out] phi Longitude of the projected point in native spherical
  340 + * coordinates [deg].
  341 + * @param[out] theta Latitude of the projected point in native spherical
  342 + * coordinates [deg].
  343 + * @param[out] stat Status return value for each vector element
  344 + *
  345 + * @exception GException::wcs_invalid_x_y
  346 + * One or more of the (x,y) coordinates were invalid, as indicated
  347 + * by the stat vector.
  348 + *
  349 + * Deproject pixel (x,y) coordinates in the plane of projection to native
  350 + * spherical coordinates (phi,theta).
  351 + *
  352 + * This method has been adapted from the wcslib function prj.c::molx2s()
  353 + * (version 4.20). The interface follows very closely that of wcslib. In
  354 + * contrast to the wcslib routine, however, the method assumes that the
  355 + * projection has been setup previously (as this will be done by the
  356 + * constructor).
  357 + ***************************************************************************/
  358 +void GWcsMOL::prj_x2s(int nx, int ny, int sxy, int spt,
  359 + const double* x, const double* y,
  360 + double* phi, double* theta, int* stat) const
  361 +{
  362 + // Set tolerance
  363 + const double tol = 1.0e-12;
  364 +
  365 + // Initialize projection if required
  366 + if (!m_prjset) {
  367 + prj_set();
  368 + }
  369 +
  370 + // Set value replication length mx,my
  371 + int mx;
  372 + int my;
  373 + if (ny > 0) {
  374 + mx = nx;
  375 + my = ny;
  376 + }
  377 + else {
  378 + mx = 1;
  379 + my = 1;
  380 + ny = nx;
  381 + }
  382 +
  383 + // Initialise status code and statistics
  384 + int status = 0;
  385 + int n_invalid = 0;
  386 +
  387 + // Do x dependence
  388 + const double* xp = x;
  389 + int rowoff = 0;
  390 + int rowlen = nx * spt;
  391 + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
  392 + double xj = *xp + m_x0;
  393 + double s = m_w[3] * xj;
  394 + double t = std::abs(xj) - tol;
  395 + double* phip = phi + rowoff;
  396 + double* thetap = theta + rowoff;
  397 + for (int iy = 0; iy < my; ++iy, phip += rowlen, thetap += rowlen) {
  398 + *phip = s;
  399 + *thetap = t;
  400 + }
  401 + }
  402 +
  403 + // Do y dependence
  404 + const double* yp = y;
  405 + double* phip = phi;
  406 + double* thetap = theta;
  407 + int* statp = stat;
  408 + for (int iy = 0; iy < ny; ++iy, yp += sxy) {
  409 +
  410 + // Initialise status
  411 + int istat = 0;
  412 +
  413 + // Compute some stuff
  414 + double yj = *yp + m_y0;
  415 + double y0 = yj / m_r0;
  416 + double r = 2.0 - y0 * y0;
  417 + double s = 0.0;
  418 +
  419 + // Check limits
  420 + if (r <= tol) {
  421 + if (r < -tol) {
  422 + istat = 1;
  423 + status = 3;
  424 + n_invalid++;
  425 + #if defined(G_DEBUG_PRJ)
  426 + std::cout << "prj_x2s(Phi)..:";
  427 + std::cout << " nx=" << nx;
  428 + std::cout << " ny=" << ny;
  429 + std::cout << " iy=" << iy;
  430 + std::cout << " yp=" << *yp;
  431 + std::cout << " yj=" << yj;
  432 + std::cout << " y0=" << y0;
  433 + std::cout << " phip=" << *phip;
  434 + std::cout << " r=" << r;
  435 + std::cout << " s=" << s;
  436 + std::cout << std::endl;
  437 + #endif
  438 + }
  439 + else {
  440 + istat = -1; // OK if fabs(x) < tol whence phi = 0.0
  441 + }
  442 + r = 0.0;
  443 + s = 0.0;
  444 + }
  445 + else {
  446 + r = std::sqrt(r);
  447 + s = 1.0 / r;
  448 + }
  449 +
  450 + // Compute z and t
  451 + double z = yj * m_w[2];
  452 + if (std::abs(z) > 1.0) {
  453 + if (std::abs(z) > 1.0+tol) {
  454 + istat = 1;
  455 + status = 3;
  456 + n_invalid++;
  457 + #if defined(G_DEBUG_PRJ)
  458 + std::cout << "prj_x2s(Phi)..:";
  459 + std::cout << " nx=" << nx;
  460 + std::cout << " ny=" << ny;
  461 + std::cout << " iy=" << iy;
  462 + std::cout << " yp=" << *yp;
  463 + std::cout << " yj=" << yj;
  464 + std::cout << " y0=" << y0;
  465 + std::cout << " phip=" << *phip;
  466 + std::cout << " r=" << r;
  467 + std::cout << " s=" << s;
  468 + std::cout << " z=" << z;
  469 + std::cout << std::endl;
  470 + #endif
  471 + z = 0.0;
  472 + }
  473 + else {
  474 + z = copysign(1.0, z) + y0 * r / gammalib::pi;
  475 + }
  476 + }
  477 + else {
  478 + z = std::asin(z) * m_w[4] + y0 * r / gammalib::pi;
  479 + }
  480 + if (std::abs(z) > 1.0) {
  481 + if (std::abs(z) > 1.0+tol) {
  482 + istat = 1;
  483 + status = 3;
  484 + n_invalid++;
  485 + #if defined(G_DEBUG_PRJ)
  486 + std::cout << "prj_x2s(Phi)..:";
  487 + std::cout << " nx=" << nx;
  488 + std::cout << " ny=" << ny;
  489 + std::cout << " iy=" << iy;
  490 + std::cout << " yp=" << *yp;
  491 + std::cout << " yj=" << yj;
  492 + std::cout << " y0=" << y0;
  493 + std::cout << " phip=" << *phip;
  494 + std::cout << " r=" << r;
  495 + std::cout << " s=" << s;
  496 + std::cout << " z=" << z;
  497 + std::cout << std::endl;
  498 + #endif
  499 + z = 0.0;
  500 + }
  501 + else {
  502 + z = copysign(1.0, z);
  503 + }
  504 + }
  505 + double t = std::asin(z) * gammalib::rad2deg;
  506 +
  507 + // Do ...
  508 + for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt, statp += spt) {
  509 +
  510 + // Set status flag
  511 + if (istat < 0 && *thetap >= 0.0) {
  512 + *statp = 1;
  513 + status = 3;
  514 + n_invalid++;
  515 + }
  516 + else {
  517 + *statp = 0;
  518 + }
  519 +
  520 + // Set values
  521 + *phip *= s;
  522 + *thetap = t;
  523 +
  524 + }
  525 +
  526 + } // endfor: did y dependence
  527 +
  528 + // Handle status code
  529 + if (status == 3) {
  530 + throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
  531 + }
  532 +
  533 + // Return
  534 + return;
  535 +}
  536 +
  537 +
  538 +/***********************************************************************//**
  539 + * @brief Generic spherical-to-pixel projection
  540 + *
  541 + * @param[in] nphi Longitude vector length.
  542 + * @param[in] ntheta Latitude vector length (0=no replication).
  543 + * @param[in] spt Input vector step.
  544 + * @param[in] sxy Output vector step.
  545 + * @param[in] phi Longitude vector of the projected point in native spherical
  546 + * coordinates [deg].
  547 + * @param[in] theta Latitude vector of the projected point in native spherical
  548 + * coordinates [deg].
  549 + * @param[out] x Vector of projected x coordinates.
  550 + * @param[out] y Vector of projected y coordinates.
  551 + * @param[out] stat Status return value for each vector element (always 0)
  552 + *
  553 + * @exception GException::wcs_invalid_phi_theta
  554 + * One or more of the (phi,theta) coordinates were invalid, as
  555 + * indicated by the stat vector.
  556 + *
  557 + * Project native spherical coordinates (phi,theta) to pixel (x,y)
  558 + * coordinates in the plane of projection.
  559 + *
  560 + * This method has been adapted from the wcslib function prj.c::mols2x()
  561 + * (version 4.20). The interface follows very closely that of wcslib. In
  562 + * contrast to the wcslib routine, however, the method assumes that the
  563 + * projection has been setup previously (as this will be done by the
  564 + * constructor).
  565 + ***************************************************************************/
  566 +void GWcsMOL::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  567 + const double* phi, const double* theta,
  568 + double* x, double* y, int* stat) const
  569 +{
  570 + // Set tolerance
  571 + const double tol = 1.0e-13;
  572 +
  573 + // Initialize projection if required
  574 + if (!m_prjset) {
  575 + prj_set();
  576 + }
  577 +
  578 + // Set value replication length mphi,mtheta
  579 + int mphi;
  580 + int mtheta;
  581 + if (ntheta > 0) {
  582 + mphi = nphi;
  583 + mtheta = ntheta;
  584 + }
  585 + else {
  586 + mphi = 1;
  587 + mtheta = 1;
  588 + ntheta = nphi;
  589 + }
  590 +
  591 + // Initialise status code and statistics
  592 + int status = 0;
  593 + int n_invalid = 0;
  594 +
  595 + // Do phi dependence
  596 + const double* phip = phi;
  597 + int rowoff = 0;
  598 + int rowlen = nphi * sxy;
  599 + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
  600 + double xi = m_w[1] * (*phip);
  601 + double* xp = x + rowoff;
  602 + for (int itheta = 0; itheta < mtheta; ++itheta) {
  603 + *xp = xi;
  604 + xp += rowlen;
  605 + }
  606 + }
  607 +
  608 + // Do theta dependence
  609 + const double* thetap = theta;
  610 + double* xp = x;
  611 + double* yp = y;
  612 + int* statp = stat;
  613 + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
  614 +
  615 + // Compute xi and eta
  616 + double xi = 0.0;
  617 + double eta = 0.0;
  618 + if (std::abs(*thetap) == 90.0) {
  619 + xi = 0.0;
  620 + eta = copysign(m_w[0], *thetap);
  621 + }
  622 + else if (*thetap == 0.0) {
  623 + xi = 1.0;
  624 + eta = 0.0;
  625 + }
  626 + else {
  627 + double u = gammalib::pi * std::sin(*thetap * gammalib::deg2rad);
  628 + double v0 = -gammalib::pi;
  629 + double v1 = gammalib::pi;
  630 + double v = u;
  631 + for (int k = 0; k < 100; ++k) {
  632 + double resid = (v - u) + std::sin(v);
  633 + if (resid < 0.0) {
  634 + if (resid > -tol) {
  635 + break;
  636 + }
  637 + v0 = v;
  638 + }
  639 + else {
  640 + if (resid < tol) {
  641 + break;
  642 + }
  643 + v1 = v;
  644 + }
  645 + v = 0.5 * (v0 + v1);
  646 + }
  647 + double gamma = 0.5 * v;
  648 + xi = std::cos(gamma);
  649 + eta = std::sin(gamma) * m_w[0];
  650 + }
  651 + eta -= m_y0;
  652 +
  653 + // Do ...
  654 + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy, statp += sxy) {
  655 + *xp = xi * (*xp) - m_x0;
  656 + *yp = eta;
  657 + *statp = 0;
  658 + }
  659 +
  660 + } // endfor: theta dependence
  661 +
  662 + // Handle status code
  663 + if (status == 4) {
  664 + throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
  665 + }
  666 +
  667 + // Return
  668 + return;
  669 +}
... ...
src/sky/Makefile.am
  View file @f048480
... ... @@ -21,6 +21,7 @@ sources = GSkyDir.cpp \
21 21 GWcsAZP.cpp \
22 22 GWcsCAR.cpp \
23 23 GWcsMER.cpp \
  24 + GWcsMOL.cpp \
24 25 GWcsTAN.cpp \
25 26 GWcsSTG.cpp \
26 27 GException_sky.cpp
... ...
test/test_GSky.cpp
  View file @f048480
... ... @@ -112,7 +112,7 @@ double TestGSky::wcs_forth_back_pixel(GWcs* wcs, int nx, int ny, double&amp; crpix1,
112 112 for (int iy = -ny; iy <= ny; ++iy) {
113 113  
114 114 // Set y value
115   - double y = double(iy) + crpix2;
  115 + double y = double(iy) + crpix2;
116 116  
117 117 // Skip pixels outside valid range (they are not expected to
118 118 // transform bijectively)
... ... @@ -138,7 +138,8 @@ double TestGSky::wcs_forth_back_pixel(GWcs* wcs, int nx, int ny, double&amp; crpix1,
138 138 #if defined(G_WCS_FORTH_BACK_PIXEL_DEBUG)
139 139 if (dist > 0.001) {
140 140 std::cout << std::endl;
141   - std::cout << "dist=" << dist;
  141 + std::cout << wcs->code() << ":";
  142 + std::cout << " dist=" << dist;
142 143 std::cout << " dx=" << dx;
143 144 std::cout << " dy=" << dy;
144 145 std::cout << " pix_in=" << pix_in;
... ...
test/test_GSky.py
  View file @f048480
1 1 # ==========================================================================
2 2 # This module performs unit tests for the GammaLib sky module.
3 3 #
4   -# Copyright (C) 2012-2014 Juergen Knoedlseder
  4 +# Copyright (C) 2012-2015 Juergen Knoedlseder
5 5 #
6 6 # This program is free software: you can redistribute it and/or modify
7 7 # it under the terms of the GNU General Public License as published by
... ... @@ -76,6 +76,7 @@ class Test(GPythonTestSuite):
76 76 self.append(self.test_skymap_azp, "Test AZP projection map")
77 77 self.append(self.test_skymap_car, "Test CAR projection map")
78 78 self.append(self.test_skymap_mer, "Test MER projection map")
  79 + self.append(self.test_skymap_mol, "Test MOL projection map")
79 80 self.append(self.test_skymap_stg, "Test STG projection map")
80 81 self.append(self.test_skymap_tan, "Test TAN projection map")
81 82 self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
... ... @@ -232,6 +233,17 @@ class Test(GPythonTestSuite):
232 233 # Return
233 234 return
234 235  
  236 + # Test MOL projection
  237 + def test_skymap_mol(self):
  238 + """
  239 + Test MOL projection.
  240 + """
  241 + # Execute generic test
  242 + self.test_skymap_proj("MOL")
  243 +
  244 + # Return
  245 + return
  246 +
235 247 # Test STG projection
236 248 def test_skymap_stg(self):
237 249 """
... ...