Michael Mayer / gammalib

Commit ba12d62e92fb3f89438f50f445970d70ccbb0152

Authored by Jürgen Knödlseder
1 parent 445cf984

Add several methods to the HealPix class. 

The following methods have been added:
- set_z_phi
- loc2dir
- xyf2loc
- pix2xyf
- max_pixrad
- boundaries
- neighbours
Inline Side-by-side
Showing 3 changed files with 455 additions and 42 deletions
include/GHealpix.hpp
  View file @ba12d62
@@ -75,41 +75,49 @@ public: @@ -75,41 +75,49 @@ public:
75 virtual std::string print(const GChatter& chatter = NORMAL) const; 75 virtual std::string print(const GChatter& chatter = NORMAL) const;
76 76
77 // Other methods 77 // Other methods
78 - const int& npix(void) const;  
79 - const int& nside(void) const;  
80 - std::string ordering(void) const;  
81 - void ordering(const std::string& ordering); 78 + const int& npix(void) const;
  79 + const int& nside(void) const;
  80 + std::string ordering(void) const;
  81 + void ordering(const std::string& ordering);
  82 + std::vector<int> neighbours(const GSkyPixel& pixel) const;
  83 + std::vector<GSkyDir> boundaries(const GSkyPixel& pixel, const int& step = 1) const;
  84 + double max_pixrad(void) const;
82 85
83 private: 86 private:
84 // Private methods 87 // Private methods
85 - void init_members(void);  
86 - void copy_members(const GHealpix& wcs);  
87 - void free_members(void);  
88 - virtual bool compare(const GSkyProjection& proj) const; 88 + void init_members(void);
  89 + void copy_members(const GHealpix& wcs);
  90 + void free_members(void);
  91 + virtual bool compare(const GSkyProjection& proj) const;
89 92
90 // Low-level HealPix methods 93 // Low-level HealPix methods
91 - int compress_bits(const int& value) const;  
92 - int spread_bits(const int& value) const;  
93 - void nest2xyf(const int& pix, int* ix, int* iy, int* face) const;  
94 - void ring2xyf(const int& pix, int* ix, int* iy, int* face) const;  
95 - int xyf2nest(const int& ix, const int& iy, const int& face) const;  
96 - int xyf2ring(const int& ix, const int& iy, const int& face) const;  
97 - int nside2order(const int& nside) const;  
98 - int nest2ring(const int& pix) const;  
99 - int ring2nest(const int& pix) const;  
100 - void pix2xy(const int& ipix, int* x, int* y) const;  
101 - int xy2pix(int x, int y) const;  
102 - void pix2ang_ring(int ipix, double* theta, double* phi) const;  
103 - void pix2ang_nest(int ipix, double* theta, double* phi) const;  
104 - int ang2pix_z_phi_ring(double z, double phi) const;  
105 - int ang2pix_z_phi_nest(double z, double phi) const;  
106 - int ring_above(const double& z) const;  
107 - void get_ring_info(const int& ring, int* startpix, int* ringpix,  
108 - bool* shifted) const;  
109 - void get_ring_info(const int& ring, int* startpix, int* ringpix,  
110 - double* theta, bool* shifted) const;  
111 - GBilinear interpolator(const double& theta, const double& phi) const;  
112 - unsigned int isqrt(unsigned int arg) const; 94 + int compress_bits(const int& value) const;
  95 + int spread_bits(const int& value) const;
  96 + void pix2xyf(const int& pix, int* ix, int* iy, int* face) const;
  97 + void nest2xyf(const int& pix, int* ix, int* iy, int* face) const;
  98 + void ring2xyf(const int& pix, int* ix, int* iy, int* face) const;
  99 + int xyf2nest(const int& ix, const int& iy, const int& face) const;
  100 + int xyf2ring(const int& ix, const int& iy, const int& face) const;
  101 + void xyf2loc(const double& x, const double& y, const int& face,
  102 + double* z, double* phi) const;
  103 + int nside2order(const int& nside) const;
  104 + int nest2ring(const int& pix) const;
  105 + int ring2nest(const int& pix) const;
  106 + void pix2xy(const int& ipix, int* x, int* y) const;
  107 + int xy2pix(int x, int y) const;
  108 + void pix2ang_ring(int ipix, double* theta, double* phi) const;
  109 + void pix2ang_nest(int ipix, double* theta, double* phi) const;
  110 + int ang2pix_z_phi_ring(double z, double phi) const;
  111 + int ang2pix_z_phi_nest(double z, double phi) const;
  112 + int ring_above(const double& z) const;
  113 + void get_ring_info(const int& ring, int* startpix, int* ringpix,
  114 + bool* shifted) const;
  115 + void get_ring_info(const int& ring, int* startpix, int* ringpix,
  116 + double* theta, bool* shifted) const;
  117 + GBilinear interpolator(const double& theta, const double& phi) const;
  118 + unsigned int isqrt(unsigned int arg) const;
  119 + GVector set_z_phi(const double& z, const double& phi) const;
  120 + GSkyDir loc2dir(const double& z, const double& phi) const;
113 121
114 // Private data area 122 // Private data area
115 int m_nside; //!< Number of divisions of each base pixel (1-8192) 123 int m_nside; //!< Number of divisions of each base pixel (1-8192)
pyext/GHealpix.i
  View file @ba12d62
@@ -27,6 +27,13 @@ @@ -27,6 +27,13 @@
27 /* Put headers and other declarations here that are needed for compilation */ 27 /* Put headers and other declarations here that are needed for compilation */
28 #include "GHealpix.hpp" 28 #include "GHealpix.hpp"
29 %} 29 %}
  30 +%include "std_vector.i"
  31 +namespace std {
  32 + %template(SkyDirVector) vector<GSkyDir>;
  33 +}
  34 +namespace std {
  35 + %template(IntVector) vector<int>;
  36 +}
30 37
31 38
32 /***********************************************************************//** 39 /***********************************************************************//**
@@ -60,10 +67,13 @@ public: @@ -60,10 +67,13 @@ public:
60 virtual GBilinear interpolator(const GSkyDir& dir) const; 67 virtual GBilinear interpolator(const GSkyDir& dir) const;
61 68
62 // Other methods 69 // Other methods
63 - const int& npix(void) const;  
64 - const int& nside(void) const;  
65 - std::string ordering(void) const;  
66 - void ordering(const std::string& ordering); 70 + const int& npix(void) const;
  71 + const int& nside(void) const;
  72 + std::string ordering(void) const;
  73 + void ordering(const std::string& ordering);
  74 + std::vector<int> neighbours(const GSkyPixel& pixel) const;
  75 + std::vector<GSkyDir> boundaries(const GSkyPixel& pixel, const int& step = 1) const;
  76 + double max_pixrad(void) const;
67 }; 77 };
68 78
69 79
src/sky/GHealpix.cpp
  View file @ba12d62
@@ -32,6 +32,7 @@ @@ -32,6 +32,7 @@
32 #include "GException.hpp" 32 #include "GException.hpp"
33 #include "GTools.hpp" 33 #include "GTools.hpp"
34 #include "GMath.hpp" 34 #include "GMath.hpp"
  35 +#include "GVector.hpp"
35 #include "GHealpix.hpp" 36 #include "GHealpix.hpp"
36 37
37 /* __ Method name definitions ____________________________________________ */ 38 /* __ Method name definitions ____________________________________________ */
@@ -45,6 +46,8 @@ @@ -45,6 +46,8 @@
45 #define G_PIX2ANG_NEST "GHealpix::pix2ang_nest(int, double*, double*)" 46 #define G_PIX2ANG_NEST "GHealpix::pix2ang_nest(int, double*, double*)"
46 #define G_ORDERING_SET "GHealpix::ordering(std::string&)" 47 #define G_ORDERING_SET "GHealpix::ordering(std::string&)"
47 #define G_INTERPOLATOR "GHealpix::interpolator(double&, double&)" 48 #define G_INTERPOLATOR "GHealpix::interpolator(double&, double&)"
  49 +#define G_NEIGHBOURS "GHealpix::neighbours(GSkyPixel&)"
  50 +#define G_BOUNDARIES "GHealpix::boundaries(GSkyPixel&, int&)"
48 51
49 /* __ Macros _____________________________________________________________ */ 52 /* __ Macros _____________________________________________________________ */
50 53
@@ -55,10 +58,30 @@ @@ -55,10 +58,30 @@
55 /* __ Local prototypes ___________________________________________________ */ 58 /* __ Local prototypes ___________________________________________________ */
56 59
57 /* __ Constants __________________________________________________________ */ 60 /* __ Constants __________________________________________________________ */
58 -const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};  
59 -const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};  
60 -const int order_max = 13;  
61 -const int ns_max = 1 << order_max; 61 +const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
  62 +const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
  63 +const int nb_xoffset[] = {-1,-1, 0, 1, 1, 1, 0,-1};
  64 +const int nb_yoffset[] = { 0, 1, 1, 1, 0,-1,-1,-1};
  65 +const int nb_facearray[][12] = {{ 8, 9,10,11,-1,-1,-1,-1,10,11, 8, 9 }, // S
  66 + { 5, 6, 7, 4, 8, 9,10,11, 9,10,11, 8 }, // SE
  67 + { -1,-1,-1,-1, 5, 6, 7, 4,-1,-1,-1,-1 }, // E
  68 + { 4, 5, 6, 7,11, 8, 9,10,11, 8, 9,10 }, // SW
  69 + { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11 }, // center
  70 + { 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 }, // NE
  71 + { -1,-1,-1,-1, 7, 4, 5, 6,-1,-1,-1,-1 }, // W
  72 + { 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 }, // NW
  73 + { 2, 3, 0, 1,-1,-1,-1,-1, 0, 1, 2, 3 }}; // N
  74 +const int nb_swaparray[][3] = {{ 0,0,3 }, // S
  75 + { 0,0,6 }, // SE
  76 + { 0,0,0 }, // E
  77 + { 0,0,5 }, // SW
  78 + { 0,0,0 }, // center
  79 + { 5,0,0 }, // NE
  80 + { 0,0,0 }, // W
  81 + { 6,0,0 }, // NW
  82 + { 3,0,0 }}; // N
  83 +const int order_max = 13;
  84 +const int ns_max = 1 << order_max;
62 85
63 /* __ Static conversion arrays ___________________________________________ */ 86 /* __ Static conversion arrays ___________________________________________ */
64 static short ctab[0x100]; 87 static short ctab[0x100];
@@ -534,6 +557,239 @@ void GHealpix::ordering(const std::string&amp; ordering) @@ -534,6 +557,239 @@ void GHealpix::ordering(const std::string&amp; ordering)
534 557
535 558
536 /***********************************************************************//** 559 /***********************************************************************//**
  560 + * @brief Return neighbouring pixels of a pixel
  561 + *
  562 + * @param[in] pix Pixel index.
  563 + * @return Array of neighbouring pixels.
  564 + *
  565 + * Returns the 8 neighbours of a given pixel. The method returns a vector
  566 + * with contains the pixel indices of the SW, W, NW, N, NE, E, SE and S
  567 + * neighbours of @p pix (in the given order). If a neighbour does not exist
  568 + * (this can only be the case for the W, N, E and S neighbors), its entry is
  569 + * set to -1.
  570 + *
  571 + * This method has been adapted from the neighbors() function located in the
  572 + * file healpix_base.cc in Healpix version 3.20.
  573 + ***************************************************************************/
  574 +std::vector<int> GHealpix::neighbours(const GSkyPixel& pixel) const
  575 +{
  576 + // Throw an exception if sky map pixel is not 1D
  577 + if (!pixel.is_1D()) {
  578 + std::string msg = "Sky map pixel "+pixel.print()+" is not"
  579 + " 1-dimensional.\n"
  580 + "Only 1-dimensional pixels are supported by the"
  581 + " Healpix projection.";
  582 + throw GException::invalid_argument(G_NEIGHBOURS, msg);
  583 + }
  584 +
  585 + // Initialise result array
  586 + std::vector<int> result;
  587 +
  588 + // Determine pixel index and face number
  589 + int ix;
  590 + int iy;
  591 + int face_num;
  592 + pix2xyf(int(pixel), &ix, &iy, &face_num);
  593 +
  594 + // ...
  595 + const int nsm1 = m_nside - 1;
  596 +
  597 + // ...
  598 + if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) {
  599 +
  600 + // Ring ordering scheme ...
  601 + if (m_ordering == 0) { // Ring?
  602 + for (int m = 0; m < 8; ++m) {
  603 + int index = xyf2ring(ix+nb_xoffset[m],iy+nb_yoffset[m],face_num);
  604 + result.push_back(index);
  605 + }
  606 + }
  607 +
  608 + // ... or nested scheme?
  609 + else {
  610 + int fpix = int(face_num) << (2*m_order);
  611 + int px0 = spread_bits(ix);
  612 + int py0 = spread_bits(iy) << 1;
  613 + int pxp = spread_bits(ix+1);
  614 + int pyp = spread_bits(iy+1) << 1;
  615 + int pxm = spread_bits(ix-1);
  616 + int pym = spread_bits(iy-1) << 1;
  617 + result.push_back(fpix + pxm + py0);
  618 + result.push_back(fpix + pxm + pyp);
  619 + result.push_back(fpix + px0 + pyp);
  620 + result.push_back(fpix + pxp + pyp);
  621 + result.push_back(fpix + pxp + py0);
  622 + result.push_back(fpix + pxp + pym);
  623 + result.push_back(fpix + px0 + pym);
  624 + result.push_back(fpix + pxm + pym);
  625 + }
  626 +
  627 + } // endif
  628 +
  629 + // ...
  630 + else {
  631 + for (int i = 0; i < 8; ++i) {
  632 + int x = ix + nb_xoffset[i];
  633 + int y = iy + nb_yoffset[i];
  634 + int nbnum = 4;
  635 + if (x < 0) {
  636 + x += m_nside;
  637 + nbnum -= 1;
  638 + }
  639 + else if (x >= m_nside) {
  640 + x -= m_nside;
  641 + nbnum += 1;
  642 + }
  643 + if (y < 0) {
  644 + y += m_nside;
  645 + nbnum -= 3;
  646 + }
  647 + else if (y >= m_nside) {
  648 + y -= m_nside;
  649 + nbnum += 3;
  650 + }
  651 +
  652 + // Compute face
  653 + int f = nb_facearray[nbnum][face_num];
  654 +
  655 + // If face is valid then compute index
  656 + if (f >= 0) {
  657 + int bits = nb_swaparray[nbnum][face_num>>2];
  658 + if (bits & 1) {
  659 + x = m_nside - x - 1;
  660 + }
  661 + if (bits & 2) {
  662 + y = m_nside - y - 1;
  663 + }
  664 + if (bits & 4) {
  665 + std::swap(x,y);
  666 + }
  667 + if (m_ordering == 0) { // Ring?
  668 + result.push_back(xyf2ring(x, y, f));
  669 + }
  670 + else {
  671 + result.push_back(xyf2nest(x, y, f));
  672 + }
  673 + }
  674 +
  675 + // ... otherwise push back an invalid pixel
  676 + else {
  677 + result.push_back(-1);
  678 + }
  679 +
  680 + } // endfor
  681 +
  682 + } // endelse
  683 +
  684 + // Return indices
  685 + return result;
  686 +}
  687 +
  688 +
  689 +/***********************************************************************//**
  690 + * @brief Return pixel boundaries
  691 + *
  692 + * @param[in] pixel Sky pixel.
  693 + * @param[in] step Number of returned points (4*step, defaults to 4).
  694 + * @return Array of pixel boundaries.
  695 + *
  696 + * The method returns a vector with sky directions along the boundary of a
  697 + * HealPix pixel. By default, the 4 corners of HealPix pixel will be
  698 + * returned. The first point corresponds to the northernmost corner, the
  699 + * subsequent points follow the pixel boundary through west, south and east
  700 + * corners. If step>1 more intermediate points will be added after the 4th
  701 + * corner of the pixel.
  702 + *
  703 + * This method has been adapted from the boundaries() function located in the
  704 + * file healpix_base.cc in Healpix version 3.20.
  705 + ***************************************************************************/
  706 +std::vector<GSkyDir> GHealpix::boundaries(const GSkyPixel& pixel,
  707 + const int& step) const
  708 +{
  709 + // Throw an exception if sky map pixel is not 1D
  710 + if (!pixel.is_1D()) {
  711 + std::string msg = "Sky map pixel "+pixel.print()+" is not"
  712 + " 1-dimensional.\n"
  713 + "Only 1-dimensional pixels are supported by the"
  714 + " Healpix projection.";
  715 + throw GException::invalid_argument(G_BOUNDARIES, msg);
  716 + }
  717 +
  718 + // Allocate boundaries
  719 + std::vector<GSkyDir> boundaries;
  720 +
  721 + // Determine pixel index and face number
  722 + int ix;
  723 + int iy;
  724 + int face;
  725 + pix2xyf(int(pixel), &ix, &iy, &face);
  726 +
  727 + // ...
  728 + double dc = 0.5 / m_nside;
  729 + double xc = (ix + 0.5) / m_nside;
  730 + double yc = (iy + 0.5) / m_nside;
  731 + double d = 1.0 / (step*m_nside);
  732 +
  733 + //
  734 + for (int i = 0; i < step; ++i) {
  735 +
  736 + // Declare local coordinates
  737 + double z;
  738 + double phi;
  739 +
  740 + // First coordinate
  741 + xyf2loc(xc+dc-i*d, yc+dc, face, &z, &phi);
  742 + boundaries.push_back(loc2dir(z, phi));
  743 +
  744 + // Second coordinate
  745 + xyf2loc(xc-dc, yc+dc-i*d, face, &z, &phi);
  746 + boundaries.push_back(loc2dir(z, phi));
  747 +
  748 + // Third coordinate
  749 + xyf2loc(xc-dc+i*d, yc-dc, face, &z, &phi);
  750 + boundaries.push_back(loc2dir(z, phi));
  751 +
  752 + // Forth coordinate
  753 + xyf2loc(xc+dc, yc-dc+i*d, face, &z, &phi);
  754 + boundaries.push_back(loc2dir(z, phi));
  755 +
  756 + } // endfor: looped over step size
  757 +
  758 + // Return boundaries
  759 + return boundaries;
  760 +}
  761 +
  762 +
  763 +/***********************************************************************//**
  764 + * @brief Return maximum angular distance between pixel centre and corners
  765 + *
  766 + * @return Maximum angular distance between pixel centre and corners (radians).
  767 + *
  768 + * Returns the maximum angular distance (in radians) between any pixel
  769 + * centre and its corners.
  770 + *
  771 + * This method has been adapted from the max_pixrad() function located in the
  772 + * file healpix_base.cc in Healpix version 3.20.
  773 + ***************************************************************************/
  774 +double GHealpix::max_pixrad(void) const
  775 +{
  776 + // Compute ...
  777 + double t1 = 1.0 - 1.0/m_nside;
  778 + t1 *= t1;
  779 +
  780 + // Get vectors
  781 + GVector va = set_z_phi(2.0/3.0, gammalib::pi/(4*m_nside));
  782 + GVector vb = set_z_phi(1.0 - t1/3, 0.0);
  783 +
  784 + // Get angle
  785 + double angle = std::atan2(norm(cross(va, vb)), va * vb);
  786 +
  787 + // Return angle
  788 + return angle;
  789 +}
  790 +
  791 +
  792 +/***********************************************************************//**
537 * @brief Print WCS information 793 * @brief Print WCS information
538 * 794 *
539 * @param[in] chatter Chattiness (defaults to NORMAL). 795 * @param[in] chatter Chattiness (defaults to NORMAL).
@@ -720,9 +976,37 @@ int GHealpix::spread_bits(const int&amp; value) const @@ -720,9 +976,37 @@ int GHealpix::spread_bits(const int&amp; value) const
720 976
721 977
722 /***********************************************************************//** 978 /***********************************************************************//**
  979 + * @brief Convert pixel number in to (x,y,face) tuple
  980 + *
  981 + * @param[in] pix Pixel number.
  982 + * @param[out] ix X index.
  983 + * @param[out] iy Y index.
  984 + * @param[out] face Face number.
  985 + *
  986 + * This method has been adapted from the pix2xyf() function located in
  987 + * the file healpix_base.cc in Healpix version 3.20.
  988 + ***************************************************************************/
  989 +void GHealpix::pix2xyf(const int& pix, int* ix, int* iy, int* face) const
  990 +{
  991 + // Handle ring scheme ...
  992 + if (m_ordering == 0) {
  993 + ring2xyf(pix, ix, iy, face);
  994 + }
  995 +
  996 + // ... or nested scheme
  997 + else {
  998 + nest2xyf(pix, ix, iy, face);
  999 + }
  1000 +
  1001 + // Return
  1002 + return;
  1003 +}
  1004 +
  1005 +
  1006 +/***********************************************************************//**
723 * @brief Convert pixel number in nested scheme to (x,y,face) tuple 1007 * @brief Convert pixel number in nested scheme to (x,y,face) tuple
724 * 1008 *
725 - * @param[in] pix Pixel number in nested scheme 1009 + * @param[in] pix Pixel number in nested scheme.
726 * @param[out] ix X index. 1010 * @param[out] ix X index.
727 * @param[out] iy Y index. 1011 * @param[out] iy Y index.
728 * @param[out] face Face number. 1012 * @param[out] face Face number.
@@ -887,6 +1171,91 @@ int GHealpix::xyf2ring(const int&amp; ix, const int&amp; iy, const int&amp; face) const @@ -887,6 +1171,91 @@ int GHealpix::xyf2ring(const int&amp; ix, const int&amp; iy, const int&amp; face) const
887 1171
888 1172
889 /***********************************************************************//** 1173 /***********************************************************************//**
  1174 + * @brief Convert (x,y,f) tuple into local coordinates
  1175 + *
  1176 + * @param[in] x X value.
  1177 + * @param[in] y Y value.
  1178 + * @param[in] face Face number.
  1179 + * @param[out] z Z value.
  1180 + * @param[out] phi Phi value.
  1181 + *
  1182 + * This method has been adapted from the xyf2loc() function located in the
  1183 + * file healpix_base.cc in Healpix version 3.20.
  1184 + ***************************************************************************/
  1185 +void GHealpix::xyf2loc(const double& x, const double& y, const int& face,
  1186 + double* z, double* phi) const
  1187 +{
  1188 + // ...
  1189 + double jr = jrll[face] - x - y;
  1190 + double nr;
  1191 +
  1192 + // Compute z
  1193 + if (jr < 1) {
  1194 + nr = jr;
  1195 + double tmp = nr*nr / 3.0;
  1196 + *z = 1 - tmp;
  1197 + }
  1198 + else if (jr > 3) {
  1199 + nr = 4 - jr;
  1200 + double tmp = nr*nr / 3.0;
  1201 + *z = tmp - 1;
  1202 + }
  1203 + else {
  1204 + nr = 1;
  1205 + *z = (2.0-jr) * 2.0/3.0;
  1206 + }
  1207 +
  1208 + // Compute Phi
  1209 + double tmp = jpll[face] * nr + x - y;
  1210 + if (tmp < 0) {
  1211 + tmp += 8;
  1212 + }
  1213 + if (tmp >= 8) {
  1214 + tmp -= 8;
  1215 + }
  1216 + *phi = (nr < 1.0e-15) ? 0.0 : (0.5 * gammalib::pihalf * tmp) / nr;
  1217 +
  1218 + // Return
  1219 + return;
  1220 +}
  1221 +
  1222 +
  1223 +/***********************************************************************//**
  1224 + * @brief Convert local coordinate into sky direction
  1225 + *
  1226 + * @param[in] z Z value.
  1227 + * @param[in] phi Phi value.
  1228 + * @return Sky direction.
  1229 + *
  1230 + * This method has been adapted from the locToVec3() function located in
  1231 + * the file healpix_base.cc in Healpix version 3.20.
  1232 + ***************************************************************************/
  1233 +GSkyDir GHealpix::loc2dir(const double& z, const double& phi) const
  1234 +{
  1235 + // Compute longitude and latitude
  1236 + double sintheta = std::sqrt((1.0 - z) * (1.0 + z));
  1237 + double longitude = std::atan2(sintheta * std::sin(phi), sintheta * std::cos(phi));
  1238 + double latitude = std::asin(z);
  1239 +
  1240 + // Set sky direction
  1241 + GSkyDir dir;
  1242 + switch (m_coordsys) {
  1243 + case 0:
  1244 + dir.radec(longitude, latitude);
  1245 + break;
  1246 + case 1:
  1247 + dir.lb(longitude, latitude);
  1248 + break;
  1249 + default:
  1250 + break;
  1251 + }
  1252 +
  1253 + // Return sky direction
  1254 + return dir;
  1255 +}
  1256 +
  1257 +
  1258 +/***********************************************************************//**
890 * @brief Convert nside to order 1259 * @brief Convert nside to order
891 * 1260 *
892 * @param[in] nside Number of sides. 1261 * @param[in] nside Number of sides.
@@ -1459,7 +1828,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const @@ -1459,7 +1828,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const
1459 // Now handle the special case that the colatitude is North of all 1828 // Now handle the special case that the colatitude is North of all
1460 // rings 1829 // rings
1461 if (ir1 == 0) { 1830 if (ir1 == 0) {
1462 - double wtheta = theta/theta2; 1831 + double wtheta = theta/theta2;
1463 interpolator.weight3() *= wtheta; 1832 interpolator.weight3() *= wtheta;
1464 interpolator.weight4() *= wtheta; 1833 interpolator.weight4() *= wtheta;
1465 double fac = (1.0-wtheta)*0.25; 1834 double fac = (1.0-wtheta)*0.25;
@@ -1473,7 +1842,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const @@ -1473,7 +1842,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const
1473 1842
1474 // ... and now the case that the colatitude is South of all rings 1843 // ... and now the case that the colatitude is South of all rings
1475 else if (ir2 == 4*m_nside) { 1844 else if (ir2 == 4*m_nside) {
1476 - double wtheta = (theta-theta1) / (gammalib::pi-theta1); 1845 + double wtheta = (theta-theta1) / (gammalib::pi-theta1);
1477 interpolator.weight1() *= (1.0 - wtheta); 1846 interpolator.weight1() *= (1.0 - wtheta);
1478 interpolator.weight2() *= (1.0 - wtheta); 1847 interpolator.weight2() *= (1.0 - wtheta);
1479 double fac = wtheta*0.25; 1848 double fac = wtheta*0.25;
@@ -1487,7 +1856,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const @@ -1487,7 +1856,7 @@ GBilinear GHealpix::interpolator(const double&amp; theta, const double&amp; phi) const
1487 1856
1488 // ... and now multiply-in the theta weights for the general case 1857 // ... and now multiply-in the theta weights for the general case
1489 else { 1858 else {
1490 - double wtheta = (theta-theta1) / (theta2-theta1); 1859 + double wtheta = (theta-theta1) / (theta2-theta1);
1491 interpolator.weight1() *= (1.0 - wtheta); 1860 interpolator.weight1() *= (1.0 - wtheta);
1492 interpolator.weight2() *= (1.0 - wtheta); 1861 interpolator.weight2() *= (1.0 - wtheta);
1493 interpolator.weight3() *= wtheta; 1862 interpolator.weight3() *= wtheta;
@@ -1520,3 +1889,29 @@ unsigned int GHealpix::isqrt(unsigned int arg) const @@ -1520,3 +1889,29 @@ unsigned int GHealpix::isqrt(unsigned int arg) const
1520 // Return 1889 // Return
1521 return unsigned(std::sqrt(arg+0.5)); 1890 return unsigned(std::sqrt(arg+0.5));
1522 } 1891 }
  1892 +
  1893 +
  1894 +/***********************************************************************//**
  1895 + * @brief Return 3D vector
  1896 + *
  1897 + * @param[in] z Z value.
  1898 + * @param[in] phi Phi value.
  1899 + * @return 3D vector
  1900 + *
  1901 + * This method has been adapted from the set_z_phi() function located in the
  1902 + * file vec3.h in Healpix version 3.20.
  1903 + ***************************************************************************/
  1904 +GVector GHealpix::set_z_phi(const double& z, const double& phi) const
  1905 +{
  1906 + // Initialise 3D vector
  1907 + GVector vector(3);
  1908 +
  1909 + // Assign elements
  1910 + double sintheta = std::sqrt((1.0 - z) * (1.0 + z));
  1911 + vector[0] = sintheta * std::cos(phi);
  1912 + vector[1] = sintheta * std::sin(phi);
  1913 + vector[2] = z;
  1914 +
  1915 + // Return vector
  1916 + return vector;
  1917 +}