diff --git a/include/GHealpix.hpp b/include/GHealpix.hpp index 1dd2920..ec2504c 100644 --- a/include/GHealpix.hpp +++ b/include/GHealpix.hpp @@ -61,16 +61,14 @@ public: // Implemented pure virtual base class methods virtual void clear(void); virtual GHealpix* clone(void) const; + virtual int size(void) const; virtual std::string code(void) const; virtual std::string name(void) const; virtual void read(const GFitsHDU* hdu); virtual void write(GFitsHDU* hdu) const; - virtual double omega(const int& pix) const; - virtual double omega(const GSkyPixel& pix) const; - virtual GSkyDir pix2dir(const int& pix) const; - virtual int dir2pix(const GSkyDir& dir) const; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const; + virtual double omega(const GSkyPixel& pixel) const; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const; virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods @@ -108,6 +106,20 @@ private: /***********************************************************************//** + * @brief Return dimension of projection + * + * @return Dimension of projection. + * + * Returns the dimension of the projection. + ***************************************************************************/ +inline +int GHealpix::size(void) const +{ + return 1; +} + + +/***********************************************************************//** * @brief Return projection code * * @return Projection code. diff --git a/include/GSkyProjection.hpp b/include/GSkyProjection.hpp index 783bfd4..aac4f56 100644 --- a/include/GSkyProjection.hpp +++ b/include/GSkyProjection.hpp @@ -62,16 +62,14 @@ public: // Pure virtual methods virtual void clear(void) = 0; virtual GSkyProjection* clone(void) const = 0; + virtual int size(void) const = 0; virtual std::string code(void) const = 0; virtual std::string name(void) const = 0; virtual void read(const GFitsHDU* hdu) = 0; virtual void write(GFitsHDU* hdu) const = 0; - virtual double omega(const int& pix) const = 0; - virtual double omega(const GSkyPixel& pix) const = 0; - virtual GSkyDir pix2dir(const int& pix) const = 0; - virtual int dir2pix(const GSkyDir& dir) const = 0; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const = 0; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const = 0; + virtual double omega(const GSkyPixel& pixel) const = 0; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const = 0; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const = 0; virtual std::string print(const GChatter& chatter = NORMAL) const = 0; // Virtual methods diff --git a/include/GSkymap.hpp b/include/GSkymap.hpp index b88e225..8deca76 100644 --- a/include/GSkymap.hpp +++ b/include/GSkymap.hpp @@ -45,20 +45,22 @@ /***********************************************************************//** * @class GSkymap * - * @brief GSkymap class interface defintion + * @brief Sky map class * - * This class implements a set of skymaps. Skymaps may either be present in + * This class implements a set of sky maps. Skymaps may either be present in * the HEALPix format or in any (supported) WCS format. Skymap pixels are * stored in a double precision array indexed by (x,y,map), with the x * axis being the most rapidely varying axis. * - * Skymap pixels may be either accessed via 1D operators (i,map), where i - * runs over the (x,y) direction of the map, or via 2D operators (pixel,map), + * Skymap pixels may be either accessed via index operators (i,map), where i + * runs over the (x,y) direction of the map, or via pixel operators (pixel,map), * where pixel implements a 2D (x,y) pixel direction. The first operator is * the preferred access method for HEALPix maps while the second operator is * the preferred access method for WCS maps. Conversion methods between the * index or sky pixel and the true physical sky direction are provided by the * pix2dir() and dir2pix() methods. + * + * @todo Rewrite class description. ***************************************************************************/ class GSkymap : public GBase { @@ -83,44 +85,38 @@ public: virtual ~GSkymap(void); // Operators - GSkymap& operator=(const GSkymap& map); - - // 1D pixel methods - double& operator()(const int& pixel, const int& map = 0); - const double& operator()(const int& pixel, const int& map = 0) const; - GSkyDir pix2dir(const int& pix) const; - int dir2pix(const GSkyDir& dir) const; - double omega(const int& pix) const; - - // 2D pixel methods + GSkymap& operator=(const GSkymap& map); + double& operator()(const int& index, const int& map = 0); + const double& operator()(const int& index, const int& map = 0) const; double& operator()(const GSkyPixel& pixel, const int& map = 0); const double& operator()(const GSkyPixel& pixel, const int& map = 0) const; - GSkyDir xy2dir(const GSkyPixel& pix) const; - GSkyPixel dir2xy(const GSkyDir& dir) const; - double omega(const GSkyPixel& pix) const; - - // Sky direction methods double operator()(const GSkyDir& dir, const int& map = 0) const; // Methods - void clear(void); - GSkymap* clone(void) const; - void load(const std::string& filename); - void save(const std::string& filename, bool clobber = false) const; - void read(const GFitsHDU* hdu); - void write(GFits* file) const; - int npix(void) const; - int nx(void) const; - int ny(void) const; - int nmaps(void) const; - int xy2pix(const GSkyPixel& pix) const; - GSkyPixel pix2xy(const int& pix) const; - GSkyProjection* projection(void) const { return m_proj; } - void projection(const GSkyProjection& proj); - double* pixels(void) const { return m_pixels; } - bool contains(const GSkyDir& dir) const; - bool contains(const GSkyPixel& pixel) const; - std::string print(const GChatter& chatter = NORMAL) const; + void clear(void); + GSkymap* clone(void) const; + const int& npix(void) const; + const int& nx(void) const; + const int& ny(void) const; + const int& nmaps(void) const; + GSkyPixel inx2pix(const int& index) const; + GSkyDir inx2dir(const int& index) const; + GSkyDir pix2dir(const GSkyPixel& pixel) const; + int pix2inx(const GSkyPixel& pixel) const; + int dir2inx(const GSkyDir& dir) const; + GSkyPixel dir2pix(const GSkyDir& dir) const; + double omega(const int& index) const; + double omega(const GSkyPixel& pixel) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyPixel& pixel) const; + const GSkyProjection* projection(void) const; + void projection(const GSkyProjection& proj); + const double* pixels(void) const; + void load(const std::string& filename); + void save(const std::string& filename, bool clobber = false) const; + void read(const GFitsHDU* hdu); + void write(GFits* file) const; + std::string print(const GChatter& chatter = NORMAL) const; private: // Private methods @@ -148,4 +144,86 @@ private: double* m_pixels; //!< Pointer to skymap pixels }; + +/***********************************************************************//** + * @brief Returns number of pixels + * + * @return Number of pixels in one sky map. + * + * Returns the number of pixels in one sky map. + ***************************************************************************/ +inline +const int& GSkymap::npix(void) const +{ + return m_num_pixels; +} + + +/***********************************************************************//** + * @brief Returns number of pixels in x coordinate + * + * @return Number of pixels in the X direction. + * + * Returns the number of pixels in the X direction. If the sky map is a + * one dimensional array (which is the case for the Healpix projection), + * the method returns 0. + ***************************************************************************/ +inline +const int& GSkymap::nx(void) const +{ + return m_num_x; +} + + +/***********************************************************************//** + * @brief Returns number of pixels in y coordinate + * + * @return Number of pixels in the Y direction. + * + * Returns the number of pixels in the Y direction. If the sky map is a + * one dimensional array (which is the case for the Healpix projection), + * the method returns 0. + ***************************************************************************/ +inline +const int& GSkymap::ny(void) const +{ + return m_num_y; +} + + +/***********************************************************************//** + * @brief Returns number of maps + * + * @return Number of maps in the sky map object. + ***************************************************************************/ +inline +const int& GSkymap::nmaps(void) const +{ + return m_num_maps; +} + + +/***********************************************************************//** + * @brief Returns pointer to sky projection + * + * @return Pointer to sky projection (NULL if no projection is defined). + ***************************************************************************/ +inline +const GSkyProjection* GSkymap::projection(void) const +{ + return m_proj; +} + + +/***********************************************************************//** + * @brief Returns pointer to pixel data + * + * @return Pointer to pixel data. + ***************************************************************************/ +inline +const double* GSkymap::pixels(void) const +{ + return m_pixels; +} + #endif /* GSKYMAP_HPP */ diff --git a/include/GWcs.hpp b/include/GWcs.hpp index cf0f90d..399c151 100644 --- a/include/GWcs.hpp +++ b/include/GWcs.hpp @@ -67,14 +67,12 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const = 0; // Implemented virtual base class methods + virtual int size(void) const; virtual void read(const GFitsHDU* hdu); virtual void write(GFitsHDU* hdu) const; - virtual double omega(const int& pix) const; - virtual double omega(const GSkyPixel& pix) const; - virtual GSkyDir pix2dir(const int& pix) const; - virtual int dir2pix(const GSkyDir& dir) const; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const; + virtual double omega(const GSkyPixel& pixel) const; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const; // Other methods void set(const std::string& coords, @@ -208,4 +206,18 @@ protected: //struct spcprm spc }; + +/***********************************************************************//** + * @brief Return dimension of projection + * + * @return Dimension of projection. + * + * Returns the dimension of the projection. + ***************************************************************************/ +inline +int GWcs::size(void) const +{ + return 2; +} + #endif /* GWCS_HPP */ diff --git a/inst/com/src/GCOMEventCube.cpp b/inst/com/src/GCOMEventCube.cpp index 2a26f10..92c9f5b 100644 --- a/inst/com/src/GCOMEventCube.cpp +++ b/inst/com/src/GCOMEventCube.cpp @@ -504,7 +504,7 @@ int GCOMEventCube::number(void) const double number = 0.0; // Get pointer on skymap pixels - double* pixels = m_map.pixels(); + const double* pixels = m_map.pixels(); // Sum event cube if (size() > 0 && pixels != NULL) { @@ -708,7 +708,7 @@ void GCOMEventCube::set_scatter_directions(void) for (int iy = 0; iy < npsi(); ++iy) { for (int ix = 0; ix < nchi(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - m_dirs.push_back(m_map.xy2dir(pixel)); + m_dirs.push_back(m_map.pix2dir(pixel)); m_omega.push_back(m_map.omega(pixel)); } } @@ -871,7 +871,7 @@ void GCOMEventCube::set_bin(const int& index) m_dir.phibar(m_phi[iphi]); // Set pointers - m_bin.m_counts = &(m_map.pixels()[index]); + m_bin.m_counts = const_cast(&(m_map.pixels()[index])); m_bin.m_omega = &(m_omega[ipix]); // Return diff --git a/inst/com/src/GCOMSupport.cpp b/inst/com/src/GCOMSupport.cpp index bc5e522..66ab2f0 100644 --- a/inst/com/src/GCOMSupport.cpp +++ b/inst/com/src/GCOMSupport.cpp @@ -54,7 +54,7 @@ void com_wcs_mer2car(GSkymap& map) { // Get WCS poiunter - GWcs* wcs = dynamic_cast(map.projection()); + const GWcs* wcs = dynamic_cast(map.projection()); // Apply kluge only if the skymap is in Mercator's projection if (wcs != NULL && wcs->code() == "MER") { diff --git a/inst/cta/src/GCTAEventCube.cpp b/inst/cta/src/GCTAEventCube.cpp index 056c984..66f9184 100644 --- a/inst/cta/src/GCTAEventCube.cpp +++ b/inst/cta/src/GCTAEventCube.cpp @@ -418,7 +418,7 @@ int GCTAEventCube::number(void) const double number = 0.0; // Get pointer on skymap pixels - double* pixels = m_map.pixels(); + const double* pixels = m_map.pixels(); // Sum event cube if (size() > 0 && pixels != NULL) { @@ -670,7 +670,7 @@ void GCTAEventCube::set_directions(void) for (int iy = 0; iy < ny(); ++iy) { for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - m_dirs.push_back(GCTAInstDir(m_map.xy2dir(pixel))); + m_dirs.push_back(GCTAInstDir(m_map.pix2dir(pixel))); m_omega.push_back(m_map.omega(pixel)); } } @@ -770,24 +770,27 @@ void GCTAEventCube::set_bin(const int& index) { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) + if (index < 0 || index >= size()) { throw GException::out_of_range(G_SET_BIN, index, 0, size()-1); + } #endif // Check for the existence of energies and energy widths - if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) + if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) { throw GCTAException::no_energies(G_SET_BIN); + } // Check for the existence of sky directions and solid angles - if (m_dirs.size() != npix() || m_omega.size() != npix()) + if (m_dirs.size() != npix() || m_omega.size() != npix()) { throw GCTAException::no_dirs(G_SET_BIN); + } // Get pixel and energy bin indices. int ipix = index % npix(); int ieng = index / npix(); // Set pointers - m_bin.m_counts = &(m_map.pixels()[index]); + m_bin.m_counts = const_cast(&(m_map.pixels()[index])); m_bin.m_energy = &(m_energies[ieng]); m_bin.m_time = &m_time; m_bin.m_dir = &(m_dirs[ipix]); diff --git a/inst/lat/src/GLATEventCube.cpp b/inst/lat/src/GLATEventCube.cpp index c128380..a5612cd 100644 --- a/inst/lat/src/GLATEventCube.cpp +++ b/inst/lat/src/GLATEventCube.cpp @@ -388,7 +388,7 @@ int GLATEventCube::number(void) const double number = 0.0; // Get pointer on skymap pixels - double* pixels = m_map.pixels(); + const double* pixels = m_map.pixels(); // Sum event cube if (size() > 0 && pixels != NULL) { @@ -561,7 +561,7 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const int iy = 0; for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - double distance = m_map.xy2dir(pixel).dist_deg(srcDir); + double distance = m_map.pix2dir(pixel).dist_deg(srcDir); if (distance < radius) { radius = distance; } @@ -571,7 +571,7 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const iy = ny()-1; for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - double distance = m_map.xy2dir(pixel).dist_deg(srcDir); + double distance = m_map.pix2dir(pixel).dist_deg(srcDir); if (distance < radius) { radius = distance; } @@ -581,7 +581,7 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const int ix = 0; for (int iy = 0; iy < ny(); ++iy) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - double distance = m_map.xy2dir(pixel).dist_deg(srcDir); + double distance = m_map.pix2dir(pixel).dist_deg(srcDir); if (distance < radius) { radius = distance; } @@ -591,7 +591,7 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const ix = nx()-1; for (int iy = 0; iy < ny(); ++iy) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - double distance = m_map.xy2dir(pixel).dist_deg(srcDir); + double distance = m_map.pix2dir(pixel).dist_deg(srcDir); if (distance < radius) { radius = distance; } @@ -841,7 +841,7 @@ void GLATEventCube::set_directions(void) for (int iy = 0; iy < ny(); ++iy) { for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - m_dirs.push_back(GLATInstDir(m_map.xy2dir(pixel))); + m_dirs.push_back(GLATInstDir(m_map.pix2dir(pixel))); m_omega.push_back(m_map.omega(pixel)); } } @@ -967,7 +967,7 @@ void GLATEventCube::set_bin(const int& index) // Set pointers m_bin.m_cube = this; - m_bin.m_counts = &(m_map.pixels()[index]); + m_bin.m_counts = const_cast(&(m_map.pixels()[index])); m_bin.m_energy = &(m_energies[m_bin.m_ieng]); m_bin.m_time = &m_time; m_bin.m_dir = &(m_dirs[m_bin.m_ipix]); diff --git a/inst/lat/src/GLATMeanPsf.cpp b/inst/lat/src/GLATMeanPsf.cpp index 18c4eed..2249767 100644 --- a/inst/lat/src/GLATMeanPsf.cpp +++ b/inst/lat/src/GLATMeanPsf.cpp @@ -738,7 +738,7 @@ void GLATMeanPsf::set_map_corrections(const GLATObservation& obs) // Compute offset angle in degrees GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); - double offset = cube->map().xy2dir(pixel).dist_deg(m_dir); + double offset = cube->map().pix2dir(pixel).dist_deg(m_dir); double omega = cube->map().omega(pixel); // Use only pixels within maximum PSF radius diff --git a/inst/lat/src/GLATResponse.cpp b/inst/lat/src/GLATResponse.cpp index 8586877..70c96e1 100644 --- a/inst/lat/src/GLATResponse.cpp +++ b/inst/lat/src/GLATResponse.cpp @@ -400,10 +400,10 @@ double GLATResponse::irf(const GLATEventBin& event, nodes.set_value(srcEng.log10MeV()); // Compute diffuse response - GSkymap* map = cube->diffrsp(idiff); - double* pixels = map->pixels() + event.ipix(); - rsp = nodes.wgt_left() * pixels[nodes.inx_left() * map->npix()] + - nodes.wgt_right() * pixels[nodes.inx_right() * map->npix()]; + GSkymap* map = cube->diffrsp(idiff); + const double* pixels = map->pixels() + event.ipix(); + rsp = nodes.wgt_left() * pixels[nodes.inx_left() * map->npix()] + + nodes.wgt_right() * pixels[nodes.inx_right() * map->npix()]; // Divide by solid angle and ontime since source maps are given in units of // counts/pixel/MeV. diff --git a/inst/lat/test/test_LAT.cpp b/inst/lat/test/test_LAT.cpp index 6f803f0..38c1980 100644 --- a/inst/lat/test/test_LAT.cpp +++ b/inst/lat/test/test_LAT.cpp @@ -287,7 +287,7 @@ void TestGLATLtCube::test_one_ltcube(const std::string& datadir, const double& r GLATLtCube ltcube(lat_ltcube); GEnergy energy; for (int i = 0; i < map.npix(); ++i) { - GSkyDir dir = map.pix2dir(i); + GSkyDir dir = map.inx2dir(i); map(i) = ltcube(dir, energy, test_fct1); } map.save(file1, true); @@ -304,7 +304,7 @@ void TestGLATLtCube::test_one_ltcube(const std::string& datadir, const double& r GLATLtCube ltcube(lat_ltcube); GEnergy energy; for (int i = 0; i < map.npix(); ++i) { - GSkyDir dir = map.pix2dir(i); + GSkyDir dir = map.inx2dir(i); map(i) = ltcube(dir, energy, test_fct2); } map.save(file2, true); diff --git a/pyext/GHealpix.i b/pyext/GHealpix.i index c20c79d..5c231a0 100644 --- a/pyext/GHealpix.i +++ b/pyext/GHealpix.i @@ -47,16 +47,14 @@ public: // Implemented pure virtual methods virtual void clear(void); virtual GHealpix* clone(void) const; + virtual int size(void) const; virtual std::string code(void) const; virtual std::string name(void) const; virtual void read(const GFitsHDU* hdu); virtual void write(GFitsHDU* hdu) const; - virtual double omega(const int& pix) const; - virtual double omega(const GSkyPixel& pix) const; - virtual GSkyDir pix2dir(const int& pix) const; - virtual int dir2pix(const GSkyDir& dir) const; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const; + virtual double omega(const GSkyPixel& pixel) const; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const; // Other methods const int& npix(void) const; diff --git a/pyext/GSkyProjection.i b/pyext/GSkyProjection.i index 2e5338d..40876c8 100644 --- a/pyext/GSkyProjection.i +++ b/pyext/GSkyProjection.i @@ -48,16 +48,14 @@ public: // Pure virtual methods (not implemented) virtual void clear(void) = 0; virtual GSkyProjection* clone(void) const = 0; + virtual int size(void) const = 0; virtual std::string code(void) const = 0; virtual std::string name(void) const = 0; virtual void read(const GFitsHDU* hdu) = 0; virtual void write(GFitsHDU* hdu) const = 0; - virtual double omega(const int& pix) const = 0; - virtual double omega(const GSkyPixel& pix) const = 0; - virtual GSkyDir pix2dir(const int& pix) const = 0; - virtual int dir2pix(const GSkyDir& dir) const = 0; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const = 0; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const = 0; + virtual double omega(const GSkyPixel& pixel) const = 0; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const = 0; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const = 0; // Virtual methods virtual std::string coordsys(void) const; diff --git a/pyext/GSkymap.i b/pyext/GSkymap.i index 776b767..6892900 100644 --- a/pyext/GSkymap.i +++ b/pyext/GSkymap.i @@ -127,34 +127,30 @@ public: GSkymap(const GSkymap& map); virtual ~GSkymap(void); - // 1D pixel methods - GSkyDir pix2dir(const int& pix) const; - int dir2pix(const GSkyDir& dir) const; - double omega(const int& pix) const; - - // 2D pixel methods - GSkyDir xy2dir(const GSkyPixel& pix) const; - GSkyPixel dir2xy(const GSkyDir& dir) const; - double omega(const GSkyPixel& pix) const; - // Methods - void clear(void); - GSkymap* clone(void) const; - void load(const std::string& filename); - void save(const std::string& filename, bool clobber = false) const; - void read(const GFitsHDU* hdu); - void write(GFits* file) const; - int npix(void) const; - int nx(void) const; - int ny(void) const; - int nmaps(void) const; - int xy2pix(const GSkyPixel& pix) const; - GSkyPixel pix2xy(const int& pix) const; - GSkyProjection* projection(void) const; - void projection(const GSkyProjection& proj); - double* pixels(void) const; - bool contains(const GSkyDir& dir) const; - bool contains(const GSkyPixel& pixel) const; + void clear(void); + GSkymap* clone(void) const; + const int& npix(void) const; + const int& nx(void) const; + const int& ny(void) const; + const int& nmaps(void) const; + GSkyPixel inx2pix(const int& index) const; + GSkyDir inx2dir(const int& index) const; + GSkyDir pix2dir(const GSkyPixel& pixel) const; + int pix2inx(const GSkyPixel& pixel) const; + int dir2inx(const GSkyDir& dir) const; + GSkyPixel dir2pix(const GSkyDir& dir) const; + double omega(const int& index) const; + double omega(const GSkyPixel& pixel) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyPixel& pixel) const; + const GSkyProjection* projection(void) const; + void projection(const GSkyProjection& proj); + const double* pixels(void) const; + void load(const std::string& filename); + void save(const std::string& filename, bool clobber = false) const; + void read(const GFitsHDU* hdu); + void write(GFits* file) const; }; diff --git a/pyext/GWcs.i b/pyext/GWcs.i index afe47f2..463738e 100644 --- a/pyext/GWcs.i +++ b/pyext/GWcs.i @@ -53,14 +53,12 @@ public: virtual std::string name(void) const = 0; // Implemented virtual methods + virtual int size(void) const; virtual void read(const GFitsHDU* hdu); virtual void write(GFitsHDU* hdu) const; - virtual double omega(const int& pix) const; - virtual double omega(const GSkyPixel& pix) const; - virtual GSkyDir pix2dir(const int& pix) const; - virtual int dir2pix(const GSkyDir& dir) const; - virtual GSkyDir xy2dir(const GSkyPixel& pix) const; - virtual GSkyPixel dir2xy(const GSkyDir& dir) const; + virtual double omega(const GSkyPixel& pixel) const; + virtual GSkyDir pix2dir(const GSkyPixel& pixel) const; + virtual GSkyPixel dir2pix(const GSkyDir& dir) const; // Other methods void set(const std::string& coords, diff --git a/src/model/GModelSpatialDiffuseCube.cpp b/src/model/GModelSpatialDiffuseCube.cpp index 8d44bb3..06a056c 100644 --- a/src/model/GModelSpatialDiffuseCube.cpp +++ b/src/model/GModelSpatialDiffuseCube.cpp @@ -446,15 +446,15 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy& energy, } } - // Convert 1D pixel index to 2D pixel index - GSkyPixel pixel = m_cube.pix2xy(low-offset); + // Convert sky map index to sky map pixel + GSkyPixel pixel = m_cube.inx2pix(low-offset); // Randomize pixel pixel.x(pixel.x() + ran.uniform() - 0.5); pixel.y(pixel.y() + ran.uniform() - 0.5); // Get sky direction - dir = m_cube.xy2dir(pixel); + dir = m_cube.pix2dir(pixel); } // endif: there were pixels in sky map diff --git a/src/model/GModelSpatialDiffuseMap.cpp b/src/model/GModelSpatialDiffuseMap.cpp index 7c1b064..890bdc8 100644 --- a/src/model/GModelSpatialDiffuseMap.cpp +++ b/src/model/GModelSpatialDiffuseMap.cpp @@ -350,15 +350,15 @@ GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy& energy, } } - // Convert 1D pixel index to 2D pixel index - GSkyPixel pixel = m_map.pix2xy(low); + // Convert sky map index to sky map pixel + GSkyPixel pixel = m_map.inx2pix(low); // Randomize pixel pixel.x(pixel.x() + ran.uniform() - 0.5); pixel.y(pixel.y() + ran.uniform() - 0.5); // Get sky direction - dir = m_map.xy2dir(pixel); + dir = m_map.pix2dir(pixel); } // endif: there were pixels in sky map diff --git a/src/sky/GHealpix.cpp b/src/sky/GHealpix.cpp index a18e1ae..d4c15ab 100644 --- a/src/sky/GHealpix.cpp +++ b/src/sky/GHealpix.cpp @@ -350,27 +350,14 @@ void GHealpix::write(GFitsHDU* hdu) const /***********************************************************************//** * @brief Returns solid angle of pixel * - * @param[in] pix Pixel number (0,1,...,num_pixels). + * @param[in] pixel Sky pixel. + * @return Solid angle of pixel. * - * HEALPix pixels have all the same solid angle, hence the pix argument is - * not used. + * Returns the solid angle of the specified @p pixel. Note that HEALPix + * pixels have all the same solid angle, hence the @p pixel argument is in + * fact not used by the method. ***************************************************************************/ -double GHealpix::omega(const int& pix) const -{ - // Return solid angle - return m_omega; -} - - -/***********************************************************************//** - * @brief Returns solid angle of pixel - * - * @param[in] pix Sky pixel. - * - * HEALPix pixels have all the same solid angle, hence the pix argument is - * not used. - ***************************************************************************/ -double GHealpix::omega(const GSkyPixel& pix) const +double GHealpix::omega(const GSkyPixel& pixel) const { // Return solid angle return m_omega; @@ -380,28 +367,35 @@ double GHealpix::omega(const GSkyPixel& pix) const /***********************************************************************//** * @brief Returns sky direction of pixel * - * @param[in] pix Pixel number (0,1,...,m_num_pixels). + * @param[in] pixel Sky map pixel. ***************************************************************************/ -GSkyDir GHealpix::pix2dir(const int& pix) const +GSkyDir GHealpix::pix2dir(const GSkyPixel& pixel) const { - // Declare result - GSkyDir result; - double theta = 0.0; - double phi = 0.0; + // Throw an exception if sky map pixel is not 1D + if (!pixel.is1D()) { + std::string msg = "Sky map pixel "+pixel.print()+" is not" + " 1-dimensional.\n" + "Only 1-dimensional pixels are supported by the" + " Healpix projection."; + throw GException::invalid_argument(G_XY2DIR, msg); + } // Perform ordering dependent conversion + double theta = 0.0; + double phi = 0.0; switch (m_ordering) { case 0: - pix2ang_ring(pix, &theta, &phi); + pix2ang_ring(int(pixel), &theta, &phi); break; case 1: - pix2ang_nest(pix, &theta, &phi); + pix2ang_nest(int(pixel), &theta, &phi); break; default: break; } // Store coordinate system dependent result + GSkyDir result; switch (m_coordsys) { case 0: result.radec(phi, gammalib::pihalf - theta); @@ -412,25 +406,23 @@ GSkyDir GHealpix::pix2dir(const int& pix) const default: break; } - - // Return result + + // Return return result; } /***********************************************************************//** - * @brief Returns pixel for a given sky direction + * @brief Returns sky map pixel of sky coordinate * - * @param[in] dir Sky direction. + * @param[in] dir Sky coordinate. + * @return Sky map pixel. ***************************************************************************/ -int GHealpix::dir2pix(const GSkyDir& dir) const +GSkyPixel GHealpix::dir2pix(const GSkyDir& dir) const { - // Declare result - int pix; + // Compute coordinate system dependent (z,phi) double z; double phi; - - // Compute coordinate system dependent (z,phi) switch (m_coordsys) { case 0: z = cos(gammalib::pihalf - dir.dec()); @@ -445,75 +437,20 @@ int GHealpix::dir2pix(const GSkyDir& dir) const } // Perform ordering dependent conversion + int index; switch (m_ordering) { case 0: - pix = ang2pix_z_phi_ring(z, phi); + index = ang2pix_z_phi_ring(z, phi); break; case 1: - pix = ang2pix_z_phi_nest(z, phi); + index = ang2pix_z_phi_nest(z, phi); break; default: break; } - // Return pixel index - return pix; -} - - -/***********************************************************************//** - * @brief Returns sky direction of pixel - * - * @param[in] pix Sky pixel. - * - * This dummy method throws an error when called as sky pixels are not - * implemented for a HealPix grid. Maybe we should implement them? - * - * @todo Think about implementation of sky pixel for HealPix data. We may - * use only the first argument, and even carry a usage flag in - * GSkyPixel - ***************************************************************************/ -GSkyDir GHealpix::xy2dir(const GSkyPixel& pix) const -{ - // Set error message - std::string message = "Method not defined for HPX projection."; - - // Throw error - throw GException::wcs(G_XY2DIR, message); - - // Set dummy return value - GSkyDir result; - - // Return - return result; -} - - -/***********************************************************************//** - * @brief Returns pixel of sky direction - * - * @param[in] dir Sky direction. - * - * This dummy method throws an error when called as sky pixels are not - * implemented for a HealPix grid. Maybe we should implement them? - * - * @todo Think about implementation of sky pixel for HealPix data. We may - * use only the first argument, and even carry a usage flag in - * GSkyPixel - ***************************************************************************/ -GSkyPixel GHealpix::dir2xy(const GSkyDir& dir) const -{ - // Set error message - std::string message = "Method not defined for HPX projection."; - - // Throw error - throw GException::wcs(G_DIR2XY2, message); - - // Set dummy return value - GSkyPixel result; - - // Return - return result; + // Return sky map pixel + return (GSkyPixel(index)); } diff --git a/src/sky/GSkymap.cpp b/src/sky/GSkymap.cpp index ce38f88..b5d6fb4 100644 --- a/src/sky/GSkymap.cpp +++ b/src/sky/GSkymap.cpp @@ -46,13 +46,13 @@ #define G_OP_ACCESS_1D "GSkymap::operator(int&, int&)" #define G_OP_ACCESS_2D "GSkymap::operator(GSkyPixel&, int&)" #define G_OP_VALUE "GSkymap::operator(GSkyDir&, int&)" -#define G_READ "GSkymap::read(const GFitsHDU*)" -#define G_PIX2DIR "GSkymap::pix2dir(int&)" +#define G_INX2DIR "GSkymap::inx2dir(int&)" +#define G_PIX2DIR "GSkymap::pix2dir(GSkyPixel&)" +#define G_DIR2INX "GSkymap::dir2inx(GSkyDir&)" #define G_DIR2PIX "GSkymap::dir2pix(GSkyDir&)" -#define G_XY2DIR "GSkymap::xy2dir(GSkyPixel&)" -#define G_DIR2XY "GSkymap::dir2xy(GSkyDir&)" #define G_OMEGA1 "GSkymap::omega(int&)" #define G_OMEGA2 "GSkymap::omega(GSkyPixel&)" +#define G_READ "GSkymap::read(const GFitsHDU*)" #define G_SET_WCS "GSkymap::set_wcs(std::string&, std::string&, double&,"\ " double&, double&, double&, double&, double&,"\ " GMatrix&, GVector&)" @@ -260,72 +260,104 @@ GSkymap::~GSkymap(void) ==========================================================================*/ /***********************************************************************//** - * @brief 1D pixel access operator + * @brief Assignment operator + * + * @param[in] map Sky map. + * @return Sky map. + ***************************************************************************/ +GSkymap& GSkymap::operator=(const GSkymap& map) +{ + // Execute only if object is not identical + if (this != &map) { + + // Free members + free_members(); + + // Initialise private members for clean destruction + init_members(); + + // Copy members + copy_members(map); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/***********************************************************************//** + * @brief Pixel index access operator * - * @param[in] pixel Pixel index (0,1,...,m_num_pixels). - * @param[in] map Map index (0,1,...,m_num_maps). + * @param[in] index Pixel index [0,...,npix()-1]. + * @param[in] map Map index [0,...,nmaps()-1]. + * @return Sky map pixel value. * * @exception GException::out_of_range * Pixel index and/or map index are outside valid range. * * Access sky map pixel by its index, where the most quickly varying axis is * the x axis of the map. - * - * @todo Implement proper skymap exception (actual is for matrix elements) ***************************************************************************/ -double& GSkymap::operator() (const int& pixel, const int& map) +double& GSkymap::operator()(const int& index, const int& map) { // Throw an error if pixel index or map index is not in valid range #if defined(G_RANGE_CHECK) - if (pixel < 0 || pixel >= m_num_pixels) { - throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1); + if (index < 0 || index >= m_num_pixels) { + throw GException::out_of_range(G_OP_ACCESS_1D, + "Sky map pixel index", + index, m_num_pixels); } if (map < 0 || map >= m_num_maps) { - throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1); + throw GException::out_of_range(G_OP_ACCESS_1D, + "Sky map map index", + map, m_num_maps); } #endif // Return reference to pixel value - return m_pixels[pixel+m_num_pixels*map]; + return m_pixels[index+m_num_pixels*map]; } /***********************************************************************//** - * @brief 1D pixel access operator (const variant) + * @brief Pixel index access operator (const variant) * - * @param[in] pixel Pixel index (0,1,...,m_num_pixels). - * @param[in] map Map index (0,1,...,m_num_maps). + * @param[in] index Pixel index [0,...,npix()-1]. + * @param[in] map Map index [0,...,nmaps()-1]. * * @exception GException::out_of_range * Pixel index and/or map index are outside valid range. * * Access sky map pixel by its index, where the most quickly varying axis is * the x axis of the map. - * - * @todo Implement proper skymap exception (actual is for matrix elements) ***************************************************************************/ -const double& GSkymap::operator() (const int& pixel, const int& map) const +const double& GSkymap::operator()(const int& index, const int& map) const { // Throw an error if pixel index or map index is not in valid range #if defined(G_RANGE_CHECK) - if (pixel < 0 || pixel >= m_num_pixels) { - throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1); + if (index < 0 || index >= m_num_pixels) { + throw GException::out_of_range(G_OP_ACCESS_1D, + "Sky map pixel index", + index, m_num_pixels); } if (map < 0 || map >= m_num_maps) { - throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1); + throw GException::out_of_range(G_OP_ACCESS_1D, + "Sky map map index", + map, m_num_maps); } #endif // Return reference to pixel value - return m_pixels[pixel+m_num_pixels*map]; + return m_pixels[index+m_num_pixels*map]; } /***********************************************************************//** - * @brief 2D pixel access operator + * @brief Sky map pixel access operator * - * @param[in] pixel Sky pixel. - * @param[in] map Map index (0,1,...,m_num_maps). + * @param[in] pixel Sky map pixel. + * @param[in] map Map index [0,...,nmaps()-1]. * * @exception GException::out_of_range * Sky pixel and/or map index are outside valid range. @@ -335,7 +367,7 @@ const double& GSkymap::operator() (const int& pixel, const int& map) const * * @todo Implement proper skymap exception (actual is for matrix elements) ***************************************************************************/ -double& GSkymap::operator() (const GSkyPixel& pixel, const int& map) +double& GSkymap::operator()(const GSkyPixel& pixel, const int& map) { // Throw an error if pixel index or map index is not in valid range #if defined(G_RANGE_CHECK) @@ -345,12 +377,14 @@ double& GSkymap::operator() (const GSkyPixel& pixel, const int& map) m_num_x-1, m_num_y-1); } if (map < 0 || map >= m_num_maps) { - throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1); + throw GException::out_of_range(G_OP_ACCESS_2D, + "Sky map map index", + map, m_num_maps); } #endif // Get pixel index - int index = xy2pix(pixel); + int index = pix2inx(pixel); // Return reference to pixel value return m_pixels[index+m_num_pixels*map]; @@ -358,10 +392,10 @@ double& GSkymap::operator() (const GSkyPixel& pixel, const int& map) /***********************************************************************//** - * @brief 2D pixel access operator + * @brief Sky map pixel access operator * - * @param[in] pixel Sky pixel. - * @param[in] map Map index (0,1,...,m_num_maps). + * @param[in] pixel Sky map pixel. + * @param[in] map Map index [0,...,nmaps()-1]. * * @exception GException::out_of_range * Sky pixel and/or map index are outside valid range. @@ -381,12 +415,14 @@ const double& GSkymap::operator()(const GSkyPixel& pixel, const int& map) const m_num_x-1, m_num_y-1); } if (map < 0 || map >= m_num_maps) { - throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1); + throw GException::out_of_range(G_OP_ACCESS_2D, + "Sky map map index", + map, m_num_maps); } #endif // Get pixel index - int index = xy2pix(pixel); + int index = pix2inx(pixel); // Return reference to pixel value return m_pixels[index+m_num_pixels*map]; @@ -397,7 +433,7 @@ const double& GSkymap::operator()(const GSkyPixel& pixel, const int& map) const * @brief Return interpolated skymap value for sky direction * * @param[in] dir Sky direction. - * @param[in] map Map index (0,1,...,m_num_maps). + * @param[in] map Map index [0,...,nmaps()-1]. * * @exception GException::out_of_range * Map index lies outside valid range. @@ -406,12 +442,14 @@ const double& GSkymap::operator()(const GSkyPixel& pixel, const int& map) const * interpolation of the neighbouring pixels. If the sky direction falls * outside the area covered by the skymap, a value of 0 is returned. ***************************************************************************/ -double GSkymap::operator() (const GSkyDir& dir, const int& map) const +double GSkymap::operator()(const GSkyDir& dir, const int& map) const { // Throw an error if the map index is not in valid range #if defined(G_RANGE_CHECK) if (map < 0 || map >= m_num_maps) { - throw GException::out_of_range(G_OP_VALUE, map, 0, m_num_maps-1); + throw GException::out_of_range(G_OP_VALUE, + "Sky map map index", + map, m_num_maps); } #endif @@ -419,7 +457,7 @@ double GSkymap::operator() (const GSkyDir& dir, const int& map) const double intensity = 0.0; // Determine sky pixel - GSkyPixel pixel = dir2xy(dir); + GSkyPixel pixel = dir2pix(dir); // Continue only if pixel is within the map if (contains(pixel)) { @@ -476,32 +514,6 @@ double GSkymap::operator() (const GSkyDir& dir, const int& map) const } -/***********************************************************************//** - * @brief Assignment operator - * - * @param[in] map Sky map. - ***************************************************************************/ -GSkymap& GSkymap::operator= (const GSkymap& map) -{ - // Execute only if object is not identical - if (this != &map) { - - // Free members - free_members(); - - // Initialise private members for clean destruction - init_members(); - - // Copy members - copy_members(map); - - } // endif: object was not identical - - // Return this object - return *this; -} - - /*========================================================================== = = = Public methods = @@ -511,7 +523,7 @@ GSkymap& GSkymap::operator= (const GSkymap& map) /***********************************************************************//** * @brief Clear instance. * - * This method properly resets the object to an initial state. + * Resets the sky map to the initial state. ***************************************************************************/ void GSkymap::clear(void) { @@ -527,7 +539,9 @@ void GSkymap::clear(void) /***********************************************************************//** - * @brief Clone instance + * @brief Clone sky map + * + * @return Pointer to deep copy of sky map. ***************************************************************************/ GSkymap* GSkymap::clone(void) const { @@ -535,6 +549,368 @@ GSkymap* GSkymap::clone(void) const } +/***********************************************************************//** + * @brief Converts pixel index into sky map pixel + * + * @param[in] index Pixel index [0,...,npix()-1]. + * @return Sky map pixel. + * + * Converts the pixel @p index into a sky map pixel (GSkyPixel).The dimension + * of GSkyPixel will be identical to the dimension of the sky map (i.e. a 1D + * sky map leads to a 1D GSkyPixel object, a 2D sky map leads to a 2D + * GSkyPixel object. + ***************************************************************************/ +GSkyPixel GSkymap::inx2pix(const int& index) const +{ + // Initialise sky map pixel + GSkyPixel pixel; + + // Get x and y indices + double x; + double y; + if (m_num_x != 0) { //!< 2D sky map + pixel.x(double(index % m_num_x)); + pixel.y(double(index / m_num_x)); + } + else { //!< 1D sky map + pixel.index(index); + } + + // Return pixel + return pixel; +} + + +/***********************************************************************//** + * @brief Returns sky direction of pixel + * + * @param[in] index Pixel index [0,...,npix()-1]. + * @return Sky direction. + * + * @exception GException::invalid_value + * No valid sky projection found. + * + * Returns sky direction for a given pixel index. + ***************************************************************************/ +GSkyDir GSkymap::inx2dir(const int& index) const +{ + // Throw error if sky projection is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_INX2DIR, msg); + } + + // Determine sky direction from pixel index. + GSkyDir dir = m_proj->pix2dir(inx2pix(index)); + + // Return sky direction + return dir; +} + + +/***********************************************************************//** + * @brief Returns sky direction of pixel + * + * @param[in] pixel Sky map pixel. + * @return Sky direction. + * + * @exception GException::invalid_value + * No valid sky projection found. + * @exception GException::invalid_argument + * 2D sky map pixel used to access 1D projection. + * + * Returns sky direction for a given sky map @p pixel. + ***************************************************************************/ +GSkyDir GSkymap::pix2dir(const GSkyPixel& pixel) const +{ + // Throw error if WCS is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_PIX2DIR, msg); + } + + // Initialise sky direction + GSkyDir dir; + + // If pixel size matches the projection size then perform a straight + // forward conversion + if (m_proj->size() == pixel.size()) { + dir = m_proj->pix2dir(pixel); + } + + // ... otherwise, if we have a 2D projection but a 1D pixel then + // interpret the pixel as the linear index in the pixel array + else if (m_proj->size() == 2) { + dir = m_proj->pix2dir(GSkyPixel(inx2pix(int(pixel)))); + } + + // ... otherwise we have a 1D projection but a 2D pixel. There is + // no unambiguous way to handle this case, hence we throw an exception + else { + std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+ + " is used to determine the sky direction for" + " the 1-dimensional sky projection \""+ + m_proj->name()+"\"\n" + "Please specify a 1-dimensional sky map pixel."; + throw GException::invalid_argument(G_PIX2DIR, msg); + } + + // Return sky direction + return dir; +} + + +/***********************************************************************//** + * @brief Converts sky map pixel into pixel index + * + * @param[in] pixel Sky map pixel. + * @return Pixel index [0,...,npix()-1]. + * + * Converts a sky map @p pixel into the pixel index. + ***************************************************************************/ +int GSkymap::pix2inx(const GSkyPixel& pixel) const +{ + // Initialise pixel index + int index = 0; + + // Handle 1D sky map pixel + if (pixel.is1D()) { + index = int(pixel); + } + + // Handle 2D sky map pixel + else if (pixel.is2D()) { + + // Get x and y indices by rounding the (x,y) values + int ix = int(pixel.x()+0.5); + int iy = int(pixel.y()+0.5); + + // Set index + index = ix + iy * m_num_x; + + } + + // Return index + return index; +} + + +/***********************************************************************//** + * @brief Returns pixel index for a given sky direction + * + * @param[in] dir Sky direction. + * @return Pixel index [0,...,npix()-1]. + * + * @exception GException::invalid_value + * No valid sky projection found. + * + * Returns sky map pixel index for a given sky direction. + ***************************************************************************/ +int GSkymap::dir2inx(const GSkyDir& dir) const +{ + // Throw error if WCS is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_DIR2INX, msg); + } + + // Determine pixel index for a given sky direction + int index = pix2inx(m_proj->dir2pix(dir)); + + // Return pixel index + return index; +} + + +/***********************************************************************//** + * @brief Returns sky map pixel for a given sky direction + * + * @param[in] dir Sky direction. + * @return Sky map pixel. + * + * @exception GException::invalid_value + * No valid sky projection found. + * + * Returns sky map pixel for a given sky direction. + ***************************************************************************/ +GSkyPixel GSkymap::dir2pix(const GSkyDir& dir) const +{ + // Throw error if WCS is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_DIR2PIX, msg); + } + + // Determine pixel for a given sky direction + GSkyPixel pixel = m_proj->dir2pix(dir); + + // Return pixel + return pixel; +} + + +/***********************************************************************//** + * @brief Returns solid angle of pixel + * + * @param[in] index Pixel index [0,...,npix()-1]. + * @return Solid angle (steradians) + * + * @exception GException::invalid_value + * No valid sky projection found. + * + * Returns the solid angle of the pixel with the specified @p index. + ***************************************************************************/ +double GSkymap::omega(const int& index) const +{ + // Throw error if WCS is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_OMEGA1, msg); + } + + // Determine solid angle from pixel index. + double omega = m_proj->omega(inx2pix(index)); + + // Return solid angle + return omega; +} + + +/***********************************************************************//** + * @brief Returns solid angle of pixel + * + * @param[in] pixel Sky map pixel. + * @return Solid angle (steradians) + * + * @exception GException::invalid_value + * No valid sky projection found. + * + * Returns the solid angle of the specified sky map @p pixel. + ***************************************************************************/ +double GSkymap::omega(const GSkyPixel& pixel) const +{ + // Throw error if WCS is not valid + if (m_proj == NULL) { + std::string msg = "Sky projection has not been defined."; + throw GException::invalid_value(G_OMEGA2, msg); + } + + // Initialise solid angle + double omega = 0.0; + + // If pixel size matches the projection size then perform a straight + // forward solid angle determination + if (m_proj->size() == pixel.size()) { + omega = m_proj->omega(pixel); + } + + // ... otherwise, if we have a 2D projection but a 1D pixel then + // interpret the pixel as the linear index in the pixel array + else if (m_proj->size() == 2) { + omega = m_proj->omega(GSkyPixel(inx2pix(int(pixel)))); + } + + // ... otherwise we have a 1D projection but a 2D pixel. There is + // no unambiguous way to handle this case, hence we throw an exception + else { + std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+ + " is used to determine the solid angle for" + " the 1-dimensional sky projection \""+ + m_proj->name()+"\"\n" + "Please specify a 1-dimensional sky map pixel."; + throw GException::invalid_argument(G_OMEGA2, msg); + } + + // Return solid angle + return omega; +} + + +/***********************************************************************//** + * @brief Set sky projection + * + * @param[in] proj Sky projection. + * + * Sets the projection from celestial to pixel coordinates. The method + * performs a deep copy of @p proj, allowing to destroy the argument after + * using the method. + * + * Warning: this method may corrupt the GSkymap object as it allows assigning + * for example a 1D projection to a 2D skymap. Please use this method only + * when you know what you're doing. + * + * @todo We may restrict this method to not allow changing the projection + * dimension. + ***************************************************************************/ +void GSkymap::projection(const GSkyProjection& proj) +{ + // Free any existing WCS + if (m_proj != NULL) delete m_proj; + + // Clone input WCS + m_proj = proj.clone(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Verifies if sky direction falls in map + * + * @param[in] dir Sky direction. + * + * This method checks if the specified sky direction falls within the pixels + * covered by the skymap. The method uses the dir2xy method to convert the + * sky direction into 2D pixel indices, and then checks whether the pixel + * indices fall in the skymap. + ***************************************************************************/ +bool GSkymap::contains(const GSkyDir& dir) const +{ + // Convert sky direction into sky pixel + GSkyPixel pixel = dir2pix(dir); + + // Return location flag + return (contains(pixel)); +} + + +/***********************************************************************//** + * @brief Checks if sky map pixel falls in map + * + * @param[in] pixel Sky map pixel. + * @return Trus if pixels is within map, false otherwise. + * + * Checks whether the specified sky map @p pixel falls within the skymap + * or not. + ***************************************************************************/ +bool GSkymap::contains(const GSkyPixel& pixel) const +{ + // Initialise containment flag + bool inmap = false; + + // Test 1D pixels + if (pixel.is1D()) { + if (pixel.index()+0.5 >= 0.0 && pixel.index()-0.5 < m_num_pixels) { + inmap = true; + } + } + + // Test 2D pixels + else if (pixel.is2D()) { + + // If pixel is in range then set containment flag to true + if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) && + (pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) { + inmap = true; + } + + } + + // Return containment flag + return inmap; +} + /***********************************************************************//** * @brief Load skymap from FITS file. @@ -544,6 +920,9 @@ GSkymap* GSkymap::clone(void) const * Loads HEALPix and non HEALPix skymaps. First searches for HEALPix map in * FITS file by scanning all HDUs for PIXTYPE=HEALPIX. If no HEALPix map has * been found then search load first non-empty image. + * + * @todo Do we have to restrict a HEALPix map to a BinTable and a WCS map + * to a Double precision image??? ***************************************************************************/ void GSkymap::load(const std::string& filename) { @@ -568,16 +947,12 @@ void GSkymap::load(const std::string& filename) // Get pointer to HDU hdu = fits.hdu(extno); - + // If PIXTYPE keyword equals "HEALPIX" then load map - try { - if (hdu->string("PIXTYPE") == "HEALPIX") { - read_healpix(static_cast(hdu)); - loaded = true; - break; - } - } - catch (GException::fits_key_not_found &e) { + if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") { + read_healpix(static_cast(hdu)); + loaded = true; + break; } } // endfor: looped over HDUs @@ -597,7 +972,7 @@ void GSkymap::load(const std::string& filename) } // Load WCS map - read_wcs(static_cast(hdu)); + read_wcs(static_cast(hdu)); loaded = true; break; @@ -674,14 +1049,10 @@ void GSkymap::read(const GFitsHDU* hdu) // Initialize load flag bool loaded = false; - // Try load as HEALPix map - try { - if (hdu->string("PIXTYPE") == "HEALPIX") { - read_healpix(static_cast(hdu)); - loaded = true; - } - } - catch (GException::fits_key_not_found &e) { + // If PIXTYPE keyword equals "HEALPIX" then load map + if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") { + read_healpix(static_cast(hdu)); + loaded = true; } // ... otherwise try loading as non HEALPix map @@ -689,7 +1060,7 @@ void GSkymap::read(const GFitsHDU* hdu) // Load only if HDU contains an image if (hdu->exttype() == 0) { - read_wcs(static_cast(hdu)); + read_wcs(static_cast(hdu)); loaded = true; } @@ -741,332 +1112,6 @@ void GSkymap::write(GFits* file) const /***********************************************************************//** - * @brief Returns sky direction of pixel - * - * @param[in] pix Pixel number (0,1,...,m_num_pixels). - * - * @exception GException::wcs - * No valid WCS found. - * - * Returns sky direction for a given sky map pixel. This methods works for - * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme. - * Automatic index conversion is provided so that 1D schemes need only to - * implement 1D conversion methods while 2D schemes need only to implement - * 2D conversion methods. - ***************************************************************************/ -GSkyDir GSkymap::pix2dir(const int& pix) const -{ - // Throw error if sky projection is not valid - if (m_proj == NULL) { - throw GException::wcs(G_PIX2DIR, "No valid sky projection found."); - } - - // Determine sky direction from pixel. Use 2D version if sky map is - // 2D, otherwise use 1D version. - GSkyDir dir = (m_num_x == 0) ? m_proj->pix2dir(pix) - : m_proj->xy2dir(pix2xy(pix)); - - // Return sky direction - return dir; -} - - -/***********************************************************************//** - * @brief Returns pixel index for a given sky direction - * - * @param[in] dir Sky direction. - * - * @exception GException::wcs - * No valid WCS found. - * - * Returns sky map pixel for a given sky direction. This methods works for - * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme. - * Automatic index conversion is provided so that 1D schemes need only to - * implement 1D conversion methods while 2D schemes need only to implement - * 2D conversion methods. - ***************************************************************************/ -int GSkymap::dir2pix(const GSkyDir& dir) const -{ - // Throw error if WCS is not valid - if (m_proj == NULL) { - throw GException::wcs(G_DIR2PIX, "No valid sky projection found."); - } - - // Determine 1D pixel index for a given sky direction - int pix = (m_num_x == 0) ? m_proj->dir2pix(dir) - : xy2pix(m_proj->dir2xy(dir)); - - // Return pixel index - return pix; -} - - -/***********************************************************************//** - * @brief Returns sky direction of pixel - * - * @param[in] pix Sky map pixel. - * - * @exception GException::wcs - * No valid WCS found. - * - * Returns sky direction for a given sky map pixel. This methods works for - * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme. - * Automatic index conversion is provided so that 1D schemes need only to - * implement 1D conversion methods while 2D schemes need only to implement - * 2D conversion methods. - ***************************************************************************/ -GSkyDir GSkymap::xy2dir(const GSkyPixel& pix) const -{ - // Throw error if WCS is not valid - if (m_proj == NULL) { - throw GException::wcs(G_XY2DIR, "No valid sky projection found."); - } - - // Determine sky direction from pixel. Use 2D version if sky map is - // 2D, otherwise use 1D version. - GSkyDir dir = (m_num_x == 0) ? m_proj->pix2dir(xy2pix(pix)) - : m_proj->xy2dir(pix); - - // Return sky direction - return dir; -} - - -/***********************************************************************//** - * @brief Returns sky map pixel for a given sky direction - * - * @param[in] dir Sky direction. - * - * @exception GException::wcs - * No valid WCS found. - * - * Returns sky map pixel for a given sky direction. This methods works for - * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme. - * Automatic index conversion is provided so that 1D schemes need only to - * implement 1D conversion methods while 2D schemes need only to implement - * 2D conversion methods. - ***************************************************************************/ -GSkyPixel GSkymap::dir2xy(const GSkyDir& dir) const -{ - // Throw error if WCS is not valid - if (m_proj == NULL) { - throw GException::wcs(G_DIR2XY, "No valid sky projection found."); - } - - // Determine 1D pixel index for a given sky direction - GSkyPixel pix = (m_num_x == 0) ? pix2xy(m_proj->dir2pix(dir)) - : m_proj->dir2xy(dir); - - // Return pixel index - return pix; -} - - -/***********************************************************************//** - * @brief Returns solid angle of pixel - * - * @param[in] pix Pixel number (0,1,...,m_num_pixels). - * - * @exception GException::wcs - * No valid WCS found. - ***************************************************************************/ -double GSkymap::omega(const int& pix) const -{ - // Throw error if WCS is not valid - if (m_proj == NULL) { - throw GException::wcs(G_OMEGA1, "No valid sky projection found."); - } - - // Determine solid angle from pixel. Use 2D version if sky map is - // 2D, otherwise use 1D version. - double omega = (m_num_x == 0) ? m_proj->omega(pix) - : m_proj->omega(pix2xy(pix)); - - // Return solid angle - return omega; -} - - -/***********************************************************************//** - * @brief Returns solid angle of pixel - * - * @param[in] pix Sky map pixel. - * - * @exception GException::wcs - * No valid WCS found. - ***************************************************************************/ -double GSkymap::omega(const GSkyPixel& pix) const -{ - // Throw error if WCS is not valid - if (m_proj == NULL) { - throw GException::wcs(G_OMEGA2, "No valid sky projection found."); - } - - // Determine solid angle from pixel. Use 2D version if sky map is - // 2D, otherwise use 1D version. - double omega = (m_num_x == 0) ? m_proj->omega(xy2pix(pix)) - : m_proj->omega(pix); - - // Return solid angle - return omega; -} - - -/***********************************************************************//** - * @brief Returns number of pixels - ***************************************************************************/ -int GSkymap::npix(void) const -{ - // Return number of pixels - return m_num_pixels; -} - - -/***********************************************************************//** - * @brief Returns number of pixels in x coordinate - ***************************************************************************/ -int GSkymap::nx(void) const -{ - // Return number of pixels - return m_num_x; -} - - -/***********************************************************************//** - * @brief Returns number of pixels in y coordinate - ***************************************************************************/ -int GSkymap::ny(void) const -{ - // Return number of pixels - return m_num_y; -} - - -/***********************************************************************//** - * @brief Returns number of maps - ***************************************************************************/ -int GSkymap::nmaps(void) const -{ - // Return number of maps - return m_num_maps; -} - - -/***********************************************************************//** - * @brief Converts 2D index (x,y) into 1D pixel index - * - * @param[in] pix 2D pixel index. - * - * The (x,y) value is rounded to nearest integers before conversion. The x - * axis is assumed as the most rapidely varying index. - ***************************************************************************/ -int GSkymap::xy2pix(const GSkyPixel& pix) const -{ - // Get x and y indices by rounding the (x,y) values - int ix = int(pix.x()+0.5); - int iy = int(pix.y()+0.5); - - // Return index - return (ix+iy*m_num_x); -} - - -/***********************************************************************//** - * @brief Converts 1D pixel index into 2D index (x,y) - * - * @param[in] pix 1D pixel index. - * - * If skymap is not 2D the pixel index is returned in the x element of the - * GSkyPixel object. - ***************************************************************************/ -GSkyPixel GSkymap::pix2xy(const int& pix) const -{ - // Get x and y indices - double x; - double y; - if (m_num_x != 0) { - x = double(pix % m_num_x); - y = double(pix / m_num_x); - } - else { - x = double(pix); - y = 0.0; - } - - // Set pixel - GSkyPixel pixel(x, y); - - // Return pixel - return pixel; -} - - -/***********************************************************************//** - * @brief Set sky projection - * - * @param[in] proj Sky projection. - * - * Sets the projection from celestial to pixel coordinates. The method - * performs a deep copy of @p proj, allowing to destroy the argument after - * using the method. - ***************************************************************************/ -void GSkymap::projection(const GSkyProjection& proj) -{ - // Free any existing WCS - if (m_proj != NULL) delete m_proj; - - // Clone input WCS - m_proj = proj.clone(); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Verifies if sky direction falls in map - * - * @param[in] dir Sky direction. - * - * This method checks if the specified sky direction falls within the pixels - * covered by the skymap. The method uses the dir2xy method to convert the - * sky direction into 2D pixel indices, and then checks whether the pixel - * indices fall in the skymap. - ***************************************************************************/ -bool GSkymap::contains(const GSkyDir& dir) const -{ - // Convert sky direction into sky pixel - GSkyPixel pixel = dir2xy(dir); - - // Return location flag - return (contains(pixel)); -} - - -/***********************************************************************//** - * @brief Verifies if sky pixel falls in map - * - * @param[in] pixel Sky pixel. - * - * This method checks if the specified sky pixel is within the skymap. - ***************************************************************************/ -bool GSkymap::contains(const GSkyPixel& pixel) const -{ - // Initialise location flag - bool inmap = false; - - // If pixel is in range then set location flag to true - if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) && - (pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) { - inmap = true; - } - - // Return location flag - return inmap; -} - - -/***********************************************************************//** * @brief Print sky map * * @param[in] chatter Chattiness (defaults to NORMAL). @@ -1088,7 +1133,7 @@ std::string GSkymap::print(const GChatter& chatter) const result.append(gammalib::str(m_num_pixels)); result.append("\n"+gammalib::parformat("Number of maps")); result.append(gammalib::str(m_num_maps)); - if (m_proj != NULL && m_proj->code() != "HPX") { + if (m_proj != NULL && m_proj->size() == 2) { result.append("\n"+gammalib::parformat("X axis dimension")); result.append(gammalib::str(m_num_x)); result.append("\n"+gammalib::parformat("Y axis dimension")); @@ -1181,7 +1226,7 @@ void GSkymap::copy_members(const GSkymap& map) // Copy pixels if (size > 0 && map.m_pixels != NULL) { alloc_pixels(); - for (int i = 0; i < size; ++i) { + for (int i = 0; i < size; ++i) { m_pixels[i] = map.m_pixels[i]; } } diff --git a/src/sky/GWcs.cpp b/src/sky/GWcs.cpp index 871f5e0..4378c63 100644 --- a/src/sky/GWcs.cpp +++ b/src/sky/GWcs.cpp @@ -436,35 +436,9 @@ void GWcs::write(GFitsHDU* hdu) const /***********************************************************************//** * @brief Returns solid angle of pixel in units of steradians * - * @param[in] pix Pixel index. + * @param[in] pixel Pixel index (x,y) * - * This dummy method throws an error when called as pixel indices are not - * implemented for a HealPix grid. Maybe we should drop this method and - * implement a general GSkyPixel class tha also handles HealPix? - * - * @todo Think about implementation of sky pixel for HealPix data. We may - * use only the first argument, and even carry a usage flag in - * GSkyPixel - ***************************************************************************/ -double GWcs::omega(const int& pix) const -{ - // Set error message - std::string message = "Method not defined for general WCS projection."; - - // Throw error - throw GException::wcs(G_OMEGA, message); - - // Return - return 0.0; -} - - -/***********************************************************************//** - * @brief Returns solid angle of pixel in units of steradians - * - * @param[in] pix Pixel index (x,y) - * - * Estimate solid angles of pixel by compuing the coordinates in the 4 pixel + * Estimate solid angles of pixel by computing the coordinates in the 4 pixel * corners. The surface is computed using a cartesian approximation: * a * 1-----5-----2 a+b @@ -475,15 +449,15 @@ double GWcs::omega(const int& pix) const * rotated sky maps. Something more intelligent should be implemented in * the future. ***************************************************************************/ -double GWcs::omega(const GSkyPixel& pix) const +double GWcs::omega(const GSkyPixel& pixel) const { // Get the sky directions of the 6 points - GSkyDir dir1 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()-0.5)); - GSkyDir dir2 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()-0.5)); - GSkyDir dir3 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()+0.5)); - GSkyDir dir4 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()+0.5)); - GSkyDir dir5 = xy2dir(GSkyPixel(pix.x(), pix.y()-0.5)); - GSkyDir dir6 = xy2dir(GSkyPixel(pix.x(), pix.y()+0.5)); + GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5)); + GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5)); + GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5)); + GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5)); + GSkyDir dir5 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5)); + GSkyDir dir6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5)); // Compute distances between sky directions double a = dir1.dist(dir2); @@ -499,64 +473,16 @@ double GWcs::omega(const GSkyPixel& pix) const /***********************************************************************//** - * @brief Returns sky direction of pixel (DUMMY METHOD) - * - * @param[in] pix Pixel number (0,1,...,m_num_pixels). - * - * @exception GException::wcs - * Method not defined for projection. - * - * This method should be overloaded by the appropriate method of the derived - * class. It thus should never get called. - ***************************************************************************/ -GSkyDir GWcs::pix2dir(const int& pix) const -{ - // Set error message - std::string message = "Method not defined for general WCS projection."; - - // Throw error - throw GException::wcs(G_PIX2DIR, message); - - // Define direction - GSkyDir dir; - - // Return - return dir; -} - - -/***********************************************************************//** - * @brief Returns pixel for sky direction (DUMMY METHOD) - * - * @param[in] dir Sky direction. - * - * @exception GException::wcs - * Method not defined for projection. - * - * This method should be overloaded by the appropriate method of the derived - * class. It thus should never get called. - ***************************************************************************/ -int GWcs::dir2pix(const GSkyDir& dir) const -{ - // Set error message - std::string message = "Method not defined for general WCS projection."; - - // Throw error - throw GException::wcs(G_DIR2PIX, message); - - // Return - return 0; -} - - -/***********************************************************************//** - * @brief Returns sky direction of pixel + * @brief Returns sky direction of sky map pixel * - * @param[in] pix Sky pixel. + * @param[in] pixel Sky map pixel. + * @return Sky direction. * - * Note that pixel indices start from 0. + * Returns the sky direction of a sky map @p pixel. Note that the sky + * map pixel values start from 0 while the WCS pixel reference starts from + * 1. ***************************************************************************/ -GSkyDir GWcs::xy2dir(const GSkyPixel& pix) const +GSkyDir GWcs::pix2dir(const GSkyPixel& pixel) const { // Allocate memory for transformation double pixcrd[2]; @@ -568,8 +494,8 @@ GSkyDir GWcs::xy2dir(const GSkyPixel& pix) const // Set sky pixel. We have to add 1.0 here as the WCS pixel reference // (CRPIX) starts from one while GSkyPixel starts from 0. - pixcrd[0] = pix.x() + 1.0; - pixcrd[1] = pix.y() + 1.0; + pixcrd[0] = pixel.x() + 1.0; + pixcrd[1] = pixel.y() + 1.0; // Trasform pixel-to-world coordinate wcs_p2s(1, 2, pixcrd, imgcrd, &phi, &theta, world, &stat); @@ -585,7 +511,7 @@ GSkyDir GWcs::xy2dir(const GSkyPixel& pix) const // Debug: Dump transformation steps #if defined(G_XY2DIR_DEBUG) - std::cout << "xy2dir: pixel=" << pix + std::cout << "xy2dir: pixel=" << pixel << " (x,y)=(" << pixcrd[0] << "," << pixcrd[1] << ")" << " (phi,theta)=(" << phi << "," << theta << ")" << " (lng,lat)=(" << world[0] << "," << world[1] << ")" << std::endl; @@ -597,14 +523,16 @@ GSkyDir GWcs::xy2dir(const GSkyPixel& pix) const /***********************************************************************//** - * @brief Returns pixel of sky direction + * @brief Returns sky map pixel of sky direction * * @param[in] dir Sky direction. + * @return Sky map pixel. * - * Note that GSkyPixel indices start from 0 while the WCS pixel reference - * starts from 1. + * Returns the sky map pixel for a given sky direction. Note that the sky + * map pixel values start from 0 while the WCS pixel reference starts from + * 1. ***************************************************************************/ -GSkyPixel GWcs::dir2xy(const GSkyDir& dir) const +GSkyPixel GWcs::dir2pix(const GSkyDir& dir) const { // Allocate memory for transformation double pixcrd[2]; diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 465b0f9..3f78196 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -106,10 +106,10 @@ double TestGSky::wcs_forth_back_pixel(GWcs* wcs, int nx, int ny, double& crpix1, GSkyPixel pix_in(x,y); // Forth: Transform to world - GSkyDir dir = wcs->xy2dir(pix_in); + GSkyDir dir = wcs->pix2dir(pix_in); // Back: Transform to pixel - GSkyPixel pix_out = wcs->dir2xy(dir); + GSkyPixel pix_out = wcs->dir2pix(dir); // Compute distance double dx = pix_in.x()-pix_out.x(); @@ -177,7 +177,7 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double& crpix1, double& crp for (int iy = -ny; iy <= ny; ++iy) { // Set y value - double y = double(iy) + crpix2; + double y = double(iy) + crpix2; // Skip pixels outside valid range (they are not expected to // transform bijectively) @@ -189,8 +189,8 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double& crpix1, double& crp GSkyPixel pix_in(x,y); // Transform to world - GSkyDir dir1 = wcs->xy2dir(pix_in); - GSkyDir dir2 = cpy->xy2dir(pix_in); + GSkyDir dir1 = wcs->pix2dir(pix_in); + GSkyDir dir2 = cpy->pix2dir(pix_in); double angle = dir1.dist_deg(dir2); if (angle > angle_max) { angle_max = angle; @@ -207,8 +207,8 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double& crpix1, double& crp #endif // Transform to pixel - GSkyPixel pix_out1 = wcs->dir2xy(dir1); - GSkyPixel pix_out2 = cpy->dir2xy(dir2); + GSkyPixel pix_out1 = wcs->dir2pix(dir1); + GSkyPixel pix_out2 = cpy->dir2pix(dir2); // Compute distance double dx = pix_out1.x()-pix_out2.x(); @@ -715,8 +715,8 @@ void TestGSky::test_GSkymap_wcs_construct(void) for (int l = -180; l < 180; ++l) { for (int b = -90; b < 90; ++b) { dir.lb_deg(double(l),double(b)); - GSkyPixel pixel = map1.dir2xy(dir); - GSkyDir dir_back = map1.xy2dir(pixel); + GSkyPixel pixel = map1.dir2pix(dir); + GSkyDir dir_back = map1.pix2dir(pixel); double dist = dir.dist_deg(dir_back); if (dist > eps) { throw exception_failure("Sky direction differs: dir="+dir.print()+" pixel="+pixel.print()+" dir_back"+ dir_back.print()+" dist="+gammalib::str(dist)+" deg"); diff --git a/test/test_GSky.py b/test/test_GSky.py index a3aea45..719ed7a 100755 --- a/test/test_GSky.py +++ b/test/test_GSky.py @@ -73,10 +73,11 @@ class Test(GPythonTestSuite): # Append tests self.append(self.test_skymap_healpix, "Test HEALPix map") self.append(self.test_skymap_ait, "Test AIT 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_tan, "Test TAN projection map") + self.append(self.test_skymap_mer, "Test MER projection map") self.append(self.test_skymap_stg, "Test STG projection map") - self.append(self.test_skymap_azp, "Test AZP projection map") + self.append(self.test_skymap_tan, "Test TAN projection map") self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion") # Return @@ -87,21 +88,21 @@ class Test(GPythonTestSuite): """ Control coordinate and pixel transformation methods. """ - # Control pix2dir / dir2pix methods for all pixels + # Control inx2dir/dir2inx methods for all pixels for i in range(pixels.npix()): - dir = pixels.pix2dir(i) - ipix = pixels.dir2pix(dir) - msg = map + " pix2dir/dir2pix check for pixel " + str(i) - err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(ipix) + \ + dir = pixels.inx2dir(i) + inx = pixels.dir2inx(dir) + msg = map + " inx2dir/dir2inx check for pixel " + str(i) + err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(inx) + \ "), RA=" + str(dir.ra() * 180 / pi) + ", Dec=" + str(dir.dec() * 180 / pi) - self.test_assert(i == ipix, msg, err) + self.test_assert(i == inx, msg, err) # Control SkyDir coordinate transformation for all pixels for i in range(pixels.npix()): - dir = pixels.pix2dir(i) + dir = pixels.inx2dir(i) dir_new = GSkyDir() dir_new.lb(dir.l(), dir.b()) - dra = abs(dir.ra() - dir_new.ra()) + dra = abs(dir.ra() - dir_new.ra()) if (dra >= 5.0): dra -= 2.0 * pi ddec = abs(dir.dec() - dir_new.dec()) @@ -196,6 +197,17 @@ class Test(GPythonTestSuite): # Return return + # Test AZP projection + def test_skymap_azp(self): + """ + Test AZP projection. + """ + # Execute generic test + self.test_skymap_proj("AZP") + + # Return + return + # Test CAR projection def test_skymap_car(self): """ @@ -207,13 +219,13 @@ class Test(GPythonTestSuite): # Return return - # Test TAN projection - def test_skymap_tan(self): + # Test MER projection + def test_skymap_mer(self): """ - Test TAN projection. + Test MER projection. """ # Execute generic test - self.test_skymap_proj("TAN") + self.test_skymap_proj("MER") # Return return @@ -229,13 +241,13 @@ class Test(GPythonTestSuite): # Return return - # Test AZP projection - def test_skymap_azp(self): + # Test TAN projection + def test_skymap_tan(self): """ - Test AZP projection. + Test TAN projection. """ # Execute generic test - self.test_skymap_proj("AZP") + self.test_skymap_proj("TAN") # Return return