Justine Devin / gammalib

Commit 0fb7d084e3b4466ce62208bf0e7df9e607b8d2eb

Authored by knodlseder
1 parent 55b62897

initial version

Inline Side-by-side
Showing 3 changed files with 753 additions and 0 deletions
include/GWcsAZP.hpp 0 → 100644
  View file @0fb7d08
  1 +/***************************************************************************
  2 + * GWcsAZP.hpp - Zenithal/azimuthal perspective (AZP) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  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 GWcsAZP.hpp
  23 + * @brief Zenithal/azimuthal perspective (AZP) projection class definition
  24 + * @author J. Knodlseder
  25 + */
  26 +
  27 +#ifndef GWCSAZP_HPP
  28 +#define GWCSAZP_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GWcslib.hpp"
  32 +
  33 +
  34 +/***********************************************************************//**
  35 + * @class GWcsAZP
  36 + *
  37 + * @brief Zenithal/azimuthal perspective (AZP) projection class definition
  38 + *
  39 + * This class implements the "zenithal/azimuthal perspective" projection for
  40 + * the World Coordinate System.
  41 + ***************************************************************************/
  42 +class GWcsAZP : public GWcslib {
  43 +
  44 +public:
  45 + // Constructors and destructors
  46 + GWcsAZP(void);
  47 + explicit GWcsAZP(const std::string& coords,
  48 + const double& crval1, const double& crval2,
  49 + const double& crpix1, const double& crpix2,
  50 + const double& cdelt1, const double& cdelt2);
  51 + GWcsAZP(const GWcsAZP& wcs);
  52 + virtual ~GWcsAZP(void);
  53 +
  54 + // Operators
  55 + GWcsAZP& operator= (const GWcsAZP& wcs);
  56 +
  57 + // Implemented pure virtual base class methods
  58 + virtual void clear(void);
  59 + virtual GWcsAZP* clone(void) const;
  60 + virtual std::string code(void) const { return "AZP"; }
  61 + virtual std::string name(void) const { return "zenithal/azimuthal perspective"; }
  62 + virtual std::string print(void) const;
  63 +
  64 +private:
  65 + // Private methods
  66 + void init_members(void);
  67 + void copy_members(const GWcsAZP& wcs);
  68 + void free_members(void);
  69 + void prj_set(void) const;
  70 + void prj_x2s(int nx, int ny, int sxy, int spt,
  71 + const double* x, const double* y,
  72 + double* phi, double* theta, int* stat) const;
  73 + void prj_s2x(int nphi, int ntheta, int spt, int sxy,
  74 + const double* phi, const double* theta,
  75 + double* x, double* y, int* stat) const;
  76 +};
  77 +
  78 +#endif /* GWCSAZP_HPP */
... ...
pyext/GWcsAZP.i 0 → 100644
  View file @0fb7d08
  1 +/***************************************************************************
  2 + * GWcsAZP.i - Zenithal/azimuthal perspective (AZP) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  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 GWcsAZP.i
  23 + * @brief Zenithal/azimuthal perspective (AZP) projection class Python interface
  24 + * @author J. Knodlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GWcsAZP.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GWcsAZP
  34 + *
  35 + * @brief Zenithal/azimuthal perspective (AZP) projection class definition
  36 + ***************************************************************************/
  37 +class GWcsAZP : public GWcslib {
  38 +public:
  39 + // Constructors and destructors
  40 + GWcsAZP(void);
  41 + explicit GWcsAZP(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 + GWcsAZP(const GWcsAZP& wcs);
  46 + virtual ~GWcsAZP(void);
  47 +
  48 + // Implemented pure virtual methods
  49 + virtual void clear(void);
  50 + virtual GWcsAZP* clone(void) const;
  51 + virtual std::string code(void) const;
  52 + virtual std::string name(void) const;
  53 +};
  54 +
  55 +
  56 +/***********************************************************************//**
  57 + * @brief GWcsAZP class extension
  58 + ***************************************************************************/
  59 +%extend GWcsAZP {
  60 + GWcsAZP copy() {
  61 + return (*self);
  62 + }
  63 +};
... ...
src/sky/GWcsAZP.cpp 0 → 100644
  View file @0fb7d08
  1 +/***************************************************************************
  2 + * GWcsAZP.cpp - Zenithal/azimuthal perspective (AZP) projection class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  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 GWcsAZP.cpp
  23 + * @brief Zenithal/azimuthal perspective (AZP) projection class implementation
  24 + * @author J. Knodlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GWcsAZP.hpp"
  34 +#include "GWcsRegistry.hpp"
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_PRJ_SET "GWcsAZP::prj_set()"
  38 +#define G_PRJ_X2S "GWcsAZP::prj_x2s(int,int,int,int,double*,double*,double*,"\
  39 + "double*,int*)"
  40 +#define G_PRJ_S2X "GWcsAZP::prj_s2x(int,int,int,int,double*,double*,double*,"\
  41 + "double*,int*)"
  42 +
  43 +/* __ Macros _____________________________________________________________ */
  44 +
  45 +/* __ Coding definitions _________________________________________________ */
  46 +
  47 +/* __ Debug definitions __________________________________________________ */
  48 +
  49 +/* __ Local prototypes ___________________________________________________ */
  50 +
  51 +/* __ Constants __________________________________________________________ */
  52 +
  53 +/* __ Globals ____________________________________________________________ */
  54 +const GWcsAZP g_wcs_azp_seed;
  55 +const GWcsRegistry g_wcs_azp_registry(&g_wcs_azp_seed);
  56 +
  57 +
  58 +/*==========================================================================
  59 + = =
  60 + = Constructors/destructors =
  61 + = =
  62 + ==========================================================================*/
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Void constructor
  66 + ***************************************************************************/
  67 +GWcsAZP::GWcsAZP(void) : GWcslib()
  68 +{
  69 + // Initialise class members
  70 + init_members();
  71 +
  72 + // Return
  73 + return;
  74 +}
  75 +
  76 +
  77 +/***********************************************************************//**
  78 + * @brief Constructor
  79 + *
  80 + * @param[in] coords Coordinate system.
  81 + * @param[in] crval1 X value of reference pixel.
  82 + * @param[in] crval2 Y value of reference pixel.
  83 + * @param[in] crpix1 X index of reference pixel (starting from 1).
  84 + * @param[in] crpix2 Y index of reference pixel (starting from 1).
  85 + * @param[in] cdelt1 Increment in x direction at reference pixel [deg].
  86 + * @param[in] cdelt2 Increment in y direction at reference pixel [deg].
  87 + ***************************************************************************/
  88 +GWcsAZP::GWcsAZP(const std::string& coords,
  89 + const double& crval1, const double& crval2,
  90 + const double& crpix1, const double& crpix2,
  91 + const double& cdelt1, const double& cdelt2) :
  92 + GWcslib(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2)
  93 +
  94 +{
  95 + // Initialise class members
  96 + init_members();
  97 +
  98 + // Return
  99 + return;
  100 +}
  101 +
  102 +
  103 +/***********************************************************************//**
  104 + * @brief Copy constructor
  105 + *
  106 + * @param[in] wcs World Coordinate System.
  107 + ***************************************************************************/
  108 +GWcsAZP::GWcsAZP(const GWcsAZP& wcs) : GWcslib(wcs)
  109 +{
  110 + // Initialise class members for clean destruction
  111 + init_members();
  112 +
  113 + // Copy members
  114 + copy_members(wcs);
  115 +
  116 + // Return
  117 + return;
  118 +}
  119 +
  120 +
  121 +/***********************************************************************//**
  122 + * @brief Destructor
  123 + ***************************************************************************/
  124 +GWcsAZP::~GWcsAZP(void)
  125 +{
  126 + // Free members
  127 + free_members();
  128 +
  129 + // Return
  130 + return;
  131 +}
  132 +
  133 +
  134 +/*==========================================================================
  135 + = =
  136 + = Operators =
  137 + = =
  138 + ==========================================================================*/
  139 +
  140 +/***********************************************************************//**
  141 + * @brief Assignment operator
  142 + *
  143 + * @param[in] wcs World Coordinate System.
  144 + ***************************************************************************/
  145 +GWcsAZP& GWcsAZP::operator= (const GWcsAZP& wcs)
  146 +{
  147 + // Execute only if object is not identical
  148 + if (this != &wcs) {
  149 +
  150 + // Copy base class members
  151 + this->GWcslib::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 instance
  177 + *
  178 + * This method properly resets the object to an initial state.
  179 + ***************************************************************************/
  180 +void GWcsAZP::clear(void)
  181 +{
  182 + // Free class members (base and derived classes, derived class first)
  183 + free_members();
  184 + this->GWcslib::free_members();
  185 + this->GWcs::free_members();
  186 +
  187 + // Initialise members
  188 + this->GWcs::init_members();
  189 + this->GWcslib::init_members();
  190 + init_members();
  191 +
  192 + // Return
  193 + return;
  194 +}
  195 +
  196 +
  197 +/***********************************************************************//**
  198 + * @brief Clone instance
  199 + ***************************************************************************/
  200 +GWcsAZP* GWcsAZP::clone(void) const
  201 +{
  202 + return new GWcsAZP(*this);
  203 +}
  204 +
  205 +
  206 +
  207 +/***********************************************************************//**
  208 + * @brief Print WCS information
  209 + ***************************************************************************/
  210 +std::string GWcsAZP::print(void) const
  211 +{
  212 + // Initialise result string
  213 + std::string result;
  214 +
  215 + // Append header
  216 + result.append("=== GWcsAZP ===\n");
  217 + result.append(wcs_print());
  218 +
  219 + // Return result
  220 + return result;
  221 +}
  222 +
  223 +
  224 +/*==========================================================================
  225 + = =
  226 + = Private methods =
  227 + = =
  228 + ==========================================================================*/
  229 +
  230 +/***********************************************************************//**
  231 + * @brief Initialise class members
  232 + *
  233 + * This method sets up the World Coordinate System by calling wcs_set().
  234 + ***************************************************************************/
  235 +void GWcsAZP::init_members(void)
  236 +{
  237 + // Setup World Coordinate System
  238 + wcs_set();
  239 +
  240 + // Return
  241 + return;
  242 +}
  243 +
  244 +
  245 +/***********************************************************************//**
  246 + * @brief Copy class members
  247 + *
  248 + * @param[in] wcs World Coordinate System.
  249 + ***************************************************************************/
  250 +void GWcsAZP::copy_members(const GWcsAZP& wcs)
  251 +{
  252 + // Return
  253 + return;
  254 +}
  255 +
  256 +
  257 +/***********************************************************************//**
  258 + * @brief Delete class members
  259 + ***************************************************************************/
  260 +void GWcsAZP::free_members(void)
  261 +{
  262 + // Return
  263 + return;
  264 +}
  265 +
  266 +
  267 +/***********************************************************************//**
  268 + * @brief Setup of projection
  269 + *
  270 + * This method sets up the projection information. The method has been
  271 + * adapted from the wcslib function prj.c::carset.
  272 + *
  273 + * @exception GException::wcs_invalid_parameter
  274 + * PV(1) or PV(2) are invalid.
  275 + *
  276 + * Given:
  277 + * m_pv[1] Distance parameter, mu in units of r0.
  278 + * m_pv[2] Tilt angle, gamma in degrees.
  279 + *
  280 + * Given and/or returned:
  281 + * m_r0 Reset to 180/pi if 0.
  282 + * m_phi0 Reset to 0.0 if undefined.
  283 + * m_theta0 Reset to 90.0 if undefined.
  284 + *
  285 + * Returned:
  286 + * m_x0 Fiducial offset in x.
  287 + * m_y0 Fiducial offset in y.
  288 + * m_w[0] r0*(mu+1)
  289 + * m_w[1] tan(gamma)
  290 + * m_w[2] sec(gamma)
  291 + * m_w[3] cos(gamma)
  292 + * m_w[4] sin(gamma)
  293 + * m_w[5] asin(-1/mu) for |mu| >= 1, -90 otherwise
  294 + * m_w[6] mu*cos(gamma)
  295 + * m_w[7] 1 if |mu*cos(gamma)| < 1, 0 otherwise
  296 + ***************************************************************************/
  297 +void GWcsAZP::prj_set(void) const
  298 +{
  299 + // Initialise projection parameters
  300 + m_w.assign(7,0.0);
  301 +
  302 + // Set undefined parameters
  303 + if (undefined(m_pv[1])) m_pv[1] = 0.0;
  304 + if (undefined(m_pv[2])) m_pv[2] = 0.0;
  305 + if (m_r0 == 0.0) m_r0 = rad2deg;
  306 +
  307 + // Precompute
  308 + m_w[0] = m_r0*(m_pv[1] + 1.0);
  309 + if (m_w[0] == 0.0) {
  310 + std::string message = "PV(1)=-1 encountered.";
  311 + throw GException::wcs_invalid_parameter(G_PRJ_SET, message);
  312 + }
  313 + m_w[3] = cosd(m_pv[2]);
  314 + if (m_w[3] == 0.0) {
  315 + std::string message = "cos(PV(2))=0 encountered.";
  316 + throw GException::wcs_invalid_parameter(G_PRJ_SET, message);
  317 + }
  318 + m_w[2] = 1.0/m_w[3];
  319 + m_w[4] = sind(m_pv[2]);
  320 + m_w[1] = m_w[4] / m_w[3];
  321 + if (std::abs(m_pv[1]) > 1.0)
  322 + m_w[5] = asind(-1.0/m_pv[1]);
  323 + else
  324 + m_w[5] = -90.0;
  325 + m_w[6] = m_pv[1] * m_w[3];
  326 + m_w[7] = (std::abs(m_w[6]) < 1.0) ? 1.0 : 0.0;
  327 +
  328 + // Compute fiducial offset
  329 + prj_off(0.0, 90.0);
  330 +
  331 + // Signal that projection has been set
  332 + m_prjset = true;
  333 +
  334 + // Return
  335 + return;
  336 +}
  337 +
  338 +
  339 +/***********************************************************************//**
  340 + * @brief Cartesian-to-spherical deprojection
  341 + *
  342 + * @param[in] nx X vector length.
  343 + * @param[in] ny Y vector length (0=no replication).
  344 + * @param[in] sxy Input vector step.
  345 + * @param[in] spt Output vector step.
  346 + * @param[in] x Vector of projected x coordinates.
  347 + * @param[in] y Vector of projected y coordinates.
  348 + * @param[out] phi Longitude of the projected point in native spherical
  349 + * coordinates [deg].
  350 + * @param[out] theta Latitude of the projected point in native spherical
  351 + * coordinates [deg].
  352 + * @param[out] stat Status return value for each vector element (always 0)
  353 + *
  354 + * @exception GException::wcs_invalid_x_y
  355 + * One or more of the (x,y) coordinates were invalid, as indicated
  356 + * by the stat vector.
  357 + *
  358 + * Deproject Cartesian (x,y) coordinates in the plane of projection to native
  359 + * spherical coordinates (phi,theta).
  360 + *
  361 + * This method has been adapted from the wcslib function prj.c::carx2s().
  362 + * The interface follows very closely that of wcslib. In contrast to the
  363 + * wcslib routine, however, the method assumes that the projection has been
  364 + * setup previsouly (as this will be done by the constructor).
  365 + ***************************************************************************/
  366 +void GWcsAZP::prj_x2s(int nx, int ny, int sxy, int spt,
  367 + const double* x, const double* y,
  368 + double* phi, double* theta, int* stat) const
  369 +{
  370 + // Set tolerance
  371 + const double tol = 1.0e-13;
  372 +
  373 + // Initialize projection if required
  374 + if (!m_prjset)
  375 + prj_set();
  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 + // Initialise status code and statistics
  391 + int status = 0;
  392 + int n_invalid = 0;
  393 +
  394 + // Do x dependence
  395 + const double* xp = x;
  396 + int rowoff = 0;
  397 + int rowlen = nx * spt;
  398 + for (int ix = 0; ix < nx; ++ix, rowoff += spt, xp += sxy) {
  399 + double xj = *xp + m_x0;
  400 + double* phip = phi + rowoff;
  401 + for (int iy = 0; iy < my; ++iy, phip += rowlen)
  402 + *phip = xj;
  403 + }
  404 +
  405 + // Do y dependence
  406 + const double* yp = y;
  407 + double* phip = phi;
  408 + double* thetap = theta;
  409 + int* statp = stat;
  410 + for (int iy = 0; iy < ny; ++iy, yp += sxy) {
  411 + double yj = *yp + m_y0;
  412 + double yc = yj * m_w[3];
  413 + double yc2 = yc*yc;
  414 + double q = m_w[0] + yj*m_w[4];
  415 + for (int ix = 0; ix < mx; ++ix, phip += spt, thetap += spt) {
  416 + double xj = *phip;
  417 + double r = sqrt(xj*xj + yc2);
  418 + if (r == 0.0) {
  419 + *phip = 0.0;
  420 + *thetap = 90.0;
  421 + *(statp++) = 0;
  422 +
  423 + }
  424 + else {
  425 + *phip = atan2d(xj, -yc);
  426 + double s = r / q;
  427 + double t = s * m_pv[1]/sqrt(s*s + 1.0);
  428 + s = atan2d(1.0, s);
  429 + if (std::abs(t) > 1.0) {
  430 + if (std::abs(t) > 1.0+tol) {
  431 + *thetap = 0.0;
  432 + *(statp++) = 1;
  433 + status = 3;
  434 + n_invalid++;
  435 + continue;
  436 + }
  437 + t = (t < 0) ? -90.0 : 90.0;
  438 + }
  439 + else
  440 + t = asind(t);
  441 + double a = s - t;
  442 + double b = s + t + 180.0;
  443 + if (a > 90.0) a -= 360.0;
  444 + if (b > 90.0) b -= 360.0;
  445 + *thetap = (a > b) ? a : b;
  446 + *(statp++) = 0;
  447 + }
  448 + }
  449 + }
  450 +
  451 + // Handle status code
  452 + if (status == 3)
  453 + throw GException::wcs_invalid_x_y(G_PRJ_X2S, n_invalid);
  454 +
  455 + // Return
  456 + return;
  457 +}
  458 +
  459 +
  460 +/***********************************************************************//**
  461 + * @brief Generic spherical-to-Cartesian projection
  462 + *
  463 + * @param[in] nphi Longitude vector length.
  464 + * @param[in] ntheta Latitude vector length (0=no replication).
  465 + * @param[in] spt Input vector step.
  466 + * @param[in] sxy Output vector step.
  467 + * @param[in] phi Longitude vector of the projected point in native spherical
  468 + * coordinates [deg].
  469 + * @param[in] theta Latitude vector of the projected point in native spherical
  470 + * coordinates [deg].
  471 + * @param[out] x Vector of projected x coordinates.
  472 + * @param[out] y Vector of projected y coordinates.
  473 + * @param[out] stat Status return value for each vector element (always 0)
  474 + *
  475 + * @exception GException::wcs_invalid_phi_theta
  476 + * One or more of the (phi,theta) coordinates were invalid, as
  477 + * indicated by the stat vector.
  478 + *
  479 + * Project native spherical coordinates (phi,theta) to Cartesian (x,y)
  480 + * coordinates in the plane of projection.
  481 + *
  482 + * This method has been adapted from the wcslib function prj.c::cars2x().
  483 + * The interface follows very closely that of wcslib. In contrast to the
  484 + * wcslib routine, however, the method assumes that the projection has been
  485 + * setup previsouly (as this will be done by the constructor).
  486 + ***************************************************************************/
  487 +void GWcsAZP::prj_s2x(int nphi, int ntheta, int spt, int sxy,
  488 + const double* phi, const double* theta,
  489 + double* x, double* y, int* stat) const
  490 +{
  491 + // Initialize projection if required
  492 + if (!m_prjset)
  493 + prj_set();
  494 +
  495 + // Set value replication length mphi,mtheta
  496 + int mphi;
  497 + int mtheta;
  498 + if (ntheta > 0) {
  499 + mphi = nphi;
  500 + mtheta = ntheta;
  501 + }
  502 + else {
  503 + mphi = 1;
  504 + mtheta = 1;
  505 + ntheta = nphi;
  506 + }
  507 +
  508 + // Initialise status code and statistics
  509 + int status = 0;
  510 + int n_invalid = 0;
  511 +
  512 + // Do phi dependence
  513 + const double* phip = phi;
  514 + int rowoff = 0;
  515 + int rowlen = nphi * sxy;
  516 + for (int iphi = 0; iphi < nphi; ++iphi, rowoff += sxy, phip += spt) {
  517 + double sinphi;
  518 + double cosphi;
  519 + sincosd(*phip, &sinphi, &cosphi);
  520 + double* xp = x + rowoff;
  521 + double* yp = y + rowoff;
  522 + for (int itheta = 0; itheta < mtheta; ++itheta) {
  523 + *xp = sinphi;
  524 + *yp = cosphi;
  525 + xp += rowlen;
  526 + yp += rowlen;
  527 + }
  528 + }
  529 +
  530 + // Do theta dependence
  531 + const double* thetap = theta;
  532 + double* xp = x;
  533 + double* yp = y;
  534 + int* statp = stat;
  535 + for (int itheta = 0; itheta < ntheta; ++itheta, thetap += spt) {
  536 +
  537 + // Compute sin(theta) and cos(theta)
  538 + double sinthe;
  539 + double costhe;
  540 + sincosd(*thetap, &sinthe, &costhe);
  541 +
  542 + // ...
  543 + for (int iphi = 0; iphi < mphi; ++iphi, xp += sxy, yp += sxy) {
  544 +
  545 + double s = m_w[1]*(*yp);
  546 + double t = (m_pv[1] + sinthe) + costhe*s;
  547 +
  548 + // ...
  549 + if (t == 0.0) {
  550 + *xp = 0.0;
  551 + *yp = 0.0;
  552 + *(statp++) = 1;
  553 + status = 4;
  554 + n_invalid++;
  555 + }
  556 +
  557 + // ...
  558 + else {
  559 + double r = m_w[0]*costhe/t;
  560 +
  561 + // Bounds checking
  562 + int istat = 0;
  563 + if (m_bounds) {
  564 +
  565 + // Is there overlap?
  566 + if (*thetap < m_w[5]) {
  567 + istat = 1;
  568 + status = 4;
  569 + n_invalid++;
  570 + }
  571 +
  572 + // Is there divergence?
  573 + else if (m_w[7] > 0.0) {
  574 + t = m_pv[1] / sqrt(1.0 + s*s);
  575 + if (std::abs(t) <= 1.0) {
  576 + s = atand(-s);
  577 + t = asind(t);
  578 + double a = s - t;
  579 + double b = s + t + 180.0;
  580 + if (a > 90.0) a -= 360.0;
  581 + if (b > 90.0) b -= 360.0;
  582 + if (*thetap < ((a > b) ? a : b)) {
  583 + istat = 1;
  584 + status = 4;
  585 + n_invalid++;
  586 + }
  587 + }
  588 + }
  589 + }
  590 +
  591 + // Set value
  592 + *xp = r*(*xp) - m_x0;
  593 + *yp = -r*(*yp) * m_w[2] - m_y0;
  594 + *(statp++) = istat;
  595 + } // endelse
  596 + } // endfor: phi
  597 + } // endfor: theta
  598 +
  599 + // Handle status code
  600 + if (status == 4)
  601 + throw GException::wcs_invalid_phi_theta(G_PRJ_S2X, n_invalid);
  602 +
  603 + // Return
  604 + return;
  605 +}
  606 +
  607 +
  608 +/*==========================================================================
  609 + = =
  610 + = Friends =
  611 + = =
  612 + ==========================================================================*/
... ...