diff --git a/ChangeLog b/ChangeLog index 65291f6..2966fbc 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-07-03 +2018-07-30 * Version 1.6.0 released ======================== + Add ARC WCS projection (#2541) Add GLog::buffer() method. Renamed GPhoton::mcid() to GPhoton::mc_id() Remove at() methods from Python interface (#2515) diff --git a/NEWS b/NEWS index c5c1f78..c8d6eb8 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -5 July 2018 +30 July 2018 1. Introduction @@ -13,6 +13,7 @@ made since the last release of GammaLib. --------------------------------- The following classes have been added: - GModelSpectralExponential +- GWcsARC The following classes have been removed: - none @@ -205,6 +206,10 @@ None 12. Sky module -------------- +The class GWcsARC was added to handle ARC WCS projections (#2541). In addition +the private GWcs::prj_bchk() method was added to implement boundary checking +similar to the checking done in WCS. + A precomputation cache was added to the GWcs::dir2pix() and GWcs::pix2dir() methods so that results are not recomputed for the same arguments. This should provide some speed-up as it avoids unnecessary coordinate transformations. diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 1cbaa4c..890342c 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -34,6 +34,8 @@ Bug fixes Improvements ------------ +* [`2541 `_] - + Add ``ARC`` WCS projection. * Add ``GLog::buffer()`` method * Renamed ``GPhoton::mcid()`` to ``GPhoton::mc_id()`` * [`2513 `_] - diff --git a/include/GWcs.hpp b/include/GWcs.hpp index fdc5564..d472134 100644 --- a/include/GWcs.hpp +++ b/include/GWcs.hpp @@ -146,6 +146,13 @@ protected: // Methods adapted from wcslib::prj.c void prj_ini(void) const; void prj_off(const double& phi0, const double& theta0) const; + int prj_bchk(const double& tol, + const int& nphi, + const int& ntheta, + const int& spt, + double* phi, + double* theta, + int* stat) const; virtual void prj_set(void) const = 0; virtual void prj_x2s(int nx, int ny, int sxy, int spt, const double* x, const double* y, diff --git a/include/GWcsARC.hpp b/include/GWcsARC.hpp new file mode 100644 index 0000000..f7eaa9c --- /dev/null +++ b/include/GWcsARC.hpp @@ -0,0 +1,120 @@ +/*************************************************************************** + * GWcsARC.hpp - Zenithal/azimuthal equidistant (ARC) projection class * + * ----------------------------------------------------------------------- * + * copyright (C) 2018 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GWcsARC.hpp + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition + * @author Juergen Knoedlseder + */ + +#ifndef GWCSARC_HPP +#define GWCSARC_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GWcs.hpp" + + +/***********************************************************************//** + * @class GWcsARC + * + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition + * + * This class implements the Zenithal/azimuthal equidistant projection for + * the World Coordinate System. + ***************************************************************************/ +class GWcsARC : public GWcs { + +public: + // Constructors and destructors + GWcsARC(void); + GWcsARC(const std::string& coords, + const double& crval1, const double& crval2, + const double& crpix1, const double& crpix2, + const double& cdelt1, const double& cdelt2); + GWcsARC(const GWcsARC& wcs); + virtual ~GWcsARC(void); + + // Operators + GWcsARC& operator=(const GWcsARC& wcs); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GWcsARC* clone(void) const; + virtual std::string classname(void) const; + virtual std::string code(void) const; + virtual std::string name(void) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + +private: + // Private methods + void init_members(void); + void copy_members(const GWcsARC& wcs); + void free_members(void); + void prj_set(void) const; + void prj_x2s(int nx, int ny, int sxy, int spt, + const double* x, const double* y, + double* phi, double* theta, int* stat) const; + void prj_s2x(int nphi, int ntheta, int spt, int sxy, + const double* phi, const double* theta, + double* x, double* y, int* stat) const; +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GWcsARC"). + ***************************************************************************/ +inline +std::string GWcsARC::classname(void) const +{ + return ("GWcsARC"); +} + + +/***********************************************************************//** + * @brief Return projection code + * + * @return Projection code. + * + * Returns the projection code "ARC". + ***************************************************************************/ +inline +std::string GWcsARC::code(void) const +{ + return "ARC"; +} + + +/***********************************************************************//** + * @brief Return projection name + * + * @return Projection name. + * + * Returns the projection name. + ***************************************************************************/ +inline +std::string GWcsARC::name(void) const +{ + return "Zenithal/azimuthal equidistant"; +} + +#endif /* GWCSARC_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index a78cfab..a422250 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GammaLib.hpp - Gamma-Ray Astronomy Library Header file * * ----------------------------------------------------------------------- * - * copyright (C) 2006-2017 by Juergen Knoedlseder * + * copyright (C) 2006-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -155,6 +155,7 @@ #include "GWcsRegistry.hpp" #include "GWcs.hpp" #include "GWcsAIT.hpp" +#include "GWcsARC.hpp" #include "GWcsAZP.hpp" #include "GWcsCAR.hpp" #include "GWcsGLS.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index 41b6963..3618f1f 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for include module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2017 Juergen Knoedlseder # +# Copyright (C) 2016-2018 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -117,6 +117,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GWcsRegistry.hpp \ GWcs.hpp \ GWcsAIT.hpp \ + GWcsARC.hpp \ GWcsAZP.hpp \ GWcsCAR.hpp \ GWcsGLS.hpp \ diff --git a/pyext/GWcsARC.i b/pyext/GWcsARC.i new file mode 100644 index 0000000..fb39541 --- /dev/null +++ b/pyext/GWcsARC.i @@ -0,0 +1,64 @@ +/*************************************************************************** + * GWcsARC.i - Zenithal/azimuthal equidistant (ARC) projection class * + * ----------------------------------------------------------------------- * + * copyright (C) 2018 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GWcsARC.i + * @brief Zenithal/azimuthal equidistant (ARC) projection class definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GWcsARC.hpp" +%} + + +/***********************************************************************//** + * @class GWcsARC + * + * @brief Zenithal/azimuthal equidistant (ARC) projection class + ***************************************************************************/ +class GWcsARC : public GWcs { +public: + // Constructors and destructors + GWcsARC(void); + GWcsARC(const std::string& coords, + const double& crval1, const double& crval2, + const double& crpix1, const double& crpix2, + const double& cdelt1, const double& cdelt2); + GWcsARC(const GWcsARC& wcs); + virtual ~GWcsARC(void); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GWcsARC* clone(void) const; + virtual std::string classname(void) const; + virtual std::string code(void) const; + virtual std::string name(void) const; +}; + + +/***********************************************************************//** + * @brief GWcsARC class extension + ***************************************************************************/ +%extend GWcsARC { + GWcsARC copy() { + return (*self); + } +}; diff --git a/pyext/gammalib/sky.i b/pyext/gammalib/sky.i index 746048d..068579e 100644 --- a/pyext/gammalib/sky.i +++ b/pyext/gammalib/sky.i @@ -1,7 +1,7 @@ /*************************************************************************** * sky.i - Sky module * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2017 by Jurgen Knodlseder * + * copyright (C) 2011-2018 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -122,6 +122,7 @@ namespace std { %include "GWcsRegistry.i" %include "GWcs.i" %include "GWcsAIT.i" +%include "GWcsARC.i" %include "GWcsAZP.i" %include "GWcsCAR.i" %include "GWcsGLS.i" diff --git a/src/sky/GWcs.cpp b/src/sky/GWcs.cpp index ff9128b..9d883f0 100644 --- a/src/sky/GWcs.cpp +++ b/src/sky/GWcs.cpp @@ -3090,3 +3090,109 @@ void GWcs::prj_off(const double& phi0, const double& theta0) const // Return return; } + + +/***********************************************************************//** + * @brief Performs bounds checking on native spherical coordinates + * + * @param[in] tol Tolerance for the bounds check [deg] + * @param[in] nphi Number of bins in phi. + * @param[in] ntheta Number of bins in theta. + * @param[in] spt Step size. + * @param[in,out] phi Pointer to phi array [deg] + * @param[in,out] theta Pointer to theta array [deg] + * @param[in,out] stat Pointer to status array. + * + * Performs bounds checking on native spherical coordinates. As returned by + * the deprojection (x2s) routines, native longitude is expected to lie in + * the closed interval [-180,180], with latitude in [-90,90]. + * + * A tolerance may be specified to provide a small allowance for numerical + * imprecision. Values that lie outside the allowed range by not more than + * the specified tolerance will be adjusted back into range. + * + * Code adapted from prj.c::prjbchk(). + ***************************************************************************/ +int GWcs::prj_bchk(const double& tol, + const int& nphi, + const int& ntheta, + const int& spt, + double* phi, + double* theta, + int* stat) const +{ + // Initialize result + int status = 0; + + // Continue only if bounds checking is requested + if (m_bounds) { + + // Initialise array pointers + register double* phip = phi; + register double* thetap = theta; + register int* statp = stat; + + // Loop over all theta bins + for (int itheta = 0; itheta < ntheta; ++itheta) { + + // Loop over all phi bins + for (int iphi = 0; iphi < nphi; ++iphi, phip += spt, thetap += spt, statp++) { + + // Skip values already marked as illegal + if (*statp == 0) { + + // Check lower phi boundary + if (*phip < -180.0) { + if (*phip < -180.0-tol) { + *statp = 1; + status = 1; + } + else { + *phip = -180.0; + } + } + + // Check upper phi boundary + else if (180.0 < *phip) { + if (180.0+tol < *phip) { + *statp = 1; + status = 1; + } + else { + *phip = 180.0; + } + } + + // Check lower theta boundary + if (*thetap < -90.0) { + if (*thetap < -90.0-tol) { + *statp = 1; + status = 1; + } + else { + *thetap = -90.0; + } + } + + // Check upper theta boundary + else if (90.0 < *thetap) { + if (90.0+tol < *thetap) { + *statp = 1; + status = 1; + } + else { + *thetap = 90.0; + } + } + + } // endif: values were marked as valid + + } // endfor: looped over phi + + } // endfor: looped over theta + + } // endif: bounds checking was requested + + // Return status + return status; +} diff --git a/src/sky/GWcsAIT.cpp b/src/sky/GWcsAIT.cpp index fffb7ff..f7c77fc 100644 --- a/src/sky/GWcsAIT.cpp +++ b/src/sky/GWcsAIT.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GWcsAIT.cpp - Aitoff (AIT) projection class * + * GWcsAIT.cpp - Zenithal/azimuthal equidistant (AIT) projection class * * ----------------------------------------------------------------------- * * copyright (C) 2013-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -306,21 +306,21 @@ void GWcsAIT::prj_set(void) const // Initialise projection parameters m_w.assign(4, 0.0); - + // Set undefined parameters if (m_r0 == 0.0) { m_r0 = gammalib::rad2deg; } - + // Precompute m_w[0] = 2.0 * m_r0 * m_r0; m_w[1] = 1.0 / (2.0 * m_w[0]); m_w[2] = m_w[1] / 4.0; m_w[3] = 1.0 / (2.0 * m_r0); - + // Compute fiducial offset prj_off(0.0, 0.0); - + // Return return; } @@ -381,7 +381,7 @@ void GWcsAIT::prj_x2s(int nx, int ny, int sxy, int spt, // Initialise status code and statistics int status = 0; int n_invalid = 0; - + // Do x dependence const double* xp = x; int rowoff = 0; @@ -480,6 +480,11 @@ void GWcsAIT::prj_x2s(int nx, int ny, int sxy, int spt, throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid); } + // Do bounds checking on the native coordinates + // if (prj->bounds&4 && prjbchk(1.0e-13, nx, my, spt, phi, theta, stat)) { + // if (!status) status = PRJERR_BAD_PIX_SET("aitx2s"); + // } + // Return return; } @@ -537,7 +542,7 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy, // Initialise status code and statistics int status = 0; int n_invalid = 0; - + // Do phi dependence const double* phip = phi; int rowoff = 0; @@ -563,7 +568,7 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy, double* yp = y; int* statp = stat; for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) { - + // Compute sin(theta) and cos(theta) double sinthe; double costhe; @@ -578,12 +583,12 @@ void GWcsAIT::prj_s2x(int nphi, int ntheta, int spt, int sxy, } // endfor: phi } // endfor: theta - + // Handle status code if (status == 4) { throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid); } - + // Return return; } diff --git a/src/sky/GWcsARC.cpp b/src/sky/GWcsARC.cpp new file mode 100644 index 0000000..a2097cb --- /dev/null +++ b/src/sky/GWcsARC.cpp @@ -0,0 +1,521 @@ +/*************************************************************************** + * GWcsARC.cpp - Zenithal/azimuthal equidistant (ARC) projection class * + * ----------------------------------------------------------------------- * + * copyright (C) 2018 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GWcsARC.cpp + * @brief Zenithal/azimuthal equidistant (ARC) projection class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GMath.hpp" +#include "GTools.hpp" +#include "GWcsARC.hpp" +#include "GWcsRegistry.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_PRJ_X2S "GWcsARC::prj_x2s(int, int, int, int, double*, double*,"\ + " double*, double*, int*)" +#define G_PRJ_S2X "GWcsARC::prj_s2x(int, int, int, int, double*, double*,"\ + " double*, double*, int*)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + +/* __ Local prototypes ___________________________________________________ */ + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GWcsARC g_wcs_arc_seed; +const GWcsRegistry g_wcs_arc_registry(&g_wcs_arc_seed); + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GWcsARC::GWcsARC(void) : GWcs() +{ + // Initialise class members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Constructor + * + * @param[in] coords Coordinate system. + * @param[in] crval1 X value of reference pixel. + * @param[in] crval2 Y value of reference pixel. + * @param[in] crpix1 X index of reference pixel (starting from 1). + * @param[in] crpix2 Y index of reference pixel (starting from 1). + * @param[in] cdelt1 Increment in x direction at reference pixel [deg]. + * @param[in] cdelt2 Increment in y direction at reference pixel [deg]. + ***************************************************************************/ +GWcsARC::GWcsARC(const std::string& coords, + const double& crval1, const double& crval2, + const double& crpix1, const double& crpix2, + const double& cdelt1, const double& cdelt2) : + GWcs(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2) + +{ + // Initialise class members + init_members(); + + // Setup WCS derived parameters + wcs_set(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] wcs World Coordinate System. + ***************************************************************************/ +GWcsARC::GWcsARC(const GWcsARC& wcs) : GWcs(wcs) +{ + // Initialise class members for clean destruction + init_members(); + + // Copy members + copy_members(wcs); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GWcsARC::~GWcsARC(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] wcs World Coordinate System. + * @return World Coordinate System. + ***************************************************************************/ +GWcsARC& GWcsARC::operator=(const GWcsARC& wcs) +{ + // Execute only if object is not identical + if (this != &wcs) { + + // Copy base class members + this->GWcs::operator=(wcs); + + // Free members + free_members(); + + // Initialise private members for clean destruction + init_members(); + + // Copy members + copy_members(wcs); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + * + * This method properly resets the object to an initial state. + ***************************************************************************/ +void GWcsARC::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GWcs::free_members(); + this->GSkyProjection::free_members(); + + // Initialise members + this->GSkyProjection::init_members(); + this->GWcs::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone instance + * + * @return Pointer to deep copy of Zenithal/azimuthal equidistant projection. + ***************************************************************************/ +GWcsARC* GWcsARC::clone(void) const +{ + return new GWcsARC(*this); +} + + + +/***********************************************************************//** + * @brief Print WCS information + * + * @param[in] chatter Chattiness. + * @return String containing WCS information. + ***************************************************************************/ +std::string GWcsARC::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GWcsARC ==="); + + // Append information + result.append(wcs_print(chatter)); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GWcsARC::init_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] wcs World Coordinate System. + ***************************************************************************/ +void GWcsARC::copy_members(const GWcsARC& wcs) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GWcsARC::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Setup of projection + * + * This method sets up the projection information. The method has been + * adapted from the wcslib function prj.c::arcset. + * + * Given and/or returned: + * m_r0 Reset to 180/pi if 0. + * m_phi0 Reset to 0.0 if undefined. + * m_theta0 Reset to 90.0 if undefined. + * + * Returned: + * m_x0 Fiducial offset in x. + * m_y0 Fiducial offset in y. + * m_w[0] r0*(pi/180) + * m_w[1] (180/pi)/r0 + ***************************************************************************/ +void GWcsARC::prj_set(void) const +{ + // Signal that projection has been set (needs to be done before calling + // the prj_off() method to avoid an endless loop) + m_prjset = true; + + // Initialise projection parameters + m_w.assign(2, 0.0); + + // Set undefined parameters + if (m_r0 == 0.0) { + m_r0 = gammalib::rad2deg; + m_w[0] = 1.0; + m_w[1] = 1.0; + } + else { + m_w[0] = m_r0 * gammalib::deg2rad; + m_w[1] = 1.0/m_w[0]; + } + + // Compute fiducial offset + prj_off(0.0, 0.0); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Pixel-to-spherical deprojection + * + * @param[in] nx X vector length. + * @param[in] ny Y vector length (0=no replication). + * @param[in] sxy Input vector step. + * @param[in] spt Output vector step. + * @param[in] x Vector of projected x coordinates. + * @param[in] y Vector of projected y coordinates. + * @param[out] phi Longitude of the projected point in native spherical + * coordinates [deg]. + * @param[out] theta Latitude of the projected point in native spherical + * coordinates [deg]. + * @param[out] stat Status return value for each vector element (always 0) + * + * @exception GException::wcs_invalid_x_y + * One or more of the (x,y) coordinates were invalid, as indicated + * by the stat vector. + * + * Deproject pixel (x,y) coordinates in the plane of projection to native + * spherical coordinates (phi,theta). + * + * This method has been adapted from the wcslib function prj.c::aitx2s(). + * The interface follows very closely that of wcslib. In contrast to the + * wcslib routine, however, the method assumes that the projection has been + * setup previously (as this will be done by the constructor). + ***************************************************************************/ +void GWcsARC::prj_x2s(int nx, int ny, int sxy, int spt, + const double* x, const double* y, + double* phi, double* theta, int* stat) const +{ + // Initialize projection if required + if (!m_prjset) { + prj_set(); + } + + // Set value replication length mx,my + int mx; + int my; + if (ny > 0) { + mx = nx; + my = ny; + } + else { + mx = 1; + my = 1; + ny = nx; + } + + // Do x dependence + const double* xp = x; + int rowoff = 0; + int rowlen = nx * spt; + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) { + double xj = *xp + m_x0; + double* phip = phi + rowoff; + for (int iy = 0; iy < my; ++iy, phip += rowlen) { + *phip = xj; + } + } + + // Do y dependence + const double* yp = y; + double* phip = phi; + double* thetap = theta; + int* statp = stat; + for (int iy = 0; iy < ny; ++iy, yp += sxy) { + double yj = *yp + m_y0; + double yj2 = yj * yj; + for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) { + double xj = *phip; + double r = std::sqrt(xj*xj + yj2); + if (r == 0.0) { + *phip = 0.0; + *thetap = 90.0; + } + else { + *phip = gammalib::atan2d(xj, -yj); + *thetap = 90.0 - r * m_w[1]; + } + *(statp++) = 0; + } + } + + // Do boundary checking + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat); + + // Handle status code + if (status != 0) { + std::string msg = "One or more of the (x, y) coordinates were invalid " + "for the "+name()+" ("+code()+") projection."; + throw GException::invalid_value(G_PRJ_X2S, msg); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generic spherical-to-pixel projection + * + * @param[in] nphi Longitude vector length. + * @param[in] ntheta Latitude vector length (0=no replication). + * @param[in] spt Input vector step. + * @param[in] sxy Output vector step. + * @param[in] phi Longitude vector of the projected point in native spherical + * coordinates [deg]. + * @param[in] theta Latitude vector of the projected point in native spherical + * coordinates [deg]. + * @param[out] x Vector of projected x coordinates. + * @param[out] y Vector of projected y coordinates. + * @param[out] stat Status return value for each vector element (always 0) + * + * @exception GException::wcs_invalid_phi_theta + * One or more of the (phi,theta) coordinates were invalid, as + * indicated by the stat vector. + * + * Project native spherical coordinates (phi,theta) to pixel (x,y) + * coordinates in the plane of projection. + * + * This method has been adapted from the wcslib function prj.c::aits2x(). + * The interface follows very closely that of wcslib. In contrast to the + * wcslib routine, however, the method assumes that the projection has been + * setup previously (as this will be done by the constructor). + ***************************************************************************/ +void GWcsARC::prj_s2x(int nphi, int ntheta, int spt, int sxy, + const double* phi, const double* theta, + double* x, double* y, int* stat) const +{ + // Initialize projection if required + if (!m_prjset) { + prj_set(); + } + + // Set value replication length mphi,mtheta + int mphi; + int mtheta; + if (ntheta > 0) { + mphi = nphi; + mtheta = ntheta; + } + else { + mphi = 1; + mtheta = 1; + ntheta = nphi; + } + + // Initialise status code and statistics + //int status = 0; + //int n_invalid = 0; + + // Do phi dependence + const double* phip = phi; + int rowoff = 0; + int rowlen = nphi * sxy; + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) { + double sinphi; + double cosphi; + gammalib::sincosd(*phip, &sinphi, &cosphi); + double* xp = x + rowoff; + double* yp = y + rowoff; + for (int itheta = 0; itheta < mtheta; ++itheta) { + *xp = sinphi; + *yp = cosphi; + xp += rowlen; + yp += rowlen; + } + } + + // Do theta dependence + const double* thetap = theta; + double* xp = x; + double* yp = y; + int* statp = stat; + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) { + + // Compute r + double r = m_w[0] * (90.0 - *thetap); + + // Do phi dependence + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) { + *xp = r * (*xp) - m_x0; + *yp = -r * (*yp) - m_y0; + *(statp++) = 0; + } // endfor: phi + + } // endfor: theta + + // Handle status code + /* + if (status == 4) { + throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid); + } + */ + + // Return + return; +} diff --git a/src/sky/GWcsCAR.cpp b/src/sky/GWcsCAR.cpp index 4681d5a..eb36a74 100644 --- a/src/sky/GWcsCAR.cpp +++ b/src/sky/GWcsCAR.cpp @@ -34,6 +34,8 @@ #include "GWcsRegistry.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_PRJ_X2S "GWcsCAR::prj_x2s(int, int, int, int, double*, double*,"\ + " double*, double*, int*)" /* __ Macros _____________________________________________________________ */ @@ -383,6 +385,16 @@ void GWcsCAR::prj_x2s(int nx, int ny, int sxy, int spt, } } + // Do boundary checking + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat); + + // Handle status code + if (status != 0) { + std::string msg = "One or more of the (x, y) coordinates were invalid " + "for the "+name()+" ("+code()+") projection."; + throw GException::invalid_value(G_PRJ_X2S, msg); + } + // Return return; } diff --git a/src/sky/GWcsMER.cpp b/src/sky/GWcsMER.cpp index e6fde11..4bf079a 100644 --- a/src/sky/GWcsMER.cpp +++ b/src/sky/GWcsMER.cpp @@ -34,6 +34,8 @@ #include "GWcsRegistry.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_PRJ_X2S "GWcsMER::prj_x2s(int, int, int, int, double*, double*,"\ + " double*, double*, int*)" /* __ Macros _____________________________________________________________ */ @@ -383,6 +385,16 @@ void GWcsMER::prj_x2s(int nx, int ny, int sxy, int spt, } } + // Do boundary checking + int status = prj_bchk(1.0e-13, nx, my, spt, phi, theta, stat); + + // Handle status code + if (status != 0) { + std::string msg = "One or more of the (x, y) coordinates were invalid " + "for the "+name()+" ("+code()+") projection."; + throw GException::invalid_value(G_PRJ_X2S, msg); + } + // Return return; } diff --git a/src/sky/GWcsMOL.cpp b/src/sky/GWcsMOL.cpp index 8c63ec3..5d1c500 100644 --- a/src/sky/GWcsMOL.cpp +++ b/src/sky/GWcsMOL.cpp @@ -529,6 +529,16 @@ void GWcsMOL::prj_x2s(int nx, int ny, int sxy, int spt, throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid); } + // Do boundary checking + status = prj_bchk(1.0e-11, nx, my, spt, phi, theta, stat); + + // Handle status code + if (status != 0) { + std::string msg = "One or more of the (x, y) coordinates were invalid " + "for the "+name()+" ("+code()+") projection."; + throw GException::invalid_value(G_PRJ_X2S, msg); + } + // Return return; } diff --git a/src/sky/GWcsSFL.cpp b/src/sky/GWcsSFL.cpp index 6b4bee4..e088635 100644 --- a/src/sky/GWcsSFL.cpp +++ b/src/sky/GWcsSFL.cpp @@ -428,6 +428,16 @@ void GWcsSFL::prj_x2s(int nx, int ny, int sxy, int spt, throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid); } + // Do boundary checking + status = prj_bchk(1.0e-12, nx, my, spt, phi, theta, stat); + + // Handle status code + if (status != 0) { + std::string msg = "One or more of the (x, y) coordinates were invalid " + "for the "+name()+" ("+code()+") projection."; + throw GException::invalid_value(G_PRJ_X2S, msg); + } + // Return return; } diff --git a/src/sky/Makefile.am b/src/sky/Makefile.am index 0cc4d25..e3faf2d 100644 --- a/src/sky/Makefile.am +++ b/src/sky/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for sky module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2017 Juergen Knoedlseder # +# Copyright (C) 2016-2018 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -38,6 +38,7 @@ sources = GSkyDir.cpp \ GWcsRegistry.cpp \ GWcs.cpp \ GWcsAIT.cpp \ + GWcsARC.cpp \ GWcsAZP.cpp \ GWcsCAR.cpp \ GWcsGLS.cpp \ diff --git a/test/test_GSky.py b/test/test_GSky.py index 1b5c735..2bf2d14 100755 --- a/test/test_GSky.py +++ b/test/test_GSky.py @@ -358,6 +358,17 @@ class Test(gammalib.GPythonTestSuite): # Return return + # Test ARC projection + def _test_skymap_arc(self): + """ + Test ARC projection. + """ + # Execute generic test + self._test_skymap_proj('ARC') + + # Return + return + # Test AZP projection def _test_skymap_azp(self): """ @@ -565,6 +576,7 @@ class Test(gammalib.GPythonTestSuite): test_support._pickeling(self, gammalib.GSkyRegionMap()) test_support._pickeling(self, gammalib.GSkyRegions()) test_support._pickeling(self, gammalib.GWcsAIT()) + test_support._pickeling(self, gammalib.GWcsARC()) test_support._pickeling(self, gammalib.GWcsAZP()) test_support._pickeling(self, gammalib.GWcsCAR()) test_support._pickeling(self, gammalib.GWcsGLS()) @@ -596,6 +608,7 @@ class Test(gammalib.GPythonTestSuite): test_support._pickeling(self, gammalib.GSkyRegionMap(map)) test_support._pickeling(self, gammalib.GSkyRegions(regions)) test_support._pickeling(self, gammalib.GWcsAIT('GAL',1.0,2.0,1.0,1.0,0.1,0.1)) + test_support._pickeling(self, gammalib.GWcsARC('GAL',1.0,2.0,1.0,1.0,0.1,0.1)) test_support._pickeling(self, gammalib.GWcsAZP('GAL',1.0,2.0,1.0,1.0,0.1,0.1)) test_support._pickeling(self, gammalib.GWcsCAR('GAL',1.0,2.0,1.0,1.0,0.1,0.1)) 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): self.append(self._test_operators, 'Test sky map operators') self.append(self._test_skymap_healpix, 'Test HEALPix map') self.append(self._test_skymap_ait, 'Test AIT projection map') + self.append(self._test_skymap_arc, 'Test ARC projection map') self.append(self._test_skymap_azp, 'Test AZP projection map') self.append(self._test_skymap_car, 'Test CAR projection map') self.append(self._test_skymap_gls, 'Test GLS projection map')