Thierry Louge / gammalib

Commit 12d0491734bf02b70083400cf9ce150c021749b0

Authored by Jürgen Knödlseder
1 parent b447670c

Add SIN WCS projection (#1798) 

Added GWcsSIN class.
Inline Side-by-side
Showing 11 changed files with 897 additions and 8 deletions
ChangeLog
  View file @12d0491
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Add GWcsSIN projection class (#1798)
6 7 CTA cube response no longer selects on RoI and energies (#1789)
7 8 Add weights to CTA event cube and event bin (#1789)
8 9 Remove GCTACubeExposure, GCTACubePsf and GCTACubeBackground elogmeans() methods
... ...
... ... @@ -15,6 +15,7 @@ made since the last release of GammaLib.
15 15 ---------------------------------
16 16 The following classes have been added:
17 17 - GCTAPsfTable
  18 +- GWcsSIN
18 19  
19 20 The following classes have been removed:
20 21 - GCTAOnOffObservations
... ... @@ -263,6 +264,9 @@ longitude, latitude and a number of maps has a dimension of 1 (#1733).
263 264 An optional extname argument has been added to the GSkyMap::save() method
264 265 so that the extension name can be set (#1787).
265 266  
  267 +The Orthographic/synthesis (SIN) projection has been added to the WCS
  268 +projections (#1798).
  269 +
266 270  
267 271 13. Support module
268 272 ------------------
... ...
include/GWcsSIN.hpp 0 → 100644
  View file @12d0491
  1 +/***************************************************************************
  2 + * GWcsSIN.hpp - Orthographic/synthesis (SIN) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 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 GWcsSIN.hpp
  23 + * @brief Orthographic/synthesis (SIN) projection class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GWcsSIN_HPP
  28 +#define GWcsSIN_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GWcs.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @class GWcsSIN
  37 + *
  38 + * @brief Orthographic/synthesis (SIN) projection class definition
  39 + *
  40 + * This class implements the "orthographic/synthesis" projection for the
  41 + * World Coordinate System.
  42 + ***************************************************************************/
  43 +class GWcsSIN : public GWcs {
  44 +
  45 +public:
  46 + // Constructors and destructors
  47 + GWcsSIN(void);
  48 + GWcsSIN(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 + GWcsSIN(const GWcsSIN& wcs);
  53 + virtual ~GWcsSIN(void);
  54 +
  55 + // Operators
  56 + GWcsSIN& operator=(const GWcsSIN& wcs);
  57 +
  58 + // Implemented pure virtual base class methods
  59 + virtual void clear(void);
  60 + virtual GWcsSIN* 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 GWcsSIN& 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 ("GWcsSIN").
  85 + ***************************************************************************/
  86 +inline
  87 +std::string GWcsSIN::classname(void) const
  88 +{
  89 + return ("GWcsSIN");
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return projection code
  95 + *
  96 + * @return Projection code.
  97 + *
  98 + * Returns the projection code "SIN".
  99 + ***************************************************************************/
  100 +inline
  101 +std::string GWcsSIN::code(void) const
  102 +{
  103 + return "SIN";
  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 GWcsSIN::name(void) const
  116 +{
  117 + return "orthographic/synthesis";
  118 +}
  119 +
  120 +#endif /* GWcsSIN_HPP */
... ...
include/GammaLib.hpp
  View file @12d0491
1 1 /***************************************************************************
2 2 * GammaLib.hpp - Gamma-Ray Astronomy Library Header file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2006-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -154,8 +154,9 @@
154 154 #include "GWcsCAR.hpp"
155 155 #include "GWcsMER.hpp"
156 156 #include "GWcsMOL.hpp"
157   -#include "GWcsTAN.hpp"
  157 +#include "GWcsSIN.hpp"
158 158 #include "GWcsSTG.hpp"
  159 +#include "GWcsTAN.hpp"
159 160  
160 161 /* __ Observation handling _______________________________________________ */
161 162 #include "GEnergy.hpp"
... ...
include/Makefile.am
  View file @12d0491
... ... @@ -93,10 +93,11 @@ pkginclude_HEADERS = GammaLib.hpp \
93 93 GWcsAIT.hpp \
94 94 GWcsAZP.hpp \
95 95 GWcsCAR.hpp \
96   - GWcsTAN.hpp \
97 96 GWcsMER.hpp \
98 97 GWcsMOL.hpp \
  98 + GWcsSIN.hpp \
99 99 GWcsSTG.hpp \
  100 + GWcsTAN.hpp \
100 101 GEnergy.hpp \
101 102 GEnergies.hpp \
102 103 GTime.hpp \
... ...
pyext/GWcsSIN.i 0 → 100644
  View file @12d0491
  1 +/***************************************************************************
  2 + * GWcsSIN.i - Orthographic/synthesis (SIN) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 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 GWcsSIN.i
  23 + * @brief Orthographic/synthesis (SIN) projection class Python interface
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GWcsSIN.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GWcsSIN
  34 + *
  35 + * @brief Orthographic/synthesis (SIN) projection class definition
  36 + ***************************************************************************/
  37 +class GWcsSIN : public GWcs {
  38 +public:
  39 + // Constructors and destructors
  40 + GWcsSIN(void);
  41 + GWcsSIN(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 + GWcsSIN(const GWcsSIN& wcs);
  46 + virtual ~GWcsSIN(void);
  47 +
  48 + // Implemented pure virtual methods
  49 + virtual void clear(void);
  50 + virtual GWcsSIN* 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 GWcsSIN class extension
  59 + ***************************************************************************/
  60 +%extend GWcsSIN {
  61 + GWcsSIN copy() {
  62 + return (*self);
  63 + }
  64 +};
... ...
pyext/gammalib/sky.i
  View file @12d0491
... ... @@ -124,5 +124,6 @@ namespace std {
124 124 %include "GWcsCAR.i"
125 125 %include "GWcsMER.i"
126 126 %include "GWcsMOL.i"
127   -%include "GWcsTAN.i"
  127 +%include "GWcsSIN.i"
128 128 %include "GWcsSTG.i"
  129 +%include "GWcsTAN.i"
... ...
src/sky/GWcsSIN.cpp 0 → 100644
  View file @12d0491
  1 +/***************************************************************************
  2 + * GWcsSIN.cpp - Orthographic/synthesis (SIN) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 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 GWcsSIN.cpp
  23 + * @brief Orthographic/synthesis (SIN) 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 "GWcsSIN.hpp"
  34 +#include "GWcsRegistry.hpp"
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_PRJ_X2S "GWcsSIN::prj_x2s(int, int, int, int, double*, double*,"\
  38 + " double*, double*, int*)"
  39 +#define G_PRJ_S2X "GWcsSIN::prj_s2x(int, int, int, int, double*, double*,"\
  40 + " double*, double*, int*)"
  41 +
  42 +/* __ Macros _____________________________________________________________ */
  43 +
  44 +/* __ Coding definitions _________________________________________________ */
  45 +
  46 +/* __ Debug definitions __________________________________________________ */
  47 +
  48 +/* __ Local prototypes ___________________________________________________ */
  49 +
  50 +/* __ Constants __________________________________________________________ */
  51 +
  52 +/* __ Globals ____________________________________________________________ */
  53 +const GWcsSIN g_wcs_sin_seed;
  54 +const GWcsRegistry g_wcs_sin_registry(&g_wcs_sin_seed);
  55 +
  56 +
  57 +/*==========================================================================
  58 + = =
  59 + = Constructors/destructors =
  60 + = =
  61 + ==========================================================================*/
  62 +
  63 +/***********************************************************************//**
  64 + * @brief Void constructor
  65 + ***************************************************************************/
  66 +GWcsSIN::GWcsSIN(void) : GWcs()
  67 +{
  68 + // Initialise class members
  69 + init_members();
  70 +
  71 + // Return
  72 + return;
  73 +}
  74 +
  75 +
  76 +/***********************************************************************//**
  77 + * @brief Projection constructor
  78 + *
  79 + * @param[in] coords Coordinate system.
  80 + * @param[in] crval1 X value of reference pixel.
  81 + * @param[in] crval2 Y value of reference pixel.
  82 + * @param[in] crpix1 X index of reference pixel (starting from 1).
  83 + * @param[in] crpix2 Y index of reference pixel (starting from 1).
  84 + * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
  85 + * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
  86 + ***************************************************************************/
  87 +GWcsSIN::GWcsSIN(const std::string& coords,
  88 + const double& crval1, const double& crval2,
  89 + const double& crpix1, const double& crpix2,
  90 + const double& cdelt1, const double& cdelt2) :
  91 + GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
  92 +
  93 +{
  94 + // Initialise class members
  95 + init_members();
  96 +
  97 + // Return
  98 + return;
  99 +}
  100 +
  101 +
  102 +/***********************************************************************//**
  103 + * @brief Copy constructor
  104 + *
  105 + * @param[in] wcs World Coordinate System.
  106 + ***************************************************************************/
  107 +GWcsSIN::GWcsSIN(const GWcsSIN& wcs) : GWcs(wcs)
  108 +{
  109 + // Initialise class members for clean destruction
  110 + init_members();
  111 +
  112 + // Copy members
  113 + copy_members(wcs);
  114 +
  115 + // Return
  116 + return;
  117 +}
  118 +
  119 +
  120 +/***********************************************************************//**
  121 + * @brief Destructor
  122 + ***************************************************************************/
  123 +GWcsSIN::~GWcsSIN(void)
  124 +{
  125 + // Free members
  126 + free_members();
  127 +
  128 + // Return
  129 + return;
  130 +}
  131 +
  132 +
  133 +/*==========================================================================
  134 + = =
  135 + = Operators =
  136 + = =
  137 + ==========================================================================*/
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Assignment operator
  141 + *
  142 + * @param[in] wcs World Coordinate System.
  143 + * @return World Coordinate System.
  144 + ***************************************************************************/
  145 +GWcsSIN& GWcsSIN::operator=(const GWcsSIN& wcs)
  146 +{
  147 + // Execute only if object is not identical
  148 + if (this != &wcs) {
  149 +
  150 + // Copy base class members
  151 + this->GWcs::operator=(wcs);
  152 +
  153 + // Free members
  154 + free_members();
  155 +
  156 + // Initialise private members for clean destruction
  157 + init_members();
  158 +
  159 + // Copy members
  160 + copy_members(wcs);
  161 +
  162 + } // endif: object was not identical
  163 +
  164 + // Return this object
  165 + return *this;
  166 +}
  167 +
  168 +
  169 +/*==========================================================================
  170 + = =
  171 + = Public methods =
  172 + = =
  173 + ==========================================================================*/
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Clear projection
  177 + *
  178 + * This method properly resets the object to an initial state.
  179 + ***************************************************************************/
  180 +void GWcsSIN::clear(void)
  181 +{
  182 + // Free class members (base and derived classes, derived class first)
  183 + free_members();
  184 + this->GWcs::free_members();
  185 + this->GSkyProjection::free_members();
  186 +
  187 + // Initialise members
  188 + this->GSkyProjection::init_members();
  189 + this->GWcs::init_members();
  190 + init_members();
  191 +
  192 + // Return
  193 + return;
  194 +}
  195 +
  196 +
  197 +/***********************************************************************//**
  198 + * @brief Clone projection
  199 + *
  200 + * @return Pointer to deep copy of projection.
  201 + ***************************************************************************/
  202 +GWcsSIN* GWcsSIN::clone(void) const
  203 +{
  204 + return new GWcsSIN(*this);
  205 +}
  206 +
  207 +
  208 +/***********************************************************************//**
  209 + * @brief Print projection information
  210 + *
  211 + * @param[in] chatter Chattiness (defaults to NORMAL).
  212 + * @return String containing projection information.
  213 + ***************************************************************************/
  214 +std::string GWcsSIN::print(const GChatter& chatter) const
  215 +{
  216 + // Initialise result string
  217 + std::string result;
  218 +
  219 + // Continue only if chatter is not silent
  220 + if (chatter != SILENT) {
  221 +
  222 + // Append header
  223 + result.append("=== GWcsSIN ===");
  224 +
  225 + // Append information
  226 + result.append(wcs_print(chatter));
  227 +
  228 + } // endif: chatter was not silent
  229 +
  230 + // Return result
  231 + return result;
  232 +}
  233 +
  234 +
  235 +/*==========================================================================
  236 + = =
  237 + = Private methods =
  238 + = =
  239 + ==========================================================================*/
  240 +
  241 +/***********************************************************************//**
  242 + * @brief Initialise class members
  243 + *
  244 + * This method sets up the World Coordinate System by calling wcs_set().
  245 + ***************************************************************************/
  246 +void GWcsSIN::init_members(void)
  247 +{
  248 + // Setup World Coordinate System
  249 + wcs_set();
  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 GWcsSIN::copy_members(const GWcsSIN& wcs)
  262 +{
  263 + // Return
  264 + return;
  265 +}
  266 +
  267 +
  268 +/***********************************************************************//**
  269 + * @brief Delete class members
  270 + ***************************************************************************/
  271 +void GWcsSIN::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::sinset.
  283 + *
  284 + * Given:
  285 + * m_pv[1] Obliqueness parameter xi.
  286 + * m_pv[2] Obliqueness parameter eta.
  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 90.0 if undefined.
  292 + *
  293 + * Returned:
  294 + * m_x0 Fiducial offset in x.
  295 + * m_y0 Fiducial offset in y.
  296 + * m_w[0] 1/r0
  297 + * m_w[1] xi**2 + eta**2
  298 + * m_w[2] xi**2 + eta**2 + 1
  299 + * m_w[3] xi**2 + eta**2 - 1
  300 + ***************************************************************************/
  301 +void GWcsSIN::prj_set(void) const
  302 +{
  303 + // Initialise projection parameters
  304 + m_w.assign(4,0.0);
  305 +
  306 + // Set undefined parameters
  307 + if (undefined(m_pv[1])) {
  308 + m_pv[1] = 0.0;
  309 + }
  310 + if (undefined(m_pv[2])) {
  311 + m_pv[2] = 0.0;
  312 + }
  313 + if (m_r0 == 0.0) {
  314 + m_r0 = gammalib::rad2deg;
  315 + }
  316 +
  317 + // Precompute
  318 + m_w[0] = 1.0 / m_r0;
  319 + m_w[1] = m_pv[1]*m_pv[1] + m_pv[2]*m_pv[2];
  320 + m_w[2] = m_w[1] + 1.0;
  321 + m_w[3] = m_w[1] - 1.0;
  322 +
  323 + // Compute fiducial offset
  324 + prj_off(0.0, 90.0);
  325 +
  326 + // Signal that projection has been set
  327 + m_prjset = true;
  328 +
  329 + // Return
  330 + return;
  331 +}
  332 +
  333 +
  334 +/***********************************************************************//**
  335 + * @brief Pixel-to-spherical deprojection
  336 + *
  337 + * @param[in] nx X vector length.
  338 + * @param[in] ny Y vector length (0=no replication).
  339 + * @param[in] sxy Input vector step.
  340 + * @param[in] spt Output vector step.
  341 + * @param[in] x Vector of projected x coordinates.
  342 + * @param[in] y Vector of projected y coordinates.
  343 + * @param[out] phi Longitude of the projected point in native spherical
  344 + * coordinates [deg].
  345 + * @param[out] theta Latitude of the projected point in native spherical
  346 + * coordinates [deg].
  347 + * @param[out] stat Status return value for each vector element (always 0)
  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::sinx2s().
  353 + * The interface follows very closely that of wcslib. In contrast to the
  354 + * wcslib routine, however, the method assumes that the projection has been
  355 + * setup previously (as this will be done by the constructor).
  356 + ***************************************************************************/
  357 +void GWcsSIN::prj_x2s(int nx, int ny, int sxy, int spt,
  358 + const double* x, const double* y,
  359 + double* phi, double* theta, int* stat) const
  360 +{
  361 + // Set tolerance
  362 + const double tol = 1.0e-13;
  363 +
  364 + // Initialize projection if required
  365 + if (!m_prjset) {
  366 + prj_set();
  367 + }
  368 +
  369 + // Initialise status code and statistics
  370 + int status = 0;
  371 + int n_invalid = 0;
  372 +
  373 + // Get obliqueness parameters
  374 + double xi = m_pv[1];
  375 + double eta = m_pv[2];
  376 +
  377 + // Set value replication length mx,my
  378 + int mx;
  379 + int my;
  380 + if (ny > 0) {
  381 + mx = nx;
  382 + my = ny;
  383 + }
  384 + else {
  385 + mx = 1;
  386 + my = 1;
  387 + ny = nx;
  388 + }
  389 +
  390 + // Do x dependence
  391 + const double* xp = x;
  392 + int rowoff = 0;
  393 + int rowlen = nx * spt;
  394 + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
  395 + double x0 = (*xp + m_x0) * m_w[0];
  396 + double* phip = phi + rowoff;
  397 + for (int iy = 0; iy < my; ++iy, phip += rowlen) {
  398 + *phip = x0;
  399 + }
  400 + }
  401 +
  402 + // Do y dependence
  403 + const double* yp = y;
  404 + double* phip = phi;
  405 + double* thetap = theta;
  406 + int* statp = stat;
  407 + for (int iy = 0; iy < ny; ++iy, yp += sxy) {
  408 + double y0 = (*yp + m_y0) * m_w[0];
  409 + double y02 = y0*y0;
  410 + for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
  411 +
  412 + // Compute intermediaries
  413 + double x0 = *phip;
  414 + double r2 = x0*x0 + y02;
  415 +
  416 + // Orthographic projection
  417 + if (m_w[1] == 0.0) {
  418 +
  419 + // Compute Phi
  420 + if (r2 != 0.0) {
  421 + *phip = gammalib::atan2d(x0, -y0);
  422 + }
  423 + else {
  424 + *phip = 0.0;
  425 + }
  426 +
  427 + // Compute Theta
  428 + if (r2 < 0.5) {
  429 + *thetap = gammalib::acosd(std::sqrt(r2));
  430 + }
  431 + else if (r2 <= 1.0) {
  432 + *thetap = gammalib::asind(std::sqrt(1.0 - r2));
  433 + }
  434 + else {
  435 + *thetap = 0.0;
  436 + *(statp++) = 1;
  437 + status = 3;
  438 + n_invalid++;
  439 + continue;
  440 + }
  441 +
  442 + } // endif: orthographic projection
  443 +
  444 + // "Synthesis" projection
  445 + else {
  446 +
  447 + // ...
  448 + double z;
  449 + double xy = x0*xi + y0*eta;
  450 +
  451 + // Use small angle formula
  452 + if (r2 < 1.0e-10) {
  453 + z = r2/2.0;
  454 + *thetap = 90.0 - gammalib::rad2deg*std::sqrt(r2/(1.0 + xy));
  455 +
  456 + }
  457 + else {
  458 + double a = m_w[2];
  459 + double b = xy - m_w[1];
  460 + double c = r2 - xy - xy + m_w[3];
  461 + double d = b*b - a*c;
  462 +
  463 + // Check for a solution
  464 + if (d < 0.0) {
  465 + *phip = 0.0;
  466 + *thetap = 0.0;
  467 + *(statp++) = 1;
  468 + status = 3;
  469 + n_invalid++;
  470 + continue;
  471 + }
  472 +
  473 + // ...
  474 + d = std::sqrt(d);
  475 +
  476 + // Choose solution closest to pole
  477 + double sinth1 = (-b + d)/a;
  478 + double sinth2 = (-b - d)/a;
  479 + double sinthe = (sinth1 > sinth2) ? sinth1 : sinth2;
  480 + if (sinthe > 1.0) {
  481 + if (sinthe-1.0 < tol) {
  482 + sinthe = 1.0;
  483 + }
  484 + else {
  485 + sinthe = (sinth1 < sinth2) ? sinth1 : sinth2;
  486 + }
  487 + }
  488 +
  489 + // ...
  490 + if (sinthe < -1.0) {
  491 + if (sinthe+1.0 > -tol) {
  492 + sinthe = -1.0;
  493 + }
  494 + }
  495 +
  496 + // ...
  497 + if (sinthe > 1.0 || sinthe < -1.0) {
  498 + *phip = 0.0;
  499 + *thetap = 0.0;
  500 + *(statp++) = 1;
  501 + status = 3;
  502 + n_invalid++;
  503 + continue;
  504 + }
  505 +
  506 + // ...
  507 + *thetap = gammalib::asind(sinthe);
  508 + z = 1.0 - sinthe;
  509 + }
  510 +
  511 + // ...
  512 + double x1 = -y0 + eta * z;
  513 + double y1 = x0 - xi * z;
  514 + if (x1 == 0.0 && y1 == 0.0) {
  515 + *phip = 0.0;
  516 + }
  517 + else {
  518 + *phip = gammalib::atan2d(y1,x1);
  519 + }
  520 +
  521 + } // endelse: Synthesis projection
  522 +
  523 + // Reset error
  524 + *(statp++) = 0;
  525 +
  526 + } // endfor
  527 +
  528 + } // endfor: do Y dependence
  529 +
  530 + // Handle status code
  531 + if (status == 3) {
  532 + throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
  533 + }
  534 +
  535 + // Return
  536 + return;
  537 +}
  538 +
  539 +
  540 +/***********************************************************************//**
  541 + * @brief Generic spherical-to-pixel projection
  542 + *
  543 + * @param[in] nphi Longitude vector length.
  544 + * @param[in] ntheta Latitude vector length (0=no replication).
  545 + * @param[in] spt Input vector step.
  546 + * @param[in] sxy Output vector step.
  547 + * @param[in] phi Longitude vector of the projected point in native spherical
  548 + * coordinates [deg].
  549 + * @param[in] theta Latitude vector of the projected point in native spherical
  550 + * coordinates [deg].
  551 + * @param[out] x Vector of projected x coordinates.
  552 + * @param[out] y Vector of projected y coordinates.
  553 + * @param[out] stat Status return value for each vector element (always 0)
  554 + *
  555 + * Project native spherical coordinates (phi,theta) to pixel (x,y)
  556 + * coordinates in the plane of projection.
  557 + *
  558 + * This method has been adapted from the wcslib function prj.c::sins2x().
  559 + * The interface follows very closely that of wcslib. In contrast to the
  560 + * wcslib routine, however, the method assumes that the projection has been
  561 + * setup previously (as this will be done by the constructor).
  562 + ***************************************************************************/
  563 +void GWcsSIN::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  564 + const double* phi, const double* theta,
  565 + double* x, double* y, int* stat) const
  566 +{
  567 + // Initialize projection if required
  568 + if (!m_prjset) {
  569 + prj_set();
  570 + }
  571 +
  572 + // Initialise status code and statistics
  573 + int status = 0;
  574 + int n_invalid = 0;
  575 +
  576 + // Set value replication length mphi,mtheta
  577 + int mphi;
  578 + int mtheta;
  579 + if (ntheta > 0) {
  580 + mphi = nphi;
  581 + mtheta = ntheta;
  582 + }
  583 + else {
  584 + mphi = 1;
  585 + mtheta = 1;
  586 + ntheta = nphi;
  587 + }
  588 +
  589 + // Do phi dependence
  590 + const double* phip = phi;
  591 + int rowoff = 0;
  592 + int rowlen = nphi * sxy;
  593 + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
  594 + double sinphi;
  595 + double cosphi;
  596 + gammalib::sincosd(*phip, &sinphi, &cosphi);
  597 + double* xp = x + rowoff;
  598 + double* yp = y + rowoff;
  599 + for (int itheta = 0; itheta < mtheta; ++itheta) {
  600 + *xp = sinphi;
  601 + *yp = cosphi;
  602 + xp += rowlen;
  603 + yp += rowlen;
  604 + }
  605 + }
  606 +
  607 + // Do theta dependence
  608 + const double* thetap = theta;
  609 + double* xp = x;
  610 + double* yp = y;
  611 + int* statp = stat;
  612 + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
  613 +
  614 + // ...
  615 + double z;
  616 + double costhe;
  617 + double t = (90.0 - std::abs(*thetap)) * gammalib::deg2rad;
  618 + if (t < 1.0e-5) {
  619 + if (*thetap > 0.0) {
  620 + z = t * t / 2.0;
  621 + }
  622 + else {
  623 + z = 2.0 - t * t / 2.0;
  624 + }
  625 + costhe = t;
  626 + }
  627 + else {
  628 + z = 1.0 - gammalib::sind(*thetap);
  629 + costhe = gammalib::cosd(*thetap);
  630 + }
  631 + double r = m_r0 * costhe;
  632 +
  633 + // Orthographic projection
  634 + if (m_w[1] == 0.0) {
  635 + int istat = 0;
  636 + if (m_bounds) {
  637 + if (*thetap < 0.0) {
  638 + istat = 1;
  639 + status = 4;
  640 + n_invalid++;
  641 + }
  642 + }
  643 + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
  644 + *xp = r*(*xp) - m_x0;
  645 + *yp = -r*(*yp) - m_y0;
  646 + *(statp++) = istat;
  647 + }
  648 + }
  649 +
  650 + // Synthesis projection
  651 + else {
  652 +
  653 + // ...
  654 + z *= m_r0;
  655 + double z1 = m_pv[1] * z - m_x0;
  656 + double z2 = m_pv[2] * z - m_y0;
  657 +
  658 + // ...
  659 + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
  660 + int istat = 0;
  661 + if (m_bounds) {
  662 + double t = -gammalib::atand(m_pv[1]*(*xp) - m_pv[2]*(*yp));
  663 + if (*thetap < t) {
  664 + istat = 1;
  665 + status = 4;
  666 + n_invalid++;
  667 + }
  668 + }
  669 + *xp = r * (*xp) + z1;
  670 + *yp = -r * (*yp) + z2;
  671 + *(statp++) = istat;
  672 + }
  673 +
  674 + } // endelse: Synthesis projection
  675 +
  676 + } // endfor: do theta dependence
  677 +
  678 + // Handle status code
  679 + if (status == 4) {
  680 + throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
  681 + }
  682 +
  683 + // Return
  684 + return;
  685 +}
... ...
src/sky/Makefile.am
  View file @12d0491
... ... @@ -22,8 +22,9 @@ sources = GSkyDir.cpp \
22 22 GWcsCAR.cpp \
23 23 GWcsMER.cpp \
24 24 GWcsMOL.cpp \
25   - GWcsTAN.cpp \
  25 + GWcsSIN.cpp \
26 26 GWcsSTG.cpp \
  27 + GWcsTAN.cpp \
27 28 GException_sky.cpp
28 29  
29 30 # Build libtool library
... ...
test/test_GSky.cpp
  View file @12d0491
... ... @@ -31,7 +31,6 @@
31 31 #include <iostream> // cout, cerr
32 32 #include <stdexcept> // std::exception
33 33 #include <cstdlib> // for system
34   -//#include <stdlib.h>
35 34 #include "test_GSky.hpp"
36 35 #include "GTools.hpp"
37 36  
... ... @@ -423,8 +422,8 @@ void TestGSky::test_GWcs(void)
423 422 double crval2 = 22.01;
424 423 double crpix1 = 100.5;
425 424 double crpix2 = 100.5;
426   - double cdelt1 = 0.5;
427   - double cdelt2 = 0.5;
  425 + double cdelt1 = 0.1;
  426 + double cdelt2 = 0.1;
428 427 int nx = 200;
429 428 int ny = 200;
430 429  
... ...
test/test_GSky.py
  View file @12d0491
... ... @@ -79,6 +79,7 @@ class Test(gammalib.GPythonTestSuite):
79 79 self.append(self.test_skymap_car, "Test CAR projection map")
80 80 self.append(self.test_skymap_mer, "Test MER projection map")
81 81 self.append(self.test_skymap_mol, "Test MOL projection map")
  82 + self.append(self.test_skymap_sin, "Test SIN projection map")
82 83 self.append(self.test_skymap_stg, "Test STG projection map")
83 84 self.append(self.test_skymap_tan, "Test TAN projection map")
84 85 self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
... ... @@ -343,6 +344,17 @@ class Test(gammalib.GPythonTestSuite):
343 344 # Return
344 345 return
345 346  
  347 + # Test SIN projection
  348 + def test_skymap_sin(self):
  349 + """
  350 + Test SIN projection.
  351 + """
  352 + # Execute generic test
  353 + self.test_skymap_proj("SIN")
  354 +
  355 + # Return
  356 + return
  357 +
346 358 # Test STG projection
347 359 def test_skymap_stg(self):
348 360 """
... ...