Bug Summary

File:src/sky/GHealpix.cpp
Location:line 421, column 17
Description:Function call argument is an uninitialized value

Annotated Source Code

1/***************************************************************************
2 * GHealpix.cpp - Healpix projection class *
3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2010-2013 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_H1
29#include <config.h>
30#endif
31#include <cmath>
32#include "GException.hpp"
33#include "GTools.hpp"
34#include "GMath.hpp"
35#include "GHealpix.hpp"
36
37/* __ Method name definitions ____________________________________________ */
38#define G_CONSTRUCT"GHealpix::GHealpix(int& ,std::string& ,std::string&)" "GHealpix::GHealpix(int& ,std::string& ,std::string&)"
39#define G_READ"GHealpix::read(GFitsHDU&)" "GHealpix::read(GFitsHDU&)"
40#define G_XY2DIR"GHealpix::xy2dir(GSkyPixel&)" "GHealpix::xy2dir(GSkyPixel&)"
41#define G_DIR2XY2"GHealpix::dir2xy(GSkyDir&)" "GHealpix::dir2xy(GSkyDir&)"
42#define G_PIX2ANG_RING"GHealpix::pix2ang_ring(int, double*, double*)" "GHealpix::pix2ang_ring(int, double*, double*)"
43#define G_PIX2ANG_NEST"GHealpix::pix2ang_nest(int, double*, double*)" "GHealpix::pix2ang_nest(int, double*, double*)"
44#define G_ORDERING_SET"GHealpix::ordering(std::string&)" "GHealpix::ordering(std::string&)"
45
46/* __ Macros _____________________________________________________________ */
47
48/* __ Coding definitions _________________________________________________ */
49
50/* __ Debug definitions __________________________________________________ */
51
52/* __ Local prototypes ___________________________________________________ */
53
54/* __ Constants __________________________________________________________ */
55const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
56const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
57const int order_max = 13;
58const int ns_max = 1 << order_max;
59
60/* __ Static conversion arrays ___________________________________________ */
61static short ctab[0x100];
62static short utab[0x100];
63
64/* __ Globals ____________________________________________________________ */
65
66
67/*==========================================================================
68 = =
69 = Constructors/destructors =
70 = =
71 ==========================================================================*/
72
73/***********************************************************************//**
74 * @brief Void constructor
75 ***************************************************************************/
76GHealpix::GHealpix(void) : GSkyProjection()
77{
78 // Initialise class members
79 init_members();
80
81 // Return
82 return;
83}
84
85
86/***********************************************************************//**
87 * @brief Constructor
88 *
89 * @param[in] nside Number of sides.
90 * @param[in] order Pixel ordering ('RING' or 'NESTED').
91 * @param[in] coords Coordinate system ('EQU' or 'GAL').
92 *
93 * @exception GException::wcs_hpx_bad_nside
94 * Invalid nside parameter.
95 * @exception GException::wcs_bad_coords
96 * Invalid coordsys parameter.
97 * @exception GException::wcs_hpx_bad_ordering
98 * Invalid ordering parameter.
99 ***************************************************************************/
100GHealpix::GHealpix(const int& nside,
101 const std::string& order,
102 const std::string& coords) : GSkyProjection()
103{
104 // Initialise class members
105 init_members();
106
107 // Check nside parameter (power of 2 between 1 and 8192)
108 if (nside2order(nside) == -1) {
109 throw GException::wcs_hpx_bad_nside(G_CONSTRUCT"GHealpix::GHealpix(int& ,std::string& ,std::string&)", nside);
110 }
111
112 // Set coordinate system
113 coordsys(coords);
114
115 // Set pixel ordering
116 ordering(order);
117
118 // Set Healpix parameters
119 m_nside = nside;
120 m_npface = m_nside * m_nside;
121 m_ncap = 2 * (m_npface - m_nside);
122 m_num_pixels = 12 * m_npface;
123 m_fact2 = 4.0 / m_num_pixels;
124 m_fact1 = 2 * m_nside * m_fact2;
125 m_omega = gammalib::fourpi / m_num_pixels;
126 m_order = nside2order(m_nside);
127
128 // Return
129 return;
130}
131
132
133/***********************************************************************//**
134 * @brief Constructor from FITS HDU table
135 *
136 * @param[in] hdu FITS HDU.
137 ***************************************************************************/
138GHealpix::GHealpix(const GFitsHDU& hdu) : GSkyProjection()
139{
140 // Initialise class members
141 init_members();
142
143 // Read Healpix defintion from FITS table
144 read(hdu);
145
146 // Return
147 return;
148}
149
150
151/***********************************************************************//**
152 * @brief Copy constructor
153 *
154 * @param[in] proj Healpix projection.
155 ***************************************************************************/
156GHealpix::GHealpix(const GHealpix& proj) : GSkyProjection(proj)
157{
158 // Initialise class members for clean destruction
159 init_members();
160
161 // Copy members
162 copy_members(proj);
163
164 // Return
165 return;
166}
167
168
169/***********************************************************************//**
170 * @brief Destructor
171 ***************************************************************************/
172GHealpix::~GHealpix()
173{
174 // Free members
175 free_members();
176
177 // Return
178 return;
179}
180
181
182/*==========================================================================
183 = =
184 = Operators =
185 = =
186 ==========================================================================*/
187
188/***********************************************************************//**
189 * @brief Assignment operator
190 *
191 * @param[in] proj Healpix projection.
192 * @return Healpix projection.
193 ***************************************************************************/
194GHealpix& GHealpix::operator=(const GHealpix& proj)
195{
196 // Execute only if object is not identical
197 if (this != &proj) {
198
199 // Copy base class members
200 this->GSkyProjection::operator=(proj);
201
202 // Free members
203 free_members();
204
205 // Initialise private members for clean destruction
206 init_members();
207
208 // Copy members
209 copy_members(proj);
210
211 } // endif: object was not identical
212
213 // Return this object
214 return *this;
215}
216
217
218/*==========================================================================
219 = =
220 = Public methods =
221 = =
222 ==========================================================================*/
223
224/***********************************************************************//**
225 * @brief Clear object
226 *
227 * This method properly resets the object to an initial state.
228 ***************************************************************************/
229void GHealpix::clear(void)
230{
231 // Free class members (base and derived classes, derived class first)
232 free_members();
233 this->GSkyProjection::free_members();
234
235 // Initialise members
236 this->GSkyProjection::init_members();
237 init_members();
238
239 // Return
240 return;
241}
242
243
244/***********************************************************************//**
245 * @brief Clone instance
246 ***************************************************************************/
247GHealpix* GHealpix::clone(void) const
248{
249 return new GHealpix(*this);
250}
251
252
253/***********************************************************************//**
254 * @brief Read Healpix definiton from FITS header
255 *
256 * @param[in] hdu FITS HDU.
257 *
258 * @exception GException::wcs
259 * Unable to load Healpix definition from HDU.
260 * @exception GException::wcs_bad_coords
261 * Invalid coordsys parameter.
262 * @exception GException::wcs_hpx_bad_ordering
263 * Invalid ordering parameter.
264 ***************************************************************************/
265void GHealpix::read(const GFitsHDU& hdu)
266{
267 // Clear object
268 clear();
269
270 // Check if we have a healpix representation
271 if (hdu.string("PIXTYPE") != "HEALPIX") {
272 throw GException::wcs(G_READ"GHealpix::read(GFitsHDU&)", "HDU does not contain Healpix data");
273 }
274
275 // Get pixel ordering
276 std::string order = hdu.string("ORDERING");
277
278 // Get coordinate system.
279 // First search for HIER_CRD keyword (this has been used in older
280 // versions of LAT exposure cubes). If not found then search for standard
281 // COORDSYS keyword.
282 std::string coords;
283 if (hdu.hascard("HIER_CRD")) {
284 coords = hdu.string("HIER_CRD");
285 }
286 else if (hdu.hascard("COORDSYS")) {
287 coords = hdu.string("COORDSYS");
288 }
289
290 // Set coordinate system
291 coordsys(coords);
292
293 // Set pixel ordering
294 ordering(order);
295
296 // Get Healpix resolution and determine number of pixels and solid angle
297 m_nside = hdu.integer("NSIDE");
298 m_npface = m_nside * m_nside;
299 m_ncap = 2 * (m_npface - m_nside);
300 m_num_pixels = 12 * m_npface;
301 m_fact2 = 4.0 / m_num_pixels;
302 m_fact1 = 2 * m_nside * m_fact2;
303 m_omega = gammalib::fourpi / m_num_pixels;
304 m_order = nside2order(m_nside);
305
306 // Return
307 return;
308}
309
310
311/***********************************************************************//**
312 * @brief Write Healpix definiton into FITS HDU
313 *
314 * @param[in] hdu FITS HDU.
315 *
316 * Writes the following keywords in the FITS HDU:
317 * EXTNAME = HEALPIX
318 * PIXTYPE = HEALPIX
319 * NSIDE = nside() = m_nside
320 * FIRSTPIX = 0
321 * LASTPIX = npix()-1 = m_num_pixels-1
322 * ORDERING = ordering()
323 * COORDSYS = coordsys()
324 ***************************************************************************/
325void GHealpix::write(GFitsHDU& hdu) const
326{
327 // Set extension name
328 hdu.extname("HEALPIX");
329
330 // Set keywords
331 hdu.card("PIXTYPE", "HEALPIX", "HEALPix pixelisation");
332 hdu.card("NSIDE", nside(), "HEALPix resolution parameter");
333 hdu.card("FIRSTPIX", 0, "Index of first pixel");
334 hdu.card("LASTPIX", npix()-1, "Index of last pixel");
335 hdu.card("ORDERING", ordering(),
336 "Pixel ordering scheme, either RING or NESTED");
337 hdu.card("COORDSYS", coordsys(),
338 "Coordinate system, either EQU or GAL");
339
340 // Return
341 return;
342}
343
344
345/***********************************************************************//**
346 * @brief Returns sky direction of pixel
347 *
348 * @param[in] pixel Sky map pixel.
349 ***************************************************************************/
350GSkyDir GHealpix::pix2dir(const GSkyPixel& pixel) const
351{
352 // Throw an exception if sky map pixel is not 1D
353 if (!pixel.is1D()) {
354 std::string msg = "Sky map pixel "+pixel.print()+" is not"
355 " 1-dimensional.\n"
356 "Only 1-dimensional pixels are supported by the"
357 " Healpix projection.";
358 throw GException::invalid_argument(G_XY2DIR"GHealpix::xy2dir(GSkyPixel&)", msg);
359 }
360
361 // Perform ordering dependent conversion
362 double theta = 0.0;
363 double phi = 0.0;
364 switch (m_ordering) {
365 case 0:
366 pix2ang_ring(int(pixel), &theta, &phi);
367 break;
368 case 1:
369 pix2ang_nest(int(pixel), &theta, &phi);
370 break;
371 default:
372 break;
373 }
374
375 // Store coordinate system dependent result
376 GSkyDir result;
377 switch (m_coordsys) {
378 case 0:
379 result.radec(phi, gammalib::pihalf - theta);
380 break;
381 case 1:
382 result.lb(phi, gammalib::pihalf - theta);
383 break;
384 default:
385 break;
386 }
387
388 // Return
389 return result;
390}
391
392
393/***********************************************************************//**
394 * @brief Returns sky map pixel of sky coordinate
395 *
396 * @param[in] dir Sky coordinate.
397 * @return Sky map pixel.
398 ***************************************************************************/
399GSkyPixel GHealpix::dir2pix(const GSkyDir& dir) const
400{
401 // Compute coordinate system dependent (z,phi)
402 double z;
1
'z' declared without an initial value
→
403 double phi;
404 switch (m_coordsys) {
2
←
Control jumps to the 'default' case at line 413
→
405 case 0:
406 z = cos(gammalib::pihalf - dir.dec());
407 phi = dir.ra();
408 break;
409 case 1:
410 z = cos(gammalib::pihalf - dir.b());
411 phi = dir.l();
412 break;
413 default:
414 break;
3
←
 Execution continues on line 418
→
415 }
416
417 // Perform ordering dependent conversion
418 int index;
419 switch (m_ordering) {
4
←
Control jumps to 'case 0:' at line 420
→
420 case 0:
421 index = ang2pix_z_phi_ring(z, phi);
5
←
Function call argument is an uninitialized value
422 break;
423 case 1:
424 index = ang2pix_z_phi_nest(z, phi);
425 break;
426 default:
427 break;
428 }
429
430 // Return sky map pixel
431 return (GSkyPixel(index));
432}
433
434
435/***********************************************************************//**
436 * @brief Returns ordering parameter.
437 ***************************************************************************/
438std::string GHealpix::ordering(void) const
439{
440 // Set pixel ordering type
441 std::string s_ordering;
442 switch (m_ordering) {
443 case 0:
444 s_ordering = "RING";
445 break;
446 case 1:
447 s_ordering = "NESTED";
448 break;
449 default:
450 s_ordering = "UNKNOWN";
451 break;
452 }
453
454 // Return ordering
455 return s_ordering;
456}
457
458
459/***********************************************************************//**
460 * @brief Set pixel ordering.
461 *
462 * @param[in] ordering Pixel ordering (RING or NEST/NESTED).
463 *
464 * @exception GException::wcs_hpx_bad_ordering
465 * Invalid ordering parameter.
466 ***************************************************************************/
467void GHealpix::ordering(const std::string& ordering)
468{
469 // Convert argument to upper case
470 std::string uordering = gammalib::toupper(ordering);
471
472 // Set pixel ordering
473 if (uordering == "RING") {
474 m_ordering = 0;
475 }
476 else if (uordering == "NESTED" || uordering == "NEST") {
477 m_ordering = 1;
478 }
479 else {
480 throw GException::wcs_hpx_bad_ordering(G_ORDERING_SET"GHealpix::ordering(std::string&)", ordering);
481 }
482
483 // Return
484 return;
485}
486
487
488/***********************************************************************//**
489 * @brief Print WCS information
490 *
491 * @param[in] chatter Chattiness (defaults to NORMAL).
492 * @return String containing WCS information.
493 ***************************************************************************/
494std::string GHealpix::print(const GChatter& chatter) const
495{
496 // Initialise result string
497 std::string result;
498
499 // Continue only if chatter is not silent
500 if (chatter != SILENT) {
501
502 // Append header
503 result.append("=== GHealpix ===");
504
505 // Append information
506 result.append("\n"+gammalib::parformat("Coordinate system"));
507 result.append(coordsys());
508 result.append("\n"+gammalib::parformat("Nside (# of divisions)"));
509 result.append(gammalib::str(m_nside));
510 result.append("\n"+gammalib::parformat("Npface (pixels per face)"));
511 result.append(gammalib::str(m_npface));
512 result.append("\n"+gammalib::parformat("Ncap (# of cap pixels)"));
513 result.append(gammalib::str(m_ncap));
514 result.append("\n"+gammalib::parformat("Npix (# of pixels)"));
515 result.append(gammalib::str(m_num_pixels));
516 result.append("\n"+gammalib::parformat("Order"));
517 result.append(gammalib::str(m_order));
518 result.append("\n"+gammalib::parformat("Solid angle per pixel"));
519 result.append(gammalib::str(m_omega)+" sr");
520 result.append("\n"+gammalib::parformat("Ordering")+ordering());
521
522 } // endif: chatter was not silent
523
524 // Return result
525 return result;
526}
527
528
529/*==========================================================================
530 = =
531 = Private methods =
532 = =
533 ==========================================================================*/
534
535/***********************************************************************//**
536 * @brief Initialise class members
537 ***************************************************************************/
538void GHealpix::init_members(void)
539{
540 // Initialise members
541 //m_type = "HPX";
542 m_nside = 0;
543 m_npface = 0;
544 m_ncap = 0;
545 m_order = 0;
546 m_ordering = 0;
547 m_num_pixels = 0;
548 m_fact1 = 0.0;
549 m_fact2 = 0.0;
550 m_omega = 0.0;
551
552 // Construct conversion arrays
553 for (int m = 0; m < 0x100; ++m) {
554 ctab[m] =
555 (m&0x1 ) | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
556 | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4);
557 utab[m] =
558 (m&0x1 ) | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
559 | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7);
560 }
561
562 // Return
563 return;
564}
565
566
567/***********************************************************************//**
568 * @brief Copy class members
569 *
570 * @param[in] proj Healpix projection.
571 ***************************************************************************/
572void GHealpix::copy_members(const GHealpix& proj)
573{
574 // Copy attributes
575 m_nside = proj.m_nside;
576 m_npface = proj.m_npface;
577 m_ncap = proj.m_ncap;
578 m_order = proj.m_order;
579 m_ordering = proj.m_ordering;
580 m_num_pixels = proj.m_num_pixels;
581 m_fact1 = proj.m_fact1;
582 m_fact2 = proj.m_fact2;
583 m_omega = proj.m_omega;
584
585 // Return
586 return;
587}
588
589
590/***********************************************************************//**
591 * @brief Delete class members
592 ***************************************************************************/
593void GHealpix::free_members(void)
594{
595 // Return
596 return;
597}
598
599
600/***********************************************************************//**
601 * @brief Returns true if argument is identical
602 *
603 * @param[in] proj Sky projection.
604 *
605 * This method is a helper for the sky projection friends.
606 ***************************************************************************/
607bool GHealpix::compare(const GSkyProjection& proj) const
608{
609 // Initialise result
610 bool result = false;
611
612 // Continue only we compare to a GHealpix object
613 const GHealpix* ptr = dynamic_cast<const GHealpix*>(&proj);
614 if (ptr != NULL__null) {
615 result = ((m_coordsys == ptr->m_coordsys) &&
616 (m_nside == ptr->m_nside) &&
617 (m_npface == ptr->m_npface) &&
618 (m_ncap == ptr->m_ncap) &&
619 (m_order == ptr->m_order) &&
620 (m_ordering == ptr->m_ordering) &&
621 (m_num_pixels == ptr->m_num_pixels));
622 }
623
624 // Return result
625 return result;
626}
627
628
629/***********************************************************************//**
630 * @brief Convert nside to order
631 *
632 * @param[in] nside Number of sides.
633 ***************************************************************************/
634int GHealpix::nside2order(int nside)
635{
636 // Initialise order
637 int order = -1;
638
639 // Determine order
640 for (int m = 0; m <= order_max; ++m) {
641 int nstest = 1 << m;
642 if (nside == nstest) {
643 order = m;
644 break;
645 }
646 if (nside < nstest)
647 break;
648 }
649
650 // Return order
651 return order;
652}
653
654
655/***********************************************************************//**
656 * @brief Convert pixel index to (x,y) coordinate
657 *
658 * @param[in] ipix Pixel index for which (x,y) are to be computed.
659 * @param[out] x Pointer to x coordinate.
660 * @param[out] y Pointer to y coordinate.
661 ***************************************************************************/
662void GHealpix::pix2xy(const int& ipix, int* x, int* y) const
663{
664 // Set x coordinate
665 int raw = (ipix & 0x5555) | ((ipix & 0x55550000) >> 15);
666 *x = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
667
668 // Set y coordinate
669 raw = ((ipix & 0xaaaa) >> 1) | ((ipix & 0xaaaa0000) >> 16);
670 *y = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
671
672 // Return
673 return;
674}
675
676
677/***********************************************************************//**
678 * @brief Convert (x,y) coordinate to pixel
679 *
680 * @param[in] x x coordinate.
681 * @param[in] y y coordinate.
682 ***************************************************************************/
683int GHealpix::xy2pix(int x, int y) const
684{
685 // Return pixel
686 return utab[x&0xff] | (utab[x>>8]<<16) | (utab[y&0xff]<<1) | (utab[y>>8]<<17);
687}
688
689
690/***********************************************************************//**
691 * @brief Convert pixel index to (theta,phi) angles for ring ordering
692 *
693 * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
694 * @param[out] theta Pointer to result zenith angle in radians.
695 * @param[out] phi Pointer to result azimuth angle in radians.
696 *
697 * @exception GException::out_of_range
698 * Pixel index is out of range.
699 ***************************************************************************/
700void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi) const
701{
702 // Check if ipix is in range
703 if (ipix < 0 || ipix >= m_num_pixels)
704 throw GException::out_of_range(G_PIX2ANG_RING"GHealpix::pix2ang_ring(int, double*, double*)", ipix, 0, m_num_pixels-1);
705
706 // Handle North Polar cap
707 if (ipix < m_ncap) {
708 int iring = int(0.5*(1+isqrt(1+2*ipix))); // counted from North pole
709 int iphi = (ipix+1) - 2*iring*(iring-1);
710 *theta = acos(1.0 - (iring*iring) * m_fact2);
711 *phi = (iphi - 0.5) * gammalib::pi/(2.0*iring);
712 }
713
714 // Handle Equatorial region
715 else if (ipix < (m_num_pixels - m_ncap)) {
716 int ip = ipix - m_ncap;
717 int iring = ip/(4*m_nside) + m_nside; // counted from North pole
718 int iphi = ip%(4*m_nside) + 1;
719 double fodd = ((iring+m_nside)&1) ? 1 : 0.5;
720 int nl2 = 2*m_nside;
721 *theta = std::acos((nl2 - iring) * m_fact1);
722 *phi = (iphi - fodd) * gammalib::pi/nl2;
723 }
724
725 // Handle South Polar cap
726 else {
727 int ip = m_num_pixels - ipix;
728 int iring = int(0.5*(1+isqrt(2*ip-1))); // Counted from South pole
729 int iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
730 *theta = std::acos(-1.0 + (iring*iring) * m_fact2);
731 *phi = (iphi - 0.5) * gammalib::pi/(2.*iring);
732 }
733
734 // Return
735 return;
736}
737
738
739/***********************************************************************//**
740 * @brief Convert pixel index to (theta,phi) angles for nested ordering
741 *
742 * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
743 * @param[out] theta Pointer to result zenith angle in radians.
744 * @param[out] phi Pointer to result azimuth angle in radians.
745 *
746 * @exception GException::out_of_range
747 * Pixel index is out of range.
748 ***************************************************************************/
749void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) const
750{
751 // Check if ipix is in range
752 if (ipix < 0 || ipix >= m_num_pixels)
753 throw GException::out_of_range(G_PIX2ANG_NEST"GHealpix::pix2ang_nest(int, double*, double*)", ipix, 0, m_num_pixels-1);
754
755 // Get face number and index in face
756 int nl4 = 4 * m_nside;
757 int face_num = ipix >> (2*m_order); // Face number in {0,11}
758 int ipf = ipix & (m_npface - 1);
759
760 // Get pixel coordinates
761 int ix;
762 int iy;
763 pix2xy(ipf, &ix, &iy);
764
765 // Computes the z coordinate on the sphere
766 int jr = (jrll[face_num] << m_order) - ix - iy - 1;
767
768 // Declare result variables
769 int nr;
770 double z;
771 int kshift;
772
773 // North pole region
774 if (jr < m_nside) {
775 nr = jr;
776 z = 1. - nr*nr*m_fact2;
777 kshift = 0;
778 }
779
780 // South pole region
781 else if (jr > 3*m_nside) {
782 nr = nl4 - jr;
783 z = nr*nr*m_fact2 - 1;
784 kshift = 0;
785 }
786
787 // Equatorial region
788 else {
789 nr = m_nside;
790 z = (2*m_nside-jr) * m_fact1;
791 kshift = (jr-m_nside) & 1;
792 }
793
794 // Computes the phi coordinate on the sphere, in [0,2Pi]
795 int jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
796 if (jp > nl4) jp -= nl4;
797 if (jp < 1) jp += nl4;
798
799 // Computes Theta and Phi
800 *theta = acos(z);
801 *phi = (jp - (kshift+1)*0.5) * (gammalib::pihalf / nr);
802
803 // Return
804 return;
805}
806
807
808/***********************************************************************//**
809 * @brief Returns pixels which contains angular coordinates (z,phi)
810 *
811 * @param[in] z Cosine of zenith angle - cos(theta).
812 * @param[in] phi Azimuth angle in radians.
813 ***************************************************************************/
814int GHealpix::ang2pix_z_phi_ring(double z, double phi) const
815{
816 // Initialise pixel
817 int ipix = 0;
818
819 // Setup
820 double za = fabs(z);
821 double tt = gammalib::modulo(phi, gammalib::twopi) * gammalib::inv_pihalf; // in [0,4)
822
823 // Equatorial region
824 if (za <= gammalib::twothird) {
825 double temp1 = m_nside*(0.5+tt);
826 double temp2 = m_nside*z*0.75;
827 int jp = int(temp1-temp2); // index of ascending edge line
828 int jm = int(temp1+temp2); // index of descending edge line
829 int ir = m_nside + 1 + jp - jm; // in {1,2n+1}
830 int kshift = 1 - (ir & 1); // kshift=1 if ir even, 0 otherwise
831 int ip = (jp+jm-m_nside+kshift+1)/2; // in {0,4n-1}
832 ip = int(gammalib::modulo(ip,4*m_nside));
833 ipix = m_ncap + (ir-1)*4*m_nside + ip;
834 }
835
836 // North & South polar caps
837 else {
838 double tp = tt - int(tt);
839 double tmp = m_nside * std::sqrt(3*(1-za));
840 int jp = int(tp*tmp); // increasing edge line index
841 int jm = int((1.0-tp)*tmp); // decreasing edge line index
842 int ir = jp + jm + 1; // ring number counted from the closest pole
843 int ip = int(tt*ir); // in {0,4*ir-1}
844 ip = int(gammalib::modulo(ip,4*ir));
845 if (z>0)
846 ipix = 2*ir*(ir-1) + ip;
847 else
848 ipix = m_num_pixels - 2*ir*(ir+1) + ip;
849 }
850
851 // Return pixel
852 return ipix;
853}
854
855
856/***********************************************************************//**
857 * @brief Returns pixels which contains angular coordinates (z,phi)
858 *
859 * @param[in] z Cosine of zenith angle - cos(theta).
860 * @param[in] phi Azimuth angle in radians.
861 ***************************************************************************/
862int GHealpix::ang2pix_z_phi_nest(double z, double phi) const
863{
864 // Initialise face and pixel numbers
865 int face_num;
866 int ix;
867 int iy;
868
869 // Setup
870 double za = fabs(z);
871 double tt = gammalib::modulo(phi, gammalib::twopi) * gammalib::inv_pihalf; // in [0,4)
872
873 // Equatorial region
874 if (za <= gammalib::twothird) {
875 double temp1 = ns_max*(0.5+tt);
876 double temp2 = ns_max*z*0.75;
877 int jp = int(temp1-temp2); // index of ascending edge line
878 int jm = int(temp1+temp2); // index of descending edge line
879 int ifp = jp >> order_max; // in {0,4}
880 int ifm = jm >> order_max;
881 if (ifp == ifm) // faces 4 to 7
882 face_num = (ifp==4) ? 4: ifp+4;
883 else if (ifp < ifm) // (half-)faces 0 to 3
884 face_num = ifp;
885 else // (half-)faces 8 to 11
886 face_num = ifm + 8;
887 ix = jm & (ns_max-1);
888 iy = ns_max - (jp & (ns_max-1)) - 1;
889 }
890
891 // Polar region, za > 2/3
892 else {
893 int ntt = int(tt);
894 double tp = tt-ntt;
895 double tmp = ns_max * std::sqrt(3*(1-za));
896 int jp = int(tp*tmp); // increasing edge line index
897 int jm = int((1.0-tp)*tmp); // decreasing edge line index
898 if (jp >= ns_max) jp = ns_max-1; // for points too close to the boundary
899 if (jm >= ns_max) jm = ns_max-1;
900 if (z >= 0) {
901 face_num = ntt; // in {0,3}
902 ix = ns_max - jm - 1;
903 iy = ns_max - jp - 1;
904 }
905 else {
906 face_num = ntt + 8; // in {8,11}
907 ix = jp;
908 iy = jm;
909 }
910 }
911
912 // Get pixel
913 int ipf = xy2pix(ix, iy);
914 ipf >>= (2*(order_max - m_order)); // in {0, nside**2 - 1}
915 return ipf + (face_num<<(2*m_order)); // in {0, 12*nside**2 - 1}
916}
917
918
919/***********************************************************************//**
920 * @brief Integer n that fulfills n*n <= arg < (n+1)*(n+1)
921 *
922 * @param[in] arg Argument.
923 *
924 * Returns the integer \a n, which fulfills \a n*n <= arg < (n+1)*(n+1).
925 ***************************************************************************/
926unsigned int GHealpix::isqrt(unsigned int arg) const
927{
928 // Return
929 return unsigned(std::sqrt(arg+0.5));
930}