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