GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GHealpix.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GHealpix.cpp - Healpix projection class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2015 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GHealpix.cpp
23  * @brief HealPix projection class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GVector.hpp"
36 #include "GHealpix.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_CONSTRUCT "GHealpix::GHealpix(int& ,std::string& ,std::string&)"
40 #define G_READ "GHealpix::read(GFitsHDU&)"
41 #define G_XY2DIR "GHealpix::xy2dir(GSkyPixel&)"
42 #define G_DIR2XY2 "GHealpix::dir2xy(GSkyDir&)"
43 #define G_NEST2RING "GHealpix::nest2ring(int&)"
44 #define G_RING2NEST "GHealpix::ring2nest(int&)"
45 #define G_PIX2ANG_RING "GHealpix::pix2ang_ring(int, double*, double*)"
46 #define G_PIX2ANG_NEST "GHealpix::pix2ang_nest(int, double*, double*)"
47 #define G_ORDERING_SET "GHealpix::ordering(std::string&)"
48 #define G_INTERPOLATOR "GHealpix::interpolator(double&, double&)"
49 #define G_NEIGHBOURS "GHealpix::neighbours(GSkyPixel&)"
50 #define G_BOUNDARIES "GHealpix::boundaries(GSkyPixel&, int&)"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 /* __ Local prototypes ___________________________________________________ */
59 
60 /* __ Constants __________________________________________________________ */
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;
85 
86 /* __ Static conversion arrays ___________________________________________ */
87 static short ctab[0x100];
88 static short utab[0x100];
89 
90 /* __ Globals ____________________________________________________________ */
91 
92 
93 /*==========================================================================
94  = =
95  = Constructors/destructors =
96  = =
97  ==========================================================================*/
98 
99 /***********************************************************************//**
100  * @brief Void constructor
101  ***************************************************************************/
103 {
104  // Initialise class members
105  init_members();
106 
107  // Return
108  return;
109 }
110 
111 
112 /***********************************************************************//**
113  * @brief Constructor
114  *
115  * @param[in] nside Number of sides.
116  * @param[in] order Pixel ordering ('RING' or 'NESTED').
117  * @param[in] coords Coordinate system ('EQU' or 'GAL').
118  *
119  * @exception GException::wcs_hpx_bad_nside
120  * Invalid nside parameter.
121  * @exception GException::wcs_bad_coords
122  * Invalid coordsys parameter.
123  * @exception GException::wcs_hpx_bad_ordering
124  * Invalid ordering parameter.
125  ***************************************************************************/
126 GHealpix::GHealpix(const int& nside,
127  const std::string& order,
128  const std::string& coords) : GSkyProjection()
129 {
130  // Initialise class members
131  init_members();
132 
133  // Check nside parameter (power of 2 between 1 and 8192)
134  if (nside2order(nside) == -1) {
136  }
137 
138  // Set coordinate system
139  coordsys(coords);
140 
141  // Set pixel ordering
142  ordering(order);
143 
144  // Set Healpix parameters
145  m_nside = nside;
147  m_ncap = 2 * (m_npface - m_nside);
148  m_num_pixels = 12 * m_npface;
149  m_fact2 = 4.0 / m_num_pixels;
150  m_fact1 = 2 * m_nside * m_fact2;
152  m_order = nside2order(m_nside);
153 
154  // Return
155  return;
156 }
157 
158 
159 /***********************************************************************//**
160  * @brief Constructor from FITS HDU table
161  *
162  * @param[in] hdu FITS HDU.
163  ***************************************************************************/
165 {
166  // Initialise class members
167  init_members();
168 
169  // Read Healpix definition from FITS table
170  read(hdu);
171 
172  // Return
173  return;
174 }
175 
176 
177 /***********************************************************************//**
178  * @brief Copy constructor
179  *
180  * @param[in] proj Healpix projection.
181  ***************************************************************************/
183 {
184  // Initialise class members for clean destruction
185  init_members();
186 
187  // Copy members
188  copy_members(proj);
189 
190  // Return
191  return;
192 }
193 
194 
195 /***********************************************************************//**
196  * @brief Destructor
197  ***************************************************************************/
199 {
200  // Free members
201  free_members();
202 
203  // Return
204  return;
205 }
206 
207 
208 /*==========================================================================
209  = =
210  = Operators =
211  = =
212  ==========================================================================*/
213 
214 /***********************************************************************//**
215  * @brief Assignment operator
216  *
217  * @param[in] proj Healpix projection.
218  * @return Healpix projection.
219  ***************************************************************************/
221 {
222  // Execute only if object is not identical
223  if (this != &proj) {
224 
225  // Copy base class members
226  this->GSkyProjection::operator=(proj);
227 
228  // Free members
229  free_members();
230 
231  // Initialise private members for clean destruction
232  init_members();
233 
234  // Copy members
235  copy_members(proj);
236 
237  } // endif: object was not identical
238 
239  // Return this object
240  return *this;
241 }
242 
243 
244 /*==========================================================================
245  = =
246  = Public methods =
247  = =
248  ==========================================================================*/
249 
250 /***********************************************************************//**
251  * @brief Clear object
252  *
253  * This method properly resets the object to an initial state.
254  ***************************************************************************/
255 void GHealpix::clear(void)
256 {
257  // Free class members (base and derived classes, derived class first)
258  free_members();
260 
261  // Initialise members
263  init_members();
264 
265  // Return
266  return;
267 }
268 
269 
270 /***********************************************************************//**
271  * @brief Clone instance
272  ***************************************************************************/
274 {
275  return new GHealpix(*this);
276 }
277 
278 
279 /***********************************************************************//**
280  * @brief Read Healpix definition from FITS header
281  *
282  * @param[in] hdu FITS HDU.
283  *
284  * @exception GException::wcs
285  * Unable to load Healpix definition from HDU.
286  * @exception GException::wcs_bad_coords
287  * Invalid coordsys parameter.
288  * @exception GException::wcs_hpx_bad_ordering
289  * Invalid ordering parameter.
290  ***************************************************************************/
291 void GHealpix::read(const GFitsHDU& hdu)
292 {
293  // Clear object
294  clear();
295 
296  // Check if we have a healpix representation
297  if (hdu.string("PIXTYPE") != "HEALPIX") {
298  throw GException::wcs(G_READ, "HDU does not contain Healpix data");
299  }
300 
301  // Get pixel ordering. First search for the ORDERING keyword, then
302  // search for the ORDER keyword
303  std::string ordering;
304  if (hdu.has_card("ORDERING")) {
305  ordering = hdu.string("ORDERING");
306  }
307  else if (hdu.has_card("ORDER")) {
308  ordering = hdu.string("ORDER");
309  }
310 
311  // Get coordinate system. First search for HIER_CRD keyword (this has been
312  // used in older versions of LAT exposure cubes). If not found then search
313  // for standard COORDSYS keyword.
314  std::string coordsys;
315  if (hdu.has_card("HIER_CRD")) {
316  coordsys = hdu.string("HIER_CRD");
317  }
318  else if (hdu.has_card("COORDSYS")) {
319  coordsys = hdu.string("COORDSYS");
320  }
321 
322  // Set coordinate system
323  this->coordsys(coordsys);
324 
325  // Set pixel ordering
326  this->ordering(ordering);
327 
328  // Get Healpix resolution and determine number of pixels and solid angle
329  m_nside = hdu.integer("NSIDE");
331  m_ncap = 2 * (m_npface - m_nside);
332  m_num_pixels = 12 * m_npface;
333  m_fact2 = 4.0 / m_num_pixels;
334  m_fact1 = 2 * m_nside * m_fact2;
336  m_order = nside2order(m_nside);
337 
338  // Return
339  return;
340 }
341 
342 
343 /***********************************************************************//**
344  * @brief Write Healpix definition into FITS HDU
345  *
346  * @param[in] hdu FITS HDU.
347  *
348  * Writes the following keywords in the FITS HDU:
349  * EXTNAME = HEALPIX
350  * PIXTYPE = HEALPIX
351  * NSIDE = nside() = m_nside
352  * FIRSTPIX = 0
353  * LASTPIX = npix()-1 = m_num_pixels-1
354  * ORDERING = ordering()
355  * COORDSYS = coordsys()
356  ***************************************************************************/
357 void GHealpix::write(GFitsHDU& hdu) const
358 {
359  // Set extension name
360  hdu.extname("HEALPIX");
361 
362  // Set pixtype (kludge to avoid that a boolean is written)
363  std::string pixtype = "HEALPIX";
364 
365  // Set keywords
366  hdu.card("PIXTYPE", pixtype, "HEALPix pixelisation");
367  hdu.card("NSIDE", nside(), "HEALPix resolution parameter");
368  hdu.card("FIRSTPIX", 0, "Index of first pixel");
369  hdu.card("LASTPIX", npix()-1, "Index of last pixel");
370  hdu.card("ORDERING", ordering(),
371  "Pixel ordering scheme, either RING or NESTED");
372  hdu.card("COORDSYS", coordsys(),
373  "Coordinate system, either EQU or GAL");
374 
375  // Return
376  return;
377 }
378 
379 
380 /***********************************************************************//**
381  * @brief Returns sky direction of pixel
382  *
383  * @param[in] pixel Sky map pixel.
384  ***************************************************************************/
386 {
387  // Throw an exception if sky map pixel is not 1D
388  if (!pixel.is_1D()) {
389  std::string msg = "Sky map pixel "+pixel.print()+" is not"
390  " 1-dimensional.\n"
391  "Only 1-dimensional pixels are supported by the"
392  " Healpix projection.";
394  }
395 
396  // Perform ordering dependent conversion
397  double theta = 0.0;
398  double phi = 0.0;
399  switch (m_ordering) {
400  case 0:
401  pix2ang_ring(int(pixel), &theta, &phi);
402  break;
403  case 1:
404  pix2ang_nest(int(pixel), &theta, &phi);
405  break;
406  default:
407  break;
408  }
409 
410  // Store coordinate system-dependent result
411  GSkyDir result;
412  switch (m_coordsys) {
413  case 0:
414  result.radec(phi, gammalib::pihalf - theta);
415  break;
416  case 1:
417  result.lb(phi, gammalib::pihalf - theta);
418  break;
419  default:
420  break;
421  }
422 
423  // Return
424  return result;
425 }
426 
427 
428 /***********************************************************************//**
429  * @brief Returns sky map pixel of sky coordinate
430  *
431  * @param[in] dir Sky coordinate.
432  * @return Sky map pixel.
433  ***************************************************************************/
435 {
436  // Compute coordinate system dependent (z,phi)
437  double z = 0;
438  double phi = 0;
439  switch (m_coordsys) {
440  case 0:
441  z = cos(gammalib::pihalf - dir.dec());
442  phi = dir.ra();
443  break;
444  case 1:
445  z = cos(gammalib::pihalf - dir.b());
446  phi = dir.l();
447  break;
448  default:
449  break;
450  }
451 
452  // Perform ordering dependent conversion
453  int index;
454  switch (m_ordering) {
455  case 0:
456  index = ang2pix_z_phi_ring(z, phi);
457  break;
458  case 1:
459  index = ang2pix_z_phi_nest(z, phi);
460  break;
461  default:
462  break;
463  }
464 
465  // Return sky map pixel
466  return (GSkyPixel(index));
467 }
468 
469 
470 /***********************************************************************//**
471  * @brief Return interpolator for given sky direction
472  *
473  * @param[in] dir Sky direction
474  ***************************************************************************/
476 {
477  // Compute coordinate system dependent theta and phi (in radians)
478  double theta = 0;
479  double phi = 0;
480  switch (m_coordsys) {
481  case 0:
482  theta = gammalib::pihalf - dir.dec();
483  phi = dir.ra();
484  break;
485  case 1:
486  theta = gammalib::pihalf - dir.b();
487  phi = dir.l();
488  break;
489  default:
490  break;
491  }
492 
493  // Normalize theta and phi value
494  theta = gammalib::modulo(theta, gammalib::twopi);
495  if (theta > gammalib::pi) {
496  phi += gammalib::pi;
497  theta = gammalib::twopi - theta;
498  }
499  phi = gammalib::modulo(phi, gammalib::twopi);
500 
501  // Get interpolator
502  GBilinear interpolator = this->interpolator(theta, phi);
503 
504  // Return interpolator
505  return (interpolator);
506 }
507 
508 
509 /***********************************************************************//**
510  * @brief Returns ordering parameter.
511  ***************************************************************************/
512 std::string GHealpix::ordering(void) const
513 {
514  // Set pixel ordering type
515  std::string ordering;
516  switch (m_ordering) {
517  case 0:
518  ordering = "RING";
519  break;
520  case 1:
521  ordering = "NESTED";
522  break;
523  default:
524  ordering = "UNKNOWN";
525  break;
526  }
527 
528  // Return ordering
529  return ordering;
530 }
531 
532 
533 /***********************************************************************//**
534  * @brief Set pixel ordering.
535  *
536  * @param[in] ordering Pixel ordering (RING or NEST/NESTED).
537  *
538  * @exception GException::wcs_hpx_bad_ordering
539  * Invalid ordering parameter.
540  ***************************************************************************/
541 void GHealpix::ordering(const std::string& ordering)
542 {
543  // Convert argument to upper case
544  std::string uordering = gammalib::toupper(ordering);
545 
546  // Set pixel ordering
547  if (uordering == "RING") {
548  m_ordering = 0;
549  }
550  else if (uordering == "NESTED" || uordering == "NEST") {
551  m_ordering = 1;
552  }
553  else {
555  }
556 
557  // Return
558  return;
559 }
560 
561 
562 /***********************************************************************//**
563  * @brief Return neighbouring pixels of a pixel
564  *
565  * @param[in] pixel Pixel index.
566  * @return Array of neighbouring pixels.
567  *
568  * Returns the 8 neighbours of a given pixel. The method returns a vector
569  * with contains the pixel indices of the SW, W, NW, N, NE, E, SE and S
570  * neighbours of @p pix (in the given order). If a neighbour does not exist
571  * (this can only be the case for the W, N, E and S neighbors), its entry is
572  * set to -1.
573  *
574  * This method has been adapted from the neighbors() function located in the
575  * file healpix_base.cc in Healpix version 3.20.
576  ***************************************************************************/
577 std::vector<int> GHealpix::neighbours(const GSkyPixel& pixel) const
578 {
579  // Throw an exception if sky map pixel is not 1D
580  if (!pixel.is_1D()) {
581  std::string msg = "Sky map pixel "+pixel.print()+" is not"
582  " 1-dimensional.\n"
583  "Only 1-dimensional pixels are supported by the"
584  " Healpix projection.";
586  }
587 
588  // Initialise result array
589  std::vector<int> result;
590 
591  // Determine pixel index and face number
592  int ix;
593  int iy;
594  int face_num;
595  pix2xyf(int(pixel), &ix, &iy, &face_num);
596 
597  // ...
598  const int nsm1 = m_nside - 1;
599 
600  // ...
601  if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) {
602 
603  // Ring ordering scheme ...
604  if (m_ordering == 0) { // Ring?
605  for (int m = 0; m < 8; ++m) {
606  int index = xyf2ring(ix+nb_xoffset[m],iy+nb_yoffset[m],face_num);
607  result.push_back(index);
608  }
609  }
610 
611  // ... or nested scheme?
612  else {
613  int fpix = int(face_num) << (2*m_order);
614  int px0 = spread_bits(ix);
615  int py0 = spread_bits(iy) << 1;
616  int pxp = spread_bits(ix+1);
617  int pyp = spread_bits(iy+1) << 1;
618  int pxm = spread_bits(ix-1);
619  int pym = spread_bits(iy-1) << 1;
620  result.push_back(fpix + pxm + py0);
621  result.push_back(fpix + pxm + pyp);
622  result.push_back(fpix + px0 + pyp);
623  result.push_back(fpix + pxp + pyp);
624  result.push_back(fpix + pxp + py0);
625  result.push_back(fpix + pxp + pym);
626  result.push_back(fpix + px0 + pym);
627  result.push_back(fpix + pxm + pym);
628  }
629 
630  } // endif
631 
632  // ...
633  else {
634  for (int i = 0; i < 8; ++i) {
635  int x = ix + nb_xoffset[i];
636  int y = iy + nb_yoffset[i];
637  int nbnum = 4;
638  if (x < 0) {
639  x += m_nside;
640  nbnum -= 1;
641  }
642  else if (x >= m_nside) {
643  x -= m_nside;
644  nbnum += 1;
645  }
646  if (y < 0) {
647  y += m_nside;
648  nbnum -= 3;
649  }
650  else if (y >= m_nside) {
651  y -= m_nside;
652  nbnum += 3;
653  }
654 
655  // Compute face
656  int f = nb_facearray[nbnum][face_num];
657 
658  // If face is valid then compute index
659  if (f >= 0) {
660  int bits = nb_swaparray[nbnum][face_num>>2];
661  if (bits & 1) {
662  x = m_nside - x - 1;
663  }
664  if (bits & 2) {
665  y = m_nside - y - 1;
666  }
667  if (bits & 4) {
668  std::swap(x,y);
669  }
670  if (m_ordering == 0) { // Ring?
671  result.push_back(xyf2ring(x, y, f));
672  }
673  else {
674  result.push_back(xyf2nest(x, y, f));
675  }
676  }
677 
678  // ... otherwise push back an invalid pixel
679  else {
680  result.push_back(-1);
681  }
682 
683  } // endfor
684 
685  } // endelse
686 
687  // Return indices
688  return result;
689 }
690 
691 
692 /***********************************************************************//**
693  * @brief Return pixel boundaries
694  *
695  * @param[in] pixel Sky pixel.
696  * @param[in] step Number of returned points (4*step, defaults to 4).
697  * @return Array of pixel boundaries.
698  *
699  * The method returns a vector with sky directions along the boundary of a
700  * HealPix pixel. By default, the 4 corners of HealPix pixel will be
701  * returned. The first point corresponds to the northernmost corner, the
702  * subsequent points follow the pixel boundary through west, south and east
703  * corners. If step>1 more intermediate points will be added after the 4th
704  * corner of the pixel.
705  *
706  * This method has been adapted from the boundaries() function located in the
707  * file healpix_base.cc in Healpix version 3.20.
708  ***************************************************************************/
709 std::vector<GSkyDir> GHealpix::boundaries(const GSkyPixel& pixel,
710  const int& step) const
711 {
712  // Throw an exception if sky map pixel is not 1D
713  if (!pixel.is_1D()) {
714  std::string msg = "Sky map pixel "+pixel.print()+" is not"
715  " 1-dimensional.\n"
716  "Only 1-dimensional pixels are supported by the"
717  " Healpix projection.";
719  }
720 
721  // Allocate boundaries
722  std::vector<GSkyDir> boundaries;
723 
724  // Determine pixel index and face number
725  int ix;
726  int iy;
727  int face;
728  pix2xyf(int(pixel), &ix, &iy, &face);
729 
730  // ...
731  double dc = 0.5 / m_nside;
732  double xc = (ix + 0.5) / m_nside;
733  double yc = (iy + 0.5) / m_nside;
734  double d = 1.0 / (step*m_nside);
735 
736  //
737  for (int i = 0; i < step; ++i) {
738 
739  // Declare local coordinates
740  double z;
741  double phi;
742 
743  // First coordinate
744  xyf2loc(xc+dc-i*d, yc+dc, face, &z, &phi);
745  boundaries.push_back(loc2dir(z, phi));
746 
747  // Second coordinate
748  xyf2loc(xc-dc, yc+dc-i*d, face, &z, &phi);
749  boundaries.push_back(loc2dir(z, phi));
750 
751  // Third coordinate
752  xyf2loc(xc-dc+i*d, yc-dc, face, &z, &phi);
753  boundaries.push_back(loc2dir(z, phi));
754 
755  // Forth coordinate
756  xyf2loc(xc+dc, yc-dc+i*d, face, &z, &phi);
757  boundaries.push_back(loc2dir(z, phi));
758 
759  } // endfor: looped over step size
760 
761  // Return boundaries
762  return boundaries;
763 }
764 
765 
766 /***********************************************************************//**
767  * @brief Return maximum angular distance between pixel centre and corners
768  *
769  * @return Maximum angular distance between pixel centre and corners (radians).
770  *
771  * Returns the maximum angular distance (in radians) between any pixel
772  * centre and its corners.
773  *
774  * This method has been adapted from the max_pixrad() function located in the
775  * file healpix_base.cc in Healpix version 3.20.
776  ***************************************************************************/
777 double GHealpix::max_pixrad(void) const
778 {
779  // Compute ...
780  double t1 = 1.0 - 1.0/m_nside;
781  t1 *= t1;
782 
783  // Get vectors
784  GVector va = set_z_phi(2.0/3.0, gammalib::pi/(4*m_nside));
785  GVector vb = set_z_phi(1.0 - t1/3, 0.0);
786 
787  // Get angle
788  double angle = std::atan2(norm(cross(va, vb)), va * vb);
789 
790  // Return angle
791  return angle;
792 }
793 
794 
795 /***********************************************************************//**
796  * @brief Print WCS information
797  *
798  * @param[in] chatter Chattiness (defaults to NORMAL).
799  * @return String containing WCS information.
800  ***************************************************************************/
801 std::string GHealpix::print(const GChatter& chatter) const
802 {
803  // Initialise result string
804  std::string result;
805 
806  // Continue only if chatter is not silent
807  if (chatter != SILENT) {
808 
809  // Append header
810  result.append("=== GHealpix ===");
811 
812  // Append information
813  result.append("\n"+gammalib::parformat("Coordinate system"));
814  result.append(coordsys());
815  result.append("\n"+gammalib::parformat("Nside (# of divisions)"));
816  result.append(gammalib::str(m_nside));
817  result.append("\n"+gammalib::parformat("Npface (pixels per face)"));
818  result.append(gammalib::str(m_npface));
819  result.append("\n"+gammalib::parformat("Ncap (# of cap pixels)"));
820  result.append(gammalib::str(m_ncap));
821  result.append("\n"+gammalib::parformat("Npix (# of pixels)"));
822  result.append(gammalib::str(m_num_pixels));
823  result.append("\n"+gammalib::parformat("Order"));
824  result.append(gammalib::str(m_order));
825  result.append("\n"+gammalib::parformat("Solid angle per pixel"));
826  result.append(gammalib::str(m_solidangle)+" sr");
827  result.append("\n"+gammalib::parformat("Ordering")+ordering());
828 
829  } // endif: chatter was not silent
830 
831  // Return result
832  return result;
833 }
834 
835 
836 /*==========================================================================
837  = =
838  = Private methods =
839  = =
840  ==========================================================================*/
841 
842 /***********************************************************************//**
843  * @brief Initialise class members
844  ***************************************************************************/
846 {
847  // Initialise members
848  m_nside = 0;
849  m_npface = 0;
850  m_ncap = 0;
851  m_order = 0;
852  m_ordering = 0;
853  m_num_pixels = 0;
854  m_fact1 = 0.0;
855  m_fact2 = 0.0;
856  m_solidangle = 0.0;
857 
858  // Construct conversion arrays
859  for (int m = 0; m < 0x100; ++m) {
860  ctab[m] =
861  (m&0x1 ) | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
862  | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4);
863  utab[m] =
864  (m&0x1 ) | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
865  | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7);
866  }
867 
868  // Return
869  return;
870 }
871 
872 
873 /***********************************************************************//**
874  * @brief Copy class members
875  *
876  * @param[in] proj Healpix projection.
877  ***************************************************************************/
879 {
880  // Copy attributes
881  m_nside = proj.m_nside;
882  m_npface = proj.m_npface;
883  m_ncap = proj.m_ncap;
884  m_order = proj.m_order;
885  m_ordering = proj.m_ordering;
886  m_num_pixels = proj.m_num_pixels;
887  m_fact1 = proj.m_fact1;
888  m_fact2 = proj.m_fact2;
889  m_solidangle = proj.m_solidangle;
890 
891  // Return
892  return;
893 }
894 
895 
896 /***********************************************************************//**
897  * @brief Delete class members
898  ***************************************************************************/
900 {
901  // Return
902  return;
903 }
904 
905 
906 /***********************************************************************//**
907  * @brief Returns true if argument is identical
908  *
909  * @param[in] proj Sky projection.
910  *
911  * This method is a helper for the sky projection friends.
912  ***************************************************************************/
913 bool GHealpix::compare(const GSkyProjection& proj) const
914 {
915  // Initialise result
916  bool result = false;
917 
918  // Continue only we compare to a GHealpix object
919  const GHealpix* ptr = dynamic_cast<const GHealpix*>(&proj);
920  if (ptr != NULL) {
921  result = ((m_coordsys == ptr->m_coordsys) &&
922  (m_nside == ptr->m_nside) &&
923  (m_npface == ptr->m_npface) &&
924  (m_ncap == ptr->m_ncap) &&
925  (m_order == ptr->m_order) &&
926  (m_ordering == ptr->m_ordering) &&
927  (m_num_pixels == ptr->m_num_pixels));
928  }
929 
930  // Return result
931  return result;
932 }
933 
934 
935 
936 /*==========================================================================
937  = =
938  = Low-level Healpix methods =
939  = =
940  ==========================================================================*/
941 
942 /***********************************************************************//**
943  * @brief Compress Bits
944  *
945  * @param[in] value Value.
946  * @return Compressed Bits.
947  *
948  * This method has been adapted from the compress_bits() function located in
949  * the file healpix_base.cc in Healpix version 3.20.
950  ***************************************************************************/
951 int GHealpix::compress_bits(const int& value) const
952 {
953  // Compress Bits
954  int raw = (value & 0x5555) | ((value & 0x55550000) >> 15);
955  int compressed = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
956 
957  // Return compressed value
958  return compressed;
959 }
960 
961 
962 /***********************************************************************//**
963  * @brief Spread Bits
964  *
965  * @param[in] value Compressed value.
966  * @return Spread Bits.
967  *
968  * This method has been adapted from the spread_bits() function located in
969  * the file healpix_base.cc in Healpix version 3.20.
970  ***************************************************************************/
971 int GHealpix::spread_bits(const int& value) const
972 {
973  // Spread bits
974  int spread = utab[value & 0xff] | (utab[(value >> 8) & 0xff] << 16);
975 
976  // Return spread value
977  return spread;
978 }
979 
980 
981 /***********************************************************************//**
982  * @brief Convert pixel number in to (x,y,face) tuple
983  *
984  * @param[in] pix Pixel number.
985  * @param[out] ix X index.
986  * @param[out] iy Y index.
987  * @param[out] face Face number.
988  *
989  * This method has been adapted from the pix2xyf() function located in
990  * the file healpix_base.cc in Healpix version 3.20.
991  ***************************************************************************/
992 void GHealpix::pix2xyf(const int& pix, int* ix, int* iy, int* face) const
993 {
994  // Handle ring scheme ...
995  if (m_ordering == 0) {
996  ring2xyf(pix, ix, iy, face);
997  }
998 
999  // ... or nested scheme
1000  else {
1001  nest2xyf(pix, ix, iy, face);
1002  }
1003 
1004  // Return
1005  return;
1006 }
1007 
1008 
1009 /***********************************************************************//**
1010  * @brief Convert pixel number in nested scheme to (x,y,face) tuple
1011  *
1012  * @param[in] pix Pixel number in nested scheme.
1013  * @param[out] ix X index.
1014  * @param[out] iy Y index.
1015  * @param[out] face Face number.
1016  *
1017  * This method has been adapted from the nest2xyf() function located in
1018  * the file healpix_base.cc in Healpix version 3.20.
1019  ***************************************************************************/
1020 void GHealpix::nest2xyf(const int& pix, int* ix, int* iy, int* face) const
1021 {
1022  // Compute face number
1023  *face = pix >> (2 * m_order);
1024 
1025  // Compute pixel
1026  int pixel = pix & (m_npface - 1);
1027 
1028  // Compute (x,y)
1029  *ix = compress_bits(pixel);
1030  *iy = compress_bits(pixel >> 1);
1031 
1032  // Return
1033  return;
1034 }
1035 
1036 
1037 /***********************************************************************//**
1038  * @brief Convert pixel number in ring scheme to (x,y,face) tuple
1039  *
1040  * @param[in] pix Pixel number in ring scheme
1041  * @param[out] ix X index.
1042  * @param[out] iy Y index.
1043  * @param[out] face Face number.
1044  *
1045  * This method has been adapted from the ring2xyf() function located in
1046  * the file healpix_base.cc in Healpix version 3.20.
1047  ***************************************************************************/
1048 void GHealpix::ring2xyf(const int& pix, int* ix, int* iy, int* face) const
1049 {
1050  // Declare some variables
1051  int nl2 = 2*m_nside;
1052  int iring;
1053  int iphi;
1054  int kshift;
1055  int nr;
1056 
1057  // Handle pixel in the North Polar cap
1058  if (pix < m_ncap) {
1059  iring = (1+isqrt(1+2*pix)) >> 1; // Counted from North pole
1060  iphi = (pix+1) - 2*iring*(iring-1);
1061  kshift = 0;
1062  nr = iring;
1063  *face = (iphi-1)/nr;
1064  }
1065 
1066  // Handle pixel in equatorial region
1067  else if (pix < (m_num_pixels-m_ncap)) {
1068  int ip = pix - m_ncap;
1069  int tmp = (m_order>=0) ? ip >> (m_order+2) : ip/(4*m_nside);
1070  iring = tmp + m_nside;
1071  iphi = ip - tmp * 4 * m_nside + 1;
1072  kshift = (iring + m_nside) & 1;
1073  nr = m_nside;
1074  int ire = iring - m_nside + 1;
1075  int irm = nl2 + 2 - ire;
1076  int ifm = iphi - ire/2 + m_nside - 1;
1077  int ifp = iphi - irm/2 + m_nside - 1;
1078  if (m_order >= 0) {
1079  ifm >>= m_order;
1080  ifp >>= m_order;
1081  }
1082  else {
1083  ifm /= m_nside;
1084  ifp /= m_nside;
1085  }
1086  *face = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8));
1087  }
1088 
1089  // Handle pixel in the South Polar cap
1090  else {
1091  int ip = m_num_pixels - pix;
1092  iring = (1+isqrt(2*ip-1))>>1; // Counted from South pole
1093  iphi = 4 * iring + 1 - (ip - 2*iring*(iring-1));
1094  kshift = 0;
1095  nr = iring;
1096  iring = 2 * nl2 - iring;
1097  *face = 8 + (iphi-1)/nr;
1098  }
1099 
1100  // Now compute the (ix,iy) values
1101  int irt = iring - (jrll[*face] * m_nside) + 1;
1102  int ipt = 2 * iphi - jpll[*face] * nr - kshift -1;
1103  if (ipt >= nl2) {
1104  ipt -= 8*m_nside;
1105  }
1106  *ix = (ipt-irt) >> 1;
1107  *iy = (-ipt-irt) >> 1;
1108 
1109  // Return
1110  return;
1111 }
1112 
1113 
1114 /***********************************************************************//**
1115  * @brief Convert (x,y,face) tuple to pixel number in nested scheme
1116  *
1117  * @param[in] ix X index
1118  * @param[in] iy Y index
1119  * @param[in] face Face number
1120  * @return Pixel number if nested scheme
1121  *
1122  * This method has been adapted from the xyf2nest() function located in
1123  * the file healpix_base.cc in Healpix version 3.20.
1124  ***************************************************************************/
1125 int GHealpix::xyf2nest(const int& ix, const int& iy, const int& face) const
1126 {
1127  // Computed pixel number
1128  int pix = (int(face) << (2 * m_order)) +
1129  spread_bits(ix) + (spread_bits(iy) << 1);
1130 
1131  // Return pixel number
1132  return pix;
1133 }
1134 
1135 
1136 /***********************************************************************//**
1137  * @brief Convert (x,y,face) tuple to pixel number in ring scheme
1138  *
1139  * @param[in] ix X index
1140  * @param[in] iy Y index
1141  * @param[in] face Face number
1142  * @return Pixel number if ring scheme
1143  *
1144  * This method has been adapted from the xyf2ring() function located in
1145  * the file healpix_base.cc in Healpix version 3.20.
1146  ***************************************************************************/
1147 int GHealpix::xyf2ring(const int& ix, const int& iy, const int& face) const
1148 {
1149  // Compute ring number
1150  int nl4 = 4 * m_nside;
1151  int jr = (jrll[face]*m_nside) - ix - iy - 1;
1152 
1153  // Get information about that ring
1154  int n_before;
1155  int nr;
1156  bool shifted;
1157  get_ring_info(jr, &n_before, &nr, &shifted);
1158 
1159  // Compute pixel number
1160  nr >>= 2;
1161  int kshift = 1-shifted;
1162  int jp = (jpll[face]*nr + ix - iy + 1 + kshift) / 2;
1163 
1164  //planck_assert(jp<=4*nr,"must not happen");
1165 
1166  // Assumption: if this triggers, then nl4==4*nr
1167  if (jp < 1) {
1168  jp += nl4;
1169  }
1170 
1171  // Return pixel number
1172  return (n_before + jp - 1);
1173 }
1174 
1175 
1176 /***********************************************************************//**
1177  * @brief Convert (x,y,f) tuple into local coordinates
1178  *
1179  * @param[in] x X value.
1180  * @param[in] y Y value.
1181  * @param[in] face Face number.
1182  * @param[out] z Z value.
1183  * @param[out] phi Phi value.
1184  *
1185  * This method has been adapted from the xyf2loc() function located in the
1186  * file healpix_base.cc in Healpix version 3.20.
1187  ***************************************************************************/
1188 void GHealpix::xyf2loc(const double& x, const double& y, const int& face,
1189  double* z, double* phi) const
1190 {
1191  // ...
1192  double jr = jrll[face] - x - y;
1193  double nr;
1194 
1195  // Compute z
1196  if (jr < 1) {
1197  nr = jr;
1198  double tmp = nr*nr / 3.0;
1199  *z = 1 - tmp;
1200  }
1201  else if (jr > 3) {
1202  nr = 4 - jr;
1203  double tmp = nr*nr / 3.0;
1204  *z = tmp - 1;
1205  }
1206  else {
1207  nr = 1;
1208  *z = (2.0-jr) * 2.0/3.0;
1209  }
1210 
1211  // Compute Phi
1212  double tmp = jpll[face] * nr + x - y;
1213  if (tmp < 0) {
1214  tmp += 8;
1215  }
1216  if (tmp >= 8) {
1217  tmp -= 8;
1218  }
1219  *phi = (nr < 1.0e-15) ? 0.0 : (0.5 * gammalib::pihalf * tmp) / nr;
1220 
1221  // Return
1222  return;
1223 }
1224 
1225 
1226 /***********************************************************************//**
1227  * @brief Convert local coordinate into sky direction
1228  *
1229  * @param[in] z Z value.
1230  * @param[in] phi Phi value.
1231  * @return Sky direction.
1232  *
1233  * This method has been adapted from the locToVec3() function located in
1234  * the file healpix_base.cc in Healpix version 3.20.
1235  ***************************************************************************/
1236 GSkyDir GHealpix::loc2dir(const double& z, const double& phi) const
1237 {
1238  // Compute longitude and latitude
1239  double sintheta = std::sqrt((1.0 - z) * (1.0 + z));
1240  double longitude = std::atan2(sintheta * std::sin(phi), sintheta * std::cos(phi));
1241  double latitude = std::asin(z);
1242 
1243  // Set sky direction
1244  GSkyDir dir;
1245  switch (m_coordsys) {
1246  case 0:
1247  dir.radec(longitude, latitude);
1248  break;
1249  case 1:
1250  dir.lb(longitude, latitude);
1251  break;
1252  default:
1253  break;
1254  }
1255 
1256  // Return sky direction
1257  return dir;
1258 }
1259 
1260 
1261 /***********************************************************************//**
1262  * @brief Convert nside to order
1263  *
1264  * @param[in] nside Number of sides.
1265  * @return Order of HealPix projection.
1266  ***************************************************************************/
1267 int GHealpix::nside2order(const int& nside) const
1268 {
1269  // Initialise order
1270  int order = -1;
1271 
1272  // Determine order
1273  for (int m = 0; m <= order_max; ++m) {
1274  int nstest = 1 << m;
1275  if (nside == nstest) {
1276  order = m;
1277  break;
1278  }
1279  if (nside < nstest)
1280  break;
1281  }
1282 
1283  // Return order
1284  return order;
1285 }
1286 
1287 
1288 /***********************************************************************//**
1289  * @brief Converts pixel number in nested indexing scheme to ring scheme
1290  *
1291  * @param[in] pix Pixel number in nested indexing scheme.
1292  * @return Pixel number in ring indexing scheme.
1293  *
1294  * @exception GException::invalid_value
1295  * Healpix projection does not represent a hiearachical map.
1296  *
1297  * This method has been adapted from the nest2ring() function located in
1298  * the file healpix_base.cc in Healpix version 3.20.
1299  ***************************************************************************/
1300 int GHealpix::nest2ring(const int& pix) const
1301 {
1302  // Throw an exception if map is not a hierachical map
1303  if (m_order < 0) {
1304  std::string msg = "A hierarchical map projection is required.";
1306  }
1307 
1308  // Convert nested index to (x,y,face) tuple
1309  int ix;
1310  int iy;
1311  int face;
1312  nest2xyf(pix, &ix, &iy, &face);
1313 
1314  // Convert (x,y,face) tuple to ring index
1315  int iring = xyf2ring(ix, iy, face);
1316 
1317  // Return ring index
1318  return iring;
1319 }
1320 
1321 
1322 /***********************************************************************//**
1323  * @brief Converts pixel number in ring indexing scheme to nested scheme
1324  *
1325  * @param[in] pix Pixel number in ring indexing scheme.
1326  * @return Pixel number in nested indexing scheme.
1327  *
1328  * @exception GException::invalid_value
1329  * Healpix projection does not represent a hiearachical map.
1330  *
1331  * This method has been adapted from the ring2nest() function located in
1332  * the file healpix_base.cc in Healpix version 3.20.
1333  ***************************************************************************/
1334 int GHealpix::ring2nest(const int& pix) const
1335 {
1336  // Throw an exception if map is not a hierachical map
1337  if (m_order < 0) {
1338  std::string msg = "A hierarchical map projection is required.";
1340  }
1341 
1342  // Convert nested index to (x,y,face) tuple
1343  int ix;
1344  int iy;
1345  int face;
1346  ring2xyf(pix, &ix, &iy, &face);
1347 
1348  // Convert (x,y,face) tuple to ring index
1349  int iring = xyf2nest(ix, iy, face);
1350 
1351  // Return ring index
1352  return iring;
1353 }
1354 
1355 
1356 /***********************************************************************//**
1357  * @brief Convert pixel index to (x,y) coordinate
1358  *
1359  * @param[in] ipix Pixel index for which (x,y) are to be computed.
1360  * @param[out] x Pointer to x coordinate.
1361  * @param[out] y Pointer to y coordinate.
1362  ***************************************************************************/
1363 void GHealpix::pix2xy(const int& ipix, int* x, int* y) const
1364 {
1365  // Set x coordinate
1366  int raw = (ipix & 0x5555) | ((ipix & 0x55550000) >> 15);
1367  *x = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
1368 
1369  // Set y coordinate
1370  raw = ((ipix & 0xaaaa) >> 1) | ((ipix & 0xaaaa0000) >> 16);
1371  *y = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
1372 
1373  // Return
1374  return;
1375 }
1376 
1377 
1378 /***********************************************************************//**
1379  * @brief Convert (x,y) coordinate to pixel
1380  *
1381  * @param[in] x x coordinate.
1382  * @param[in] y y coordinate.
1383  ***************************************************************************/
1384 int GHealpix::xy2pix(int x, int y) const
1385 {
1386  // Return pixel
1387  return utab[x&0xff] | (utab[x>>8]<<16) | (utab[y&0xff]<<1) | (utab[y>>8]<<17);
1388 }
1389 
1390 
1391 /***********************************************************************//**
1392  * @brief Convert pixel index to (theta,phi) angles for ring ordering
1393  *
1394  * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
1395  * @param[out] theta Pointer to result zenith angle in radians.
1396  * @param[out] phi Pointer to result azimuth angle in radians.
1397  *
1398  * @exception GException::out_of_range
1399  * Pixel index is out of range.
1400  ***************************************************************************/
1401 void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi) const
1402 {
1403  // Check if ipix is in range
1404  if (ipix < 0 || ipix >= m_num_pixels) {
1406  }
1407 
1408  // Handle North Polar cap
1409  if (ipix < m_ncap) {
1410  int iring = int(0.5*(1+isqrt(1+2*ipix))); // counted from North pole
1411  int iphi = (ipix+1) - 2*iring*(iring-1);
1412  *theta = std::acos(1.0 - (iring*iring) * m_fact2);
1413  *phi = (iphi - 0.5) * gammalib::pi/(2.0*iring);
1414  }
1415 
1416  // Handle Equatorial region
1417  else if (ipix < (m_num_pixels - m_ncap)) {
1418  int ip = ipix - m_ncap;
1419  int iring = ip/(4*m_nside) + m_nside; // counted from North pole
1420  int iphi = ip%(4*m_nside) + 1;
1421  double fodd = ((iring+m_nside)&1) ? 1 : 0.5;
1422  int nl2 = 2*m_nside;
1423  *theta = std::acos((nl2 - iring) * m_fact1);
1424  *phi = (iphi - fodd) * gammalib::pi/nl2;
1425  }
1426 
1427  // Handle South Polar cap
1428  else {
1429  int ip = m_num_pixels - ipix;
1430  int iring = int(0.5*(1+isqrt(2*ip-1))); // Counted from South pole
1431  int iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
1432  *theta = std::acos(-1.0 + (iring*iring) * m_fact2);
1433  *phi = (iphi - 0.5) * gammalib::pi/(2.*iring);
1434  }
1435 
1436  // Return
1437  return;
1438 }
1439 
1440 
1441 /***********************************************************************//**
1442  * @brief Convert pixel index to (theta,phi) angles for nested ordering
1443  *
1444  * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
1445  * @param[out] theta Pointer to result zenith angle in radians.
1446  * @param[out] phi Pointer to result azimuth angle in radians.
1447  *
1448  * @exception GException::out_of_range
1449  * Pixel index is out of range.
1450  ***************************************************************************/
1451 void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) const
1452 {
1453  // Check if ipix is in range
1454  if (ipix < 0 || ipix >= m_num_pixels) {
1456  }
1457 
1458  // Get face number and index in face
1459  int nl4 = 4 * m_nside;
1460  int face_num = ipix >> (2*m_order); // Face number in {0,11}
1461  int ipf = ipix & (m_npface - 1);
1462 
1463  // Get pixel coordinates
1464  int ix;
1465  int iy;
1466  pix2xy(ipf, &ix, &iy);
1467 
1468  // Computes the z coordinate on the sphere
1469  int jr = (jrll[face_num] << m_order) - ix - iy - 1;
1470 
1471  // Declare result variables
1472  int nr;
1473  double z;
1474  int kshift;
1475 
1476  // North pole region
1477  if (jr < m_nside) {
1478  nr = jr;
1479  z = 1. - nr*nr*m_fact2;
1480  kshift = 0;
1481  }
1482 
1483  // South pole region
1484  else if (jr > 3*m_nside) {
1485  nr = nl4 - jr;
1486  z = nr*nr*m_fact2 - 1;
1487  kshift = 0;
1488  }
1489 
1490  // Equatorial region
1491  else {
1492  nr = m_nside;
1493  z = (2*m_nside-jr) * m_fact1;
1494  kshift = (jr-m_nside) & 1;
1495  }
1496 
1497  // Computes the phi coordinate on the sphere, in [0,2Pi]
1498  int jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
1499  if (jp > nl4) jp -= nl4;
1500  if (jp < 1) jp += nl4;
1501 
1502  // Computes Theta and Phi
1503  *theta = std::acos(z);
1504  *phi = (jp - (kshift+1)*0.5) * (gammalib::pihalf / nr);
1505 
1506  // Return
1507  return;
1508 }
1509 
1510 
1511 /***********************************************************************//**
1512  * @brief Returns pixel which contains angular coordinates (z,phi)
1513  *
1514  * @param[in] z Cosine of zenith angle - cos(theta).
1515  * @param[in] phi Azimuth angle in radians.
1516  ***************************************************************************/
1517 int GHealpix::ang2pix_z_phi_ring(double z, double phi) const
1518 {
1519  // Initialise pixel
1520  int ipix = 0;
1521 
1522  // Setup
1523  double za = fabs(z);
1524  double tt = gammalib::modulo(phi, gammalib::twopi) * gammalib::inv_pihalf; // in [0,4)
1525 
1526  // Equatorial region
1527  if (za <= gammalib::twothird) {
1528  double temp1 = m_nside*(0.5+tt);
1529  double temp2 = m_nside*z*0.75;
1530  int jp = int(temp1-temp2); // index of ascending edge line
1531  int jm = int(temp1+temp2); // index of descending edge line
1532  int ir = m_nside + 1 + jp - jm; // in {1,2n+1}
1533  int kshift = 1 - (ir & 1); // kshift=1 if ir even, 0 otherwise
1534  int ip = (jp+jm-m_nside+kshift+1)/2; // in {0,4n-1}
1535  ip = int(gammalib::modulo(ip,4*m_nside));
1536  ipix = m_ncap + (ir-1)*4*m_nside + ip;
1537  }
1538 
1539  // North & South polar caps
1540  else {
1541  double tp = tt - int(tt);
1542  double tmp = m_nside * std::sqrt(3*(1-za));
1543  int jp = int(tp*tmp); // increasing edge line index
1544  int jm = int((1.0-tp)*tmp); // decreasing edge line index
1545  int ir = jp + jm + 1; // ring number counted from the closest pole
1546  int ip = int(tt*ir); // in {0,4*ir-1}
1547  ip = int(gammalib::modulo(ip,4*ir));
1548  if (z>0)
1549  ipix = 2*ir*(ir-1) + ip;
1550  else
1551  ipix = m_num_pixels - 2*ir*(ir+1) + ip;
1552  }
1553 
1554  // Return pixel
1555  return ipix;
1556 }
1557 
1558 
1559 /***********************************************************************//**
1560  * @brief Returns pixel which contains angular coordinates (z,phi)
1561  *
1562  * @param[in] z Cosine of zenith angle - cos(theta).
1563  * @param[in] phi Azimuth angle in radians.
1564  ***************************************************************************/
1565 int GHealpix::ang2pix_z_phi_nest(double z, double phi) const
1566 {
1567  // Initialise face and pixel numbers
1568  int face_num;
1569  int ix;
1570  int iy;
1571 
1572  // Setup
1573  double za = fabs(z);
1574  double tt = gammalib::modulo(phi, gammalib::twopi) * gammalib::inv_pihalf; // in [0,4)
1575 
1576  // Equatorial region
1577  if (za <= gammalib::twothird) {
1578  double temp1 = ns_max*(0.5+tt);
1579  double temp2 = ns_max*z*0.75;
1580  int jp = int(temp1-temp2); // index of ascending edge line
1581  int jm = int(temp1+temp2); // index of descending edge line
1582  int ifp = jp >> order_max; // in {0,4}
1583  int ifm = jm >> order_max;
1584  if (ifp == ifm) // faces 4 to 7
1585  face_num = (ifp==4) ? 4: ifp+4;
1586  else if (ifp < ifm) // (half-)faces 0 to 3
1587  face_num = ifp;
1588  else // (half-)faces 8 to 11
1589  face_num = ifm + 8;
1590  ix = jm & (ns_max-1);
1591  iy = ns_max - (jp & (ns_max-1)) - 1;
1592  }
1593 
1594  // Polar region, za > 2/3
1595  else {
1596  int ntt = int(tt);
1597  double tp = tt-ntt;
1598  double tmp = ns_max * std::sqrt(3*(1-za));
1599  int jp = int(tp*tmp); // increasing edge line index
1600  int jm = int((1.0-tp)*tmp); // decreasing edge line index
1601  if (jp >= ns_max) jp = ns_max-1; // for points too close to the boundary
1602  if (jm >= ns_max) jm = ns_max-1;
1603  if (z >= 0) {
1604  face_num = ntt; // in {0,3}
1605  ix = ns_max - jm - 1;
1606  iy = ns_max - jp - 1;
1607  }
1608  else {
1609  face_num = ntt + 8; // in {8,11}
1610  ix = jp;
1611  iy = jm;
1612  }
1613  }
1614 
1615  // Get pixel
1616  int ipf = xy2pix(ix, iy);
1617  ipf >>= (2*(order_max - m_order)); // in {0, nside**2 - 1}
1618  return ipf + (face_num<<(2*m_order)); // in {0, 12*nside**2 - 1}
1619 }
1620 
1621 
1622 /***********************************************************************//**
1623  * @brief Get ring index north of cos(theta)
1624  *
1625  * @param[in] costheta Cosine of colatitude theta
1626  *
1627  * Returns the number of the next ring to the north of @a cos(theta). It may
1628  * return 0; in this case @a costheta lies north of all rings.
1629  *
1630  * This method has been adapted from the ring_above() function located in
1631  * the file healpix_base.cc in Healpix version 3.20.
1632  ***************************************************************************/
1633 int GHealpix::ring_above(const double& costheta) const
1634 {
1635  // Initialise ring index
1636  int iring;
1637 
1638  // Get absolute cosine
1639  double acostheta = std::abs(costheta);
1640 
1641  // Are we in the equatorial region
1642  if (acostheta <= gammalib::twothird) {
1643  iring = int(m_nside * (2.0 - 1.5 * costheta));
1644  }
1645 
1646  // ... otherwise we're in the pole region
1647  else {
1648  iring = int(m_nside * std::sqrt(3.0 * (1.0 - acostheta)));
1649  if (costheta <= 0.0) {
1650  iring = 4 * m_nside - iring - 1;
1651  }
1652  }
1653 
1654  // Return
1655  return iring;
1656 }
1657 
1658 
1659 /***********************************************************************//**
1660  * @brief Returns useful information about a given ring of the projection
1661  *
1662  * @param[in] ring The ring number (the number of the first ring is 1)
1663  * @param[out] startpix The number of the first pixel in the ring
1664  * @param[out] ringpix The number of pixels in the ring
1665  * @param[out] shifted If @a true, the center of the first pixel is not at
1666  * @a phi=0
1667  *
1668  * This method has been adapted from the get_ring_info_small() function
1669  * located in the file healpix_base.cc in Healpix version 3.20.
1670  ***************************************************************************/
1671 void GHealpix::get_ring_info(const int& ring,
1672  int* startpix,
1673  int* ringpix,
1674  bool* shifted) const
1675 {
1676  // Handle ring in North polar cap
1677  if (ring < m_nside) {
1678  *shifted = true;
1679  *ringpix = 4 * ring;
1680  *startpix = 2 * ring * (ring-1);
1681  }
1682 
1683  // Handle ring in equatorial region
1684  else if (ring < 3*m_nside) {
1685  *shifted = ((ring-m_nside) & 1) == 0;
1686  *ringpix = 4 * m_nside;
1687  *startpix = m_ncap + (ring-m_nside) * *ringpix;
1688  }
1689 
1690  // Handle ring in South polar cap
1691  else {
1692  int nr = 4 * m_nside - ring;
1693  *shifted = true;
1694  *ringpix = 4 * nr;
1695  *startpix = m_num_pixels - 2 * nr * (nr+1);
1696  }
1697 
1698  // Return
1699  return;
1700 }
1701 
1702 
1703 /***********************************************************************//**
1704  * @brief Returns useful information about a given ring of the projection
1705  *
1706  * @param[in] ring The ring number (the number of the first ring is 1)
1707  * @param[out] startpix The number of the first pixel in the ring
1708  * @param[out] ringpix The number of pixels in the ring
1709  * @param[out] theta (radians)
1710  * @param[out] shifted If @a true, the center of the first pixel is not at
1711  * @a phi=0
1712  *
1713  * This method has been adapted from the get_ring_info2() function
1714  * located in the file healpix_base.cc in Healpix version 3.20.
1715  ***************************************************************************/
1716 void GHealpix::get_ring_info(const int& ring,
1717  int* startpix,
1718  int* ringpix,
1719  double* theta,
1720  bool* shifted) const
1721 {
1722  //
1723  int northring = (ring > 2 * m_nside) ? 4 * m_nside - ring : ring;
1724 
1725  // Are we in the North?
1726  if (northring < m_nside) {
1727  double tmp = northring * northring * m_fact2;
1728  double costheta = 1.0 - tmp;
1729  double sintheta = std::sqrt(tmp * (2.0-tmp));
1730  *startpix = 2 * northring * (northring - 1);
1731  *ringpix = 4 * northring;
1732  *theta = std::atan2(sintheta, costheta);
1733  *shifted = true;
1734  }
1735  else {
1736  *theta = std::acos((2.0 * m_nside-northring) * m_fact1);
1737  *ringpix = 4 * m_nside;
1738  *shifted = ((northring - m_nside) & 1) == 0;
1739  *startpix = m_ncap + (northring - m_nside) * *ringpix;
1740  }
1741 
1742  // Are we in the southern hemisphere?
1743  if (northring != ring) {
1744  *theta = gammalib::pi - *theta;
1745  *startpix = m_num_pixels - *startpix - *ringpix;
1746  }
1747 
1748  // Return
1749  return;
1750 }
1751 
1752 
1753 /***********************************************************************//**
1754  * @brief Return interpolator
1755  *
1756  * @param[in] theta Colatitude of direction (radian, the North pole is at theta=0)
1757  * @param[in] phi Longitude of direction (radians)
1758  *
1759  * @exception GException::invalid_argument
1760  * Invalid @p theta argument
1761  *
1762  * Returns a bilinear pixel interpolator for a given sky direction.
1763  *
1764  * This method has been adapted from the get_interpol() function
1765  * located in the file healpix_base.cc in Healpix version 3.20.
1766  ***************************************************************************/
1767 GBilinear GHealpix::interpolator(const double& theta, const double& phi) const
1768 {
1769  // Allocate interpolator
1771 
1772  // Check that theta is valid
1773  if (theta < 0.0 || theta > gammalib::pi) {
1774  std::string msg = "Colatitude "+gammalib::str(theta)+" is outside "
1775  "valid range [0,pi].";
1777  }
1778 
1779  // Prepare computation
1780  double costheta = std::cos(theta);
1781  int ir1 = ring_above(costheta); // Ring above actual colatitude
1782  int ir2 = ir1 + 1; // Ring below actual colatitude
1783  int sp; // Start pixel in ring
1784  int nr; // Number of pixels in ring
1785  double theta1; // Colatitude of ring above
1786  double theta2; // Colatitude of ring below
1787  bool shift;
1788 
1789  // Compute interpolating pixels and phi weights if the colatitude is not
1790  // North of all rings (if we're North of all rings, ir1=0)
1791  if (ir1 > 0) {
1792  get_ring_info(ir1, &sp, &nr, &theta1, &shift);
1793  double dphi = gammalib::twopi / nr; // Phi spacing of pixels
1794  double tmp = (phi/dphi - 0.5*shift);
1795  int i1 = (tmp < 0.0) ? int(tmp)-1 : int(tmp);
1796  double w1 = (phi - (i1+0.5*shift)*dphi) / dphi;
1797  int i2 = i1 + 1;
1798  if (i1 < 0) {
1799  i1 += nr;
1800  }
1801  if (i2 >= nr) {
1802  i2 -= nr;
1803  }
1804  interpolator.index1() = sp + i1;
1805  interpolator.index2() = sp + i2;
1806  interpolator.weight1() = 1.0 - w1;
1807  interpolator.weight2() = w1;
1808  }
1809 
1810  // Compute interpolating pixels and phi weights is the colatitude is not
1811  // South of all rings (if we're South of all rings, ir2=4*m_nside)
1812  if (ir2 < (4*m_nside)) {
1813  get_ring_info(ir2, &sp, &nr, &theta2, &shift);
1814  double dphi = gammalib::twopi / nr; // Phi spacing of pixels
1815  double tmp = (phi/dphi - 0.5*shift);
1816  int i1 = (tmp < 0.0) ? int(tmp)-1 : int(tmp);
1817  double w1 = (phi - (i1+0.5*shift)*dphi) / dphi;
1818  int i2 = i1 + 1;
1819  if (i1 < 0) {
1820  i1 += nr;
1821  }
1822  if (i2 >= nr) {
1823  i2 -= nr;
1824  }
1825  interpolator.index3() = sp + i1;
1826  interpolator.index4() = sp + i2;
1827  interpolator.weight3() = 1.0 - w1;
1828  interpolator.weight4() = w1;
1829  }
1830 
1831  // Now handle the special case that the colatitude is North of all
1832  // rings
1833  if (ir1 == 0) {
1834  double wtheta = theta/theta2;
1835  interpolator.weight3() *= wtheta;
1836  interpolator.weight4() *= wtheta;
1837  double fac = (1.0-wtheta)*0.25;
1838  interpolator.weight1() = fac;
1839  interpolator.weight2() = fac;
1840  interpolator.weight3() += fac;
1841  interpolator.weight4() += fac;
1842  interpolator.index1() = (interpolator.index3() + 2) & 3;
1843  interpolator.index2() = (interpolator.index4() + 2) & 3;
1844  }
1845 
1846  // ... and now the case that the colatitude is South of all rings
1847  else if (ir2 == 4*m_nside) {
1848  double wtheta = (theta-theta1) / (gammalib::pi-theta1);
1849  interpolator.weight1() *= (1.0 - wtheta);
1850  interpolator.weight2() *= (1.0 - wtheta);
1851  double fac = wtheta*0.25;
1852  interpolator.weight1() += fac;
1853  interpolator.weight2() += fac;
1854  interpolator.weight3() = fac;
1855  interpolator.weight4() = fac;
1856  interpolator.index3() = ((interpolator.index1() + 2) & 3) + m_num_pixels - 4;
1857  interpolator.index4() = ((interpolator.index2() + 2) & 3) + m_num_pixels - 4;
1858  }
1859 
1860  // ... and now multiply-in the theta weights for the general case
1861  else {
1862  double wtheta = (theta-theta1) / (theta2-theta1);
1863  interpolator.weight1() *= (1.0 - wtheta);
1864  interpolator.weight2() *= (1.0 - wtheta);
1865  interpolator.weight3() *= wtheta;
1866  interpolator.weight4() *= wtheta;
1867  }
1868 
1869  // If we have a nested pixel scheme then convert now the ring
1870  // indices into nested indices
1871  if (m_ordering == 1) {
1872  interpolator.index1() = ring2nest(interpolator.index1());
1873  interpolator.index2() = ring2nest(interpolator.index2());
1874  interpolator.index3() = ring2nest(interpolator.index3());
1875  interpolator.index4() = ring2nest(interpolator.index4());
1876  }
1877 
1878  // Return interpolator
1879  return interpolator;
1880 }
1881 
1882 
1883 /***********************************************************************//**
1884  * @brief Integer n that fulfills n*n <= arg < (n+1)*(n+1)
1885  *
1886  * @param[in] arg Argument.
1887  *
1888  * Returns the integer @a n, which fulfills @a n*n <= arg < (n+1)*(n+1).
1889  ***************************************************************************/
1890 unsigned int GHealpix::isqrt(unsigned int arg) const
1891 {
1892  // Return
1893  return unsigned(std::sqrt(arg+0.5));
1894 }
1895 
1896 
1897 /***********************************************************************//**
1898  * @brief Return 3D vector
1899  *
1900  * @param[in] z Z value.
1901  * @param[in] phi Phi value.
1902  * @return 3D vector
1903  *
1904  * This method has been adapted from the set_z_phi() function located in the
1905  * file vec3.h in Healpix version 3.20.
1906  ***************************************************************************/
1907 GVector GHealpix::set_z_phi(const double& z, const double& phi) const
1908 {
1909  // Initialise 3D vector
1910  GVector vector(3);
1911 
1912  // Assign elements
1913  double sintheta = std::sqrt((1.0 - z) * (1.0 + z));
1914  vector[0] = sintheta * std::cos(phi);
1915  vector[1] = sintheta * std::sin(phi);
1916  vector[2] = z;
1917 
1918  // Return vector
1919  return vector;
1920 }
#define G_NEST2RING
Definition: GHealpix.cpp:43
static short ctab[0x100]
Definition: GHealpix.cpp:87
double m_solidangle
Solid angle of pixel.
Definition: GHealpix.hpp:131
void nest2xyf(const int &pix, int *ix, int *iy, int *face) const
Convert pixel number in nested scheme to (x,y,face) tuple.
Definition: GHealpix.cpp:1020
const int nb_yoffset[]
Definition: GHealpix.cpp:64
std::string print(const GChatter &chatter=NORMAL) const
Print pixel.
Definition: GSkyPixel.cpp:292
void init_members(void)
Initialise class members.
Definition: GHealpix.cpp:845
virtual void clear(void)
Clear object.
Definition: GHealpix.cpp:255
GHealpix(void)
Void constructor.
Definition: GHealpix.cpp:102
double m_fact1
Definition: GHealpix.hpp:129
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
bool is_1D(void) const
Check if pixel is 1D.
Definition: GSkyPixel.hpp:157
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const int & nside(void) const
Returns number of divisions of the side of each base pixel.
Definition: GHealpix.hpp:204
std::vector< int > neighbours(const GSkyPixel &pixel) const
Return neighbouring pixels of a pixel.
Definition: GHealpix.cpp:577
const double pi
Definition: GMath.hpp:35
void copy_members(const GHealpix &wcs)
Copy class members.
Definition: GHealpix.cpp:878
GHealpix & operator=(const GHealpix &wcs)
Assignment operator.
Definition: GHealpix.cpp:220
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
void lb(const double &l, const double &b)
Set galactic sky direction (radians)
Definition: GSkyDir.cpp:229
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
virtual GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel of sky coordinate.
Definition: GHealpix.cpp:434
int & index2(void)
Access index 2.
Definition: GBilinear.hpp:114
const double & b(void) const
Return galactic latitude in radians.
Definition: GSkyDir.hpp:181
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
int m_order
Order.
Definition: GHealpix.hpp:126
virtual void read(const GFitsHDU &hdu)
Read Healpix definition from FITS header.
Definition: GHealpix.cpp:291
int m_nside
Number of divisions of each base pixel (1-8192)
Definition: GHealpix.hpp:123
unsigned int isqrt(unsigned int arg) const
Integer n that fulfills n*n &lt;= arg &lt; (n+1)*(n+1)
Definition: GHealpix.cpp:1890
double modulo(const double &v1, const double &v2)
Returns the remainder of the division.
Definition: GMath.cpp:526
double max_pixrad(void) const
Return maximum angular distance between pixel centre and corners.
Definition: GHealpix.cpp:777
virtual GBilinear interpolator(const GSkyDir &dir) const
Return interpolator for given sky direction.
Definition: GHealpix.cpp:475
virtual std::string print(const GChatter &chatter=NORMAL) const
Print WCS information.
Definition: GHealpix.cpp:801
void pix2ang_ring(int ipix, double *theta, double *phi) const
Convert pixel index to (theta,phi) angles for ring ordering.
Definition: GHealpix.cpp:1401
double m_fact2
Definition: GHealpix.hpp:130
#define G_INTERPOLATOR
Definition: GHealpix.cpp:48
static short utab[0x100]
Definition: GHealpix.cpp:88
virtual GSkyProjection & operator=(const GSkyProjection &proj)
Assignment operator.
void ring2xyf(const int &pix, int *ix, int *iy, int *face) const
Convert pixel number in ring scheme to (x,y,face) tuple.
Definition: GHealpix.cpp:1048
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
void get_ring_info(const int &ring, int *startpix, int *ringpix, bool *shifted) const
Returns useful information about a given ring of the projection.
Definition: GHealpix.cpp:1671
double & weight2(void)
Access weight 2.
Definition: GBilinear.hpp:162
const int nb_swaparray[][3]
Definition: GHealpix.cpp:74
Gammalib tools definition.
int nside2order(const int &nside) const
Convert nside to order.
Definition: GHealpix.cpp:1267
#define G_CONSTRUCT
Definition: GHealpix.cpp:39
int m_ordering
Pixel ordering (0=ring, 1=nested, -1=?)
Definition: GHealpix.hpp:127
#define G_BOUNDARIES
Definition: GHealpix.cpp:50
int xyf2ring(const int &ix, const int &iy, const int &face) const
Convert (x,y,face) tuple to pixel number in ring scheme.
Definition: GHealpix.cpp:1147
#define G_ORDERING_SET
Definition: GHealpix.cpp:47
void pix2ang_nest(int ipix, double *theta, double *phi) const
Convert pixel index to (theta,phi) angles for nested ordering.
Definition: GHealpix.cpp:1451
HealPix projection class definition.
const double & ra(void) const
Return Right Ascension in radians.
Definition: GSkyDir.hpp:208
int m_ncap
Number of pixels in polar cap.
Definition: GHealpix.hpp:125
const int jpll[12]
Definition: GHealpix.cpp:62
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
Definition: GSkyDir.cpp:175
Sky map pixel class.
Definition: GSkyPixel.hpp:74
const double pihalf
Definition: GMath.hpp:38
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_PIX2ANG_NEST
Definition: GHealpix.cpp:46
#define G_PIX2ANG_RING
Definition: GHealpix.cpp:45
#define G_NEIGHBOURS
Definition: GHealpix.cpp:49
virtual GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GHealpix.cpp:385
HealPix projection class interface definition.
Definition: GHealpix.hpp:47
const int order_max
Definition: GHealpix.cpp:83
int spread_bits(const int &value) const
Spread Bits.
Definition: GHealpix.cpp:971
std::vector< GSkyDir > boundaries(const GSkyPixel &pixel, const int &step=1) const
Return pixel boundaries.
Definition: GHealpix.cpp:709
int ang2pix_z_phi_nest(double z, double phi) const
Returns pixel which contains angular coordinates (z,phi)
Definition: GHealpix.cpp:1565
GVector set_z_phi(const double &z, const double &phi) const
Return 3D vector.
Definition: GHealpix.cpp:1907
GVector cross(const GVector &a, const GVector &b)
Vector cross product.
Definition: GVector.cpp:762
Bilinear interpolator class.
Definition: GBilinear.hpp:40
const double inv_pihalf
Definition: GMath.hpp:39
const double twothird
Definition: GMath.hpp:49
int & index1(void)
Access index 1.
Definition: GBilinear.hpp:102
const int nb_facearray[][12]
Definition: GHealpix.cpp:65
#define G_RING2NEST
Definition: GHealpix.cpp:44
int ring2nest(const int &pix) const
Converts pixel number in ring indexing scheme to nested scheme.
Definition: GHealpix.cpp:1334
int & index4(void)
Access index 4.
Definition: GBilinear.hpp:138
int m_coordsys
0=EQU, 1=GAL
void pix2xy(const int &ipix, int *x, int *y) const
Convert pixel index to (x,y) coordinate.
Definition: GHealpix.cpp:1363
GChatter
Definition: GTypemaps.hpp:33
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
const int ns_max
Definition: GHealpix.cpp:84
const double & l(void) const
Return galactic longitude in radians.
Definition: GSkyDir.hpp:154
const double & dec(void) const
Return Declination in radians.
Definition: GSkyDir.hpp:235
#define G_READ
Definition: GHealpix.cpp:40
int compress_bits(const int &value) const
Compress Bits.
Definition: GHealpix.cpp:951
const int & npix(void) const
Returns number of pixels.
Definition: GHealpix.hpp:193
int nest2ring(const int &pix) const
Converts pixel number in nested indexing scheme to ring scheme.
Definition: GHealpix.cpp:1300
Vector class interface definition.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
int xyf2nest(const int &ix, const int &iy, const int &face) const
Convert (x,y,face) tuple to pixel number in nested scheme.
Definition: GHealpix.cpp:1125
GVector asin(const GVector &vector)
Computes arcsin of vector elements.
Definition: GVector.cpp:1016
double & weight3(void)
Access weight 3.
Definition: GBilinear.hpp:174
virtual std::string coordsys(void) const
Returns coordinate system.
void free_members(void)
Delete class members.
const int nb_xoffset[]
Definition: GHealpix.cpp:63
int xy2pix(int x, int y) const
Convert (x,y) coordinate to pixel.
Definition: GHealpix.cpp:1384
int & index3(void)
Access index 3.
Definition: GBilinear.hpp:126
virtual bool compare(const GSkyProjection &proj) const
Returns true if argument is identical.
Definition: GHealpix.cpp:913
GSkyDir loc2dir(const double &z, const double &phi) const
Convert local coordinate into sky direction.
Definition: GHealpix.cpp:1236
double & weight4(void)
Access weight 4.
Definition: GBilinear.hpp:186
int ang2pix_z_phi_ring(double z, double phi) const
Returns pixel which contains angular coordinates (z,phi)
Definition: GHealpix.cpp:1517
void xyf2loc(const double &x, const double &y, const int &face, double *z, double *phi) const
Convert (x,y,f) tuple into local coordinates.
Definition: GHealpix.cpp:1188
double & weight1(void)
Access weight 1.
Definition: GBilinear.hpp:150
int m_npface
Definition: GHealpix.hpp:124
std::string ordering(void) const
Returns ordering parameter.
Definition: GHealpix.cpp:512
virtual ~GHealpix(void)
Destructor.
Definition: GHealpix.cpp:198
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
void pix2xyf(const int &pix, int *ix, int *iy, int *face) const
Convert pixel number in to (x,y,face) tuple.
Definition: GHealpix.cpp:992
void init_members(void)
Initialise class members.
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:820
double atan2(const double &y, const double &x)
Compute arc tangens in radians.
Definition: GMath.cpp:102
const double fourpi
Definition: GMath.hpp:37
const double twopi
Definition: GMath.hpp:36
Vector class.
Definition: GVector.hpp:46
Abstract sky projection base class.
int m_num_pixels
Number of pixels in projection.
Definition: GHealpix.hpp:128
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
virtual GHealpix * clone(void) const
Clone instance.
Definition: GHealpix.cpp:273
int ring_above(const double &z) const
Get ring index north of cos(theta)
Definition: GHealpix.cpp:1633
virtual void write(GFitsHDU &hdu) const
Write Healpix definition into FITS HDU.
Definition: GHealpix.cpp:357
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
const int jrll[12]
Definition: GHealpix.cpp:61
#define G_XY2DIR
Definition: GHealpix.cpp:41
Mathematical function definitions.
void free_members(void)
Delete class members.
Definition: GHealpix.cpp:899
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413