Peter Deiml / gammalib

Commit fdfb897d4246eadfef1a1ef56d5b4e0fef4daaa8

Authored by knodlseder
1 parent dea0af5f

enhance GHealpix and GSkyDir classes

Inline Side-by-side
Showing 12 changed files with 749 additions and 130 deletions
include/GException.hpp
  View file @fdfb897
... ... @@ -233,6 +233,12 @@ public:
233 233 rsp_invalid_type(string origin, string type);
234 234 };
235 235  
  236 + // General Healpix error
  237 + class healpix : public GExceptionHandler {
  238 + public:
  239 + healpix(string origin, string message);
  240 + };
  241 +
236 242 };
237 243  
238 244 #endif /* GEXCEPTION_HPP */
... ...
include/GFitsHDU.hpp
  View file @fdfb897
... ... @@ -85,7 +85,7 @@ private:
85 85  
86 86 // Private data area
87 87 __fitsfile m_fitsfile; //!< FITS file pointer
88   - int m_hdunum; //!< HDU number (starting from 1)
  88 + int m_hdunum; //!< HDU number (starting from 0)
89 89 std::string m_name; //!< HDU name (extname)
90 90 int m_type; //!< HDU type
91 91 GFitsHeader* m_header; //!< HDU header
... ...
include/GHealpix.hpp
  View file @fdfb897
... ... @@ -21,6 +21,7 @@
21 21  
22 22 /* __ Includes ___________________________________________________________ */
23 23 #include "GFitsHDU.hpp"
  24 +#include "GSkyDir.hpp"
24 25  
25 26 /* __ Namespaces _________________________________________________________ */
26 27  
... ... @@ -32,6 +33,9 @@
32 33 ***************************************************************************/
33 34 class GHealpix {
34 35  
  36 + // I/O friends
  37 + friend ostream& operator<< (ostream& os, const GHealpix& pixels);
  38 +
35 39 public:
36 40 // Constructors and destructors
37 41 GHealpix();
... ... @@ -42,10 +46,11 @@ public:
42 46 GHealpix& operator= (const GHealpix& pixels);
43 47  
44 48 // Methods
45   - void load(const GFitsHDU* hdu);
46   - int nside(void) const;
47   - int num_pixels(void) const;
48   - double omega(void) const;
  49 + void load(const GFitsHDU* hdu);
  50 + int nside(void) const;
  51 + int num_pixels(void) const;
  52 + double omega(void) const;
  53 + GSkyDir pix2ang(const int& ipix);
49 54  
50 55 private:
51 56 // Private methods
... ... @@ -53,16 +58,19 @@ private:
53 58 void copy_members(const GHealpix& pixels);
54 59 void free_members(void);
55 60 GHealpix* clone(void) const;
56   -
  61 + void mk_pix2xy(void);
  62 + void pix2ang_ring(int ipix, double* theta, double* phi);
  63 + void pix2ang_nest(int ipix, double* theta, double* phi);
  64 +
57 65 // Private data area
58   - int m_nside; //!< Number of divisions of each base pixel (1,2,..)
59   - int m_order; //!< Ordering (0=ring, 1=nested)
60   - int m_coordsys; //!< Coordinate system (0=equatorial, 1=galactic)
61   - int m_num_pixels; //!< Number of pixels
62   - int m_size_pixels; //!< Vector size of each pixel
63   - double* m_pixels; //!< Pixel array
64   - double* m_lon; //!< RA or GLON values for each pixel
65   - double* m_lat; //!< DEC or GLAT values for each pixel
  66 + int m_nside; //!< Number of divisions of each base pixel (1-8192)
  67 + int m_order; //!< Ordering (0=ring, 1=nested, -1=?)
  68 + int m_coordsys; //!< Coordinate system (0=equatorial, 1=galactic, -1=?)
  69 + int m_num_pixels; //!< Number of pixels
  70 + int m_size_pixels; //!< Vector size of each pixel
  71 + double m_omega; //!< Solid angle of pixel
  72 + double* m_pixels; //!< Pixel array
  73 + GSkyDir* m_dir; //!< Sky direction for each pixel
66 74 };
67 75  
68 76 #endif /* GHEALPIX_HPP */
... ...
include/GSkyDir.hpp
  View file @fdfb897
... ... @@ -24,54 +24,48 @@
24 24 /* __ Namespaces _________________________________________________________ */
25 25  
26 26  
27   -/***************************************************************************
28   - * GSkyDir class definition *
29   - ***************************************************************************/
30   -/**
  27 +/***********************************************************************//**
31 28 * @class GSkyDir
32 29 *
33   - * @brief Interface for the sky direction classes.
34   - *
35   - * @author J. Knodlseder
36   - */
  30 + * @brief GSkyDir class interface defintion
  31 + ***************************************************************************/
37 32 class GSkyDir {
38 33  
  34 + // I/O friends
  35 + friend ostream& operator<< (ostream& os, const GSkyDir& dir);
  36 +
39 37 public:
40   - /// Constructor
  38 + // Constructors and destructors
41 39 GSkyDir();
42   -
43   - /// Copy constructor
44 40 GSkyDir(const GSkyDir& dir);
45   -
46   - /// Destructor
47 41 virtual ~GSkyDir();
48 42  
49   - /// Assignment operator
  43 + // Operators
50 44 GSkyDir& operator= (const GSkyDir& dir);
51 45  
52   -protected:
53   - /// Longitude
54   - double m_lon;
55   -
56   - /// Latitude
57   - double m_lat;
  46 + // Methods
  47 + void radec(const double& ra, const double& dec);
  48 + void lb(const double& l, const double& b);
  49 + double l(void);
  50 + double b(void);
  51 + double ra(void);
  52 + double dec(void);
58 53  
59   - /// Protected method that initialises protected members of the class
  54 +private:
  55 + // Private methods
60 56 void init_members(void);
61   -
62   - /// Protected method that copies protected members of the class
63 57 void copy_members(const GSkyDir& dir);
64   -
65   - /// Protected method that frees memory allocated by instances of the class
66 58 void free_members(void);
67   -
68   -private:
69   -
  59 + void equ2gal(void);
  60 + void gal2equ(void);
  61 +
  62 + // Private data area
  63 + int m_has_lb; //!< Has galactic coordinates
  64 + int m_has_radec; //!< Has equatorial coordinates
  65 + double m_l; //!< Galactic longitude in radians
  66 + double m_b; //!< Galactic latitude in radians
  67 + double m_ra; //!< Right Ascension in radians
  68 + double m_dec; //!< Declination in radians
70 69 };
71 70  
72   -
73   -/***************************************************************************
74   - * Inline methods *
75   - ***************************************************************************/
76   -
77 71 #endif /* GSKYDIR_HPP */
... ...
pyext/GHealpix.i 0 → 100644
  View file @fdfb897
  1 +/***************************************************************************
  2 + * GHealpix.i - Healpix sky representation class SWIG definition *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright : (C) 2008 by Jurgen Knodlseder *
  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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GHealpix.i
  15 + * @brief GHealpix class SWIG file.
  16 + * @author J. Knodlseder
  17 + */
  18 +%{
  19 +/* Put headers and other declarations here that are needed for compilation */
  20 +#include "GHealpix.hpp"
  21 +%}
  22 +
  23 +
  24 +/***********************************************************************//**
  25 + * @class GHealpix
  26 + *
  27 + * @brief GHealpix class interface defintion
  28 + ***************************************************************************/
  29 +class GHealpix {
  30 +public:
  31 + // Constructors and destructors
  32 + GHealpix();
  33 + GHealpix(const GHealpix& pixels);
  34 + virtual ~GHealpix();
  35 +
  36 + // Methods
  37 + void load(const GFitsHDU* hdu);
  38 + int nside(void) const;
  39 + int num_pixels(void) const;
  40 + double omega(void) const;
  41 + GSkyDir pix2ang(const int& ipix);
  42 +};
  43 +
  44 +
  45 +/***********************************************************************//**
  46 + * @brief GHealpix class extension
  47 + ***************************************************************************/
  48 +%extend GHealpix {
  49 + char *__str__() {
  50 + static char str_buffer[1001];
  51 + std::ostringstream buffer;
  52 + buffer << *self;
  53 + std::string str = buffer.str();
  54 + strncpy(str_buffer, (char*)str.c_str(), 1001);
  55 + str_buffer[1000] = '\0';
  56 + return str_buffer;
  57 + }
  58 + GHealpix copy() {
  59 + return (*self);
  60 + }
  61 +};
... ...
pyext/GLATResponse.i
  View file @fdfb897
1 1 /***************************************************************************
2   - * GLATResponse.hpp - GLAST LAT Response class SWIG definition *
  2 + * GLATResponse.i - GLAST LAT Response class SWIG definition *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright : (C) 2008 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
... ...
pyext/GSkyDir.i 0 → 100644
  View file @fdfb897
  1 +/***************************************************************************
  2 + * GSkyDir.i - Sky direction class SWIG file *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright : (C) 2008 by Jurgen Knodlseder *
  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 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GSkyDir.i
  15 + * @brief GSkyDir class SWIG file.
  16 + * @author J. Knodlseder
  17 + */
  18 +%{
  19 +/* Put headers and other declarations here that are needed for compilation */
  20 +#include "GSkyDir.hpp"
  21 +%}
  22 +
  23 +
  24 +/***********************************************************************//**
  25 + * @class GSkyDir
  26 + *
  27 + * @brief GSkyDir class interface defintion
  28 + ***************************************************************************/
  29 +class GSkyDir {
  30 +public:
  31 + // Constructors and destructors
  32 + GSkyDir();
  33 + GSkyDir(const GSkyDir& dir);
  34 + virtual ~GSkyDir();
  35 +
  36 + // Methods
  37 + void radec(const double& ra, const double& dec);
  38 + void lb(const double& l, const double& b);
  39 + double l(void);
  40 + double b(void);
  41 + double ra(void);
  42 + double dec(void);
  43 +};
  44 +
  45 +
  46 +/***********************************************************************//**
  47 + * @brief GSkyDir class extension
  48 + ***************************************************************************/
  49 +%extend GSkyDir {
  50 + char *__str__() {
  51 + static char str_buffer[1001];
  52 + std::ostringstream buffer;
  53 + buffer << *self;
  54 + std::string str = buffer.str();
  55 + strncpy(str_buffer, (char*)str.c_str(), 1001);
  56 + str_buffer[1000] = '\0';
  57 + return str_buffer;
  58 + }
  59 + GSkyDir copy() {
  60 + return (*self);
  61 + }
  62 +};
... ...
pyext/gammalib.i
  View file @fdfb897
... ... @@ -18,6 +18,8 @@
18 18 /* General */
19 19 #include "GException.i"
20 20 #include "GNodeArray.i"
  21 +#include "GHealpix.i"
  22 +#include "GSkyDir.i"
21 23  
22 24 /* Response classes */
23 25 #include "GResponse.i"
... ...
src/GException.cpp
  View file @fdfb897
... ... @@ -519,3 +519,16 @@ GException::rsp_invalid_type::rsp_invalid_type(string origin,
519 519 m_origin = origin;
520 520 m_message = "Invalid response type '"+type+"' specified";
521 521 }
  522 +
  523 +
  524 +/***********************************************************************//**
  525 + * @brief General Healpix error
  526 + *
  527 + * @param[in] origin Method that throws the error.
  528 + * @param[in] message Error message.
  529 + ***************************************************************************/
  530 +GException::healpix::healpix(string origin, string message)
  531 +{
  532 + m_origin = origin;
  533 + m_message = message;
  534 +}
... ...
src/GHealpix.cpp
  View file @fdfb897
... ... @@ -18,6 +18,8 @@
18 18 */
19 19  
20 20 /* __ Includes ___________________________________________________________ */
  21 +#include <iostream>
  22 +#include <math.h>
21 23 #include "GException.hpp"
22 24 #include "GTools.hpp"
23 25 #include "GHealpix.hpp"
... ... @@ -25,6 +27,9 @@
25 27 /* __ Namespaces _________________________________________________________ */
26 28  
27 29 /* __ Method name definitions ____________________________________________ */
  30 +#define G_LOAD "GHealpix::load(const GFitsHDU*)"
  31 +#define G_PIX2ANG_RING "GHealpix::pix2ang_ring(int,double*,double*)"
  32 +#define G_PIX2ANG_NEST "GHealpix::pix2ang_nest(int,double*,double*)"
28 33  
29 34 /* __ Macros _____________________________________________________________ */
30 35  
... ... @@ -32,6 +37,13 @@
32 37  
33 38 /* __ Debug definitions __________________________________________________ */
34 39  
  40 +/* __ Constants __________________________________________________________ */
  41 +const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
  42 +const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
  43 +
  44 +/* __ Static conversion arrays ___________________________________________ */
  45 +static int pix2x[1024];
  46 +static int pix2y[1024];
35 47  
36 48  
37 49 /*==========================================================================
... ... @@ -125,7 +137,10 @@ GHealpix&amp; GHealpix::operator= (const GHealpix&amp; pixels)
125 137 /***********************************************************************//**
126 138 * @brief Load Healpix data from table.
127 139 *
128   - * @param hdu FITS HDU containing the Healpix data
  140 + * @param[in] hdu FITS HDU containing the Healpix data
  141 + *
  142 + * @exception GException::healpix
  143 + * Unable to load Healpix data from table.
129 144 ***************************************************************************/
130 145 void GHealpix::load(const GFitsHDU* hdu)
131 146 {
... ... @@ -133,56 +148,66 @@ void GHealpix::load(const GFitsHDU* hdu)
133 148 free_members();
134 149 init_members();
135 150  
136   - // Get HDU header
137   - GFitsHeader* hdr = hdu->header();
138   - if (hdr == NULL)
139   - throw;
140   -
141 151 // Check if we have a healpix representation
142   - if (hdr->card("PIXTYPE")->string() != "HEALPIX")
143   - throw;
  152 + if (hdu->card("PIXTYPE")->string() != "HEALPIX")
  153 + throw GException::healpix(G_LOAD, "HDU does not contain Healpix data");
144 154  
145   - // Get Healpix resolution and determine number of pixels
146   - m_nside = hdr->card("NSIDE")->integer();
147   - m_num_pixels = 12*m_nside*m_nside;
  155 + // Get Healpix resolution and determine number of pixels and solid angle
  156 + m_nside = hdu->card("NSIDE")->integer();
  157 + m_num_pixels = 12 * m_nside * m_nside;
  158 + m_omega = fourpi / m_num_pixels;
148 159  
149   - // Get ordering
  160 + // Get ordering scheme from ORDERING keyword
  161 + std::string ordering;
150 162 try {
151   - std::string ordering = hdr->card("ORDERING")->string();
152   - if (ordering == "RING")
153   - m_order = 0;
154   - else if (ordering == "NESTED")
155   - m_order = 1;
156   - else
157   - throw;
  163 + ordering = hdu->card("ORDERING")->string();
158 164 }
159   - catch (exception &e) {
160   - m_order = 1; // Default is NESTED
  165 + catch (GException::fits_key_not_found &e) {
  166 + m_order = -1; // Flag unknown ordering
161 167 }
162   -
163   - // Get coordinate system
  168 +
  169 + // Decode ordering scheme string
  170 + if (ordering == "RING")
  171 + m_order = 0;
  172 + else if (ordering == "NESTED")
  173 + m_order = 1;
  174 + else
  175 + throw GException::healpix(G_LOAD, "Invalid ordering");
  176 +
  177 + // Get coordinate system. First search for the COORDSYS keyword.
  178 + // If not found then search for the HIER_CRD keyword. This has
  179 + // been used in older versions of LAT exposure cubes ...
  180 + std::string coordsys;
164 181 try {
165   - std::string coordsys = hdr->card("COORDSYS")->string();
166   - if (coordsys == "EQUATORIAL")
167   - m_coordsys = 0;
168   - else if (coordsys == "GALACTIC")
169   - m_coordsys = 1;
170   - else
171   - throw;
  182 + coordsys = hdu->card("COORDSYS")->string();
172 183 }
173   - catch (exception &e) {
174   - m_coordsys = 0; // Default is EQUATORIAL
  184 + catch (GException::fits_key_not_found &e) {
  185 + try {
  186 + coordsys = hdu->card("HIER_CRD")->string();
  187 + }
  188 + catch (GException::fits_key_not_found &e) {
  189 + m_coordsys = -1; // Flag coordinate system as unknown
  190 + }
175 191 }
176 192  
  193 + // Decode coordinate system string
  194 + if (coordsys.find("EQU") == 0)
  195 + m_coordsys = 0;
  196 + else if (coordsys.find("GAL") == 0)
  197 + m_coordsys = 1;
  198 + else
  199 + throw GException::healpix(G_LOAD, "Invalid coordinate system");
  200 +
177 201 // Continue only of we have pixels
178 202 if (m_num_pixels > 0) {
179 203  
180 204 // Get first column
181   - GFitsTableCol* col = hdu->column(1);
  205 + GFitsTableCol* col = hdu->column(0);
182 206  
183 207 // Check column consistency
184 208 if (col->length() != m_num_pixels)
185   - throw;
  209 + throw GException::healpix(G_LOAD,
  210 + "NSIDE inconsistent with number of rows");
186 211  
187 212 // Extract vector size of each pixel
188 213 m_size_pixels = col->number();
... ... @@ -198,7 +223,10 @@ void GHealpix::load(const GFitsHDU* hdu)
198 223 }
199 224 } // endif: there were pixels to load
200 225  
201   - // Calculate longitude and latitude for each pixel
  226 + // Compute sky direction for each pixel
  227 + m_dir = new GSkyDir[m_num_pixels];
  228 + for (int i = 0; i < m_num_pixels; ++i)
  229 + m_dir[i] = pix2ang(i);
202 230  
203 231 } // endif: we had pixels
204 232  
... ... @@ -207,7 +235,6 @@ void GHealpix::load(const GFitsHDU* hdu)
207 235 }
208 236  
209 237  
210   -
211 238 /***********************************************************************//**
212 239 * @brief Returns number of divisions of the side of each base pixel.
213 240 ***************************************************************************/
... ... @@ -233,11 +260,49 @@ int GHealpix::num_pixels(void) const
233 260 ***************************************************************************/
234 261 double GHealpix::omega(void) const
235 262 {
236   - // Set solid angle
237   - double omega = fourpi / m_num_pixels;
238   -
239 263 // Return solid angle
240   - return omega;
  264 + return m_omega;
  265 +}
  266 +
  267 +
  268 +/***********************************************************************//**
  269 + * @brief Returns sky direction of pixel
  270 + *
  271 + * @param[in] ipix Pixel number (0,1,...,num_pixels)
  272 + ***************************************************************************/
  273 +GSkyDir GHealpix::pix2ang(const int& ipix)
  274 +{
  275 + // Declare result
  276 + GSkyDir result;
  277 + double theta = 0.0;
  278 + double phi = 0.0;
  279 +
  280 + // Perform ordering dependent conversion
  281 + switch (m_order) {
  282 + case 0:
  283 + pix2ang_ring(ipix, &theta, &phi);
  284 + break;
  285 + case 1:
  286 + pix2ang_nest(ipix, &theta, &phi);
  287 + break;
  288 + default:
  289 + break;
  290 + }
  291 +
  292 + // Store coordinate system dependent result
  293 + switch (m_coordsys) {
  294 + case 0:
  295 + result.radec(phi, pihalf-theta);
  296 + break;
  297 + case 1:
  298 + result.lb(phi, pihalf-theta);
  299 + break;
  300 + default:
  301 + break;
  302 + }
  303 +
  304 + // Return result
  305 + return result;
241 306 }
242 307  
243 308  
... ... @@ -258,9 +323,12 @@ void GHealpix::init_members(void)
258 323 m_coordsys = 0;
259 324 m_num_pixels = 0;
260 325 m_size_pixels = 0;
  326 + m_omega = 0.0;
261 327 m_pixels = NULL;
262   - m_lon = NULL;
263   - m_lat = NULL;
  328 + m_dir = NULL;
  329 +
  330 + // Initiates the array for the pixel number -> (x,y) mapping
  331 + mk_pix2xy();
264 332  
265 333 // Return
266 334 return;
... ... @@ -280,6 +348,7 @@ void GHealpix::copy_members(const GHealpix&amp; pixels)
280 348 m_coordsys = pixels.m_coordsys;
281 349 m_num_pixels = pixels.m_num_pixels;
282 350 m_size_pixels = pixels.m_size_pixels;
  351 + m_omega = pixels.m_omega;
283 352  
284 353 // Copy arrays
285 354 if (m_num_pixels > 0) {
... ... @@ -288,13 +357,10 @@ void GHealpix::copy_members(const GHealpix&amp; pixels)
288 357 m_pixels = new double[size];
289 358 memcpy(m_pixels, pixels.m_pixels, size*sizeof(double));
290 359 }
291   - if (pixels.m_lon != NULL) {
292   - m_lon = new double[m_num_pixels];
293   - memcpy(m_lon, pixels.m_lon, m_num_pixels*sizeof(double));
294   - }
295   - if (pixels.m_lat != NULL) {
296   - m_lat = new double[m_num_pixels];
297   - memcpy(m_lat, pixels.m_lat, m_num_pixels*sizeof(double));
  360 + if (pixels.m_dir != NULL) {
  361 + m_dir = new GSkyDir[m_num_pixels];
  362 + for (int i = 0; i < m_num_pixels; ++i)
  363 + m_dir[i] = pixels.m_dir[i];
298 364 }
299 365 }
300 366  
... ... @@ -310,13 +376,11 @@ void GHealpix::free_members(void)
310 376 {
311 377 // Free memory
312 378 if (m_pixels != NULL) delete [] m_pixels;
313   - if (m_lon != NULL) delete [] m_lon;
314   - if (m_lat != NULL) delete [] m_lat;
  379 + if (m_dir != NULL) delete [] m_dir;
315 380  
316 381 // Mark memory as free
317 382 m_pixels = NULL;
318   - m_lon = NULL;
319   - m_lat = NULL;
  383 + m_dir = NULL;
320 384  
321 385 // Return
322 386 return;
... ... @@ -332,12 +396,247 @@ GHealpix* GHealpix::clone(void) const
332 396 }
333 397  
334 398  
  399 +/***********************************************************************//**
  400 + * @brief Constructs pixel to (x,y) conversion array
  401 + *
  402 + * Constructs the array giving x and y in the face from pixel number for the
  403 + * nested (quad-cube like) ordering of pixels. The bits corresponding to x
  404 + * and y are interleaved in the pixel number one breaks up the pixel number
  405 + * by even and odd bits.
  406 + ***************************************************************************/
  407 +void GHealpix::mk_pix2xy(void)
  408 +{
  409 + // Loop over all pixels of array
  410 + for (int kpix = 0; kpix < 1024; ++kpix) {
  411 +
  412 + // Setup
  413 + int jpix = kpix;
  414 + int IX = 0;
  415 + int IY = 0;
  416 + int IP = 1; // Bit position (in x and y)
  417 +
  418 + // Go through all the bits
  419 + while (jpix != 0) {
  420 + int ID = (int)fmod(jpix,2); // Bit value (in kpix), goes in ix
  421 + jpix = jpix/2;
  422 + IX = ID*IP+IX;
  423 + ID = (int)fmod(jpix,2); // Bit value (in kpix), goes in iy
  424 + jpix = jpix/2;
  425 + IY = ID*IP+IY;
  426 + IP = 2*IP; // Next bit (in x and y)
  427 + }
  428 +
  429 + // Set array elements
  430 + pix2x[kpix] = IX; // in {0,31}
  431 + pix2y[kpix] = IY; // in {0,31}
  432 +
  433 + }
  434 +
  435 + // Return
  436 + return;
  437 +}
  438 +
  439 +
  440 +/***********************************************************************//**
  441 + * @brief Convert pixel index to (theta,phi) angles for ring ordering
  442 + *
  443 + * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
  444 + * @param[out] theta Pointer to result zenith angle in radians.
  445 + * @param[out] phi Pointer to result azimuth angle in radians.
  446 + *
  447 + * @exception GException::out_of_range Pixel index is out of range.
  448 + ***************************************************************************/
  449 +void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi)
  450 +{
  451 + // Check if ipix is in range
  452 + if (ipix < 0 || ipix >= m_num_pixels)
  453 + throw GException::out_of_range(G_PIX2ANG_RING, ipix, 0, m_num_pixels-1);
  454 +
  455 + // Setup
  456 + int ipix1 = ipix + 1; // in {1, m_num_pixels}
  457 + int nl2 = 2 * m_nside;
  458 + int nl4 = 4 * m_nside;
  459 +
  460 + // Determine points in each polar cap (0 for nside=1)
  461 + int ncap = 2 * m_nside * (m_nside-1);
  462 +
  463 + double fact1 = 1.5 * m_nside;
  464 + double fact2 = 3.0 * m_nside * m_nside;
  465 +
  466 + // Handle North Polar cap
  467 + if (ipix1 <= ncap) {
  468 + double hip = ipix1/2.;
  469 + double fihip = floor(hip);
  470 + int iring = (int)floor(sqrt(hip - sqrt(fihip)))+1;
  471 + int iphi = ipix1 - 2*iring*(iring - 1);
  472 + *theta = acos(1. - iring*iring / fact2);
  473 + *phi = (1.*iphi - 0.5) * pi/(2.*iring);
  474 + }
  475 +
  476 + // Handle Equatorial region
  477 + else if(ipix1 <= nl2*(5*m_nside+1)) {
  478 + int ip = ipix1 - ncap - 1;
  479 + int iring = (int)floor(ip/nl4) + m_nside;
  480 + int iphi = (int)fmod(ip,nl4) + 1;
  481 + double fodd = 0.5 * (1 + fmod((double)(iring+m_nside),2));
  482 + *theta = acos((nl2 - iring) / fact1);
  483 + *phi = (1.*iphi - fodd) * pi/(2.*m_nside);
  484 + }
  485 +
  486 + // Handle South Polar cap
  487 + else {
  488 + int ip = m_num_pixels - ipix1 + 1;
  489 + double hip = ip/2.;
  490 + double fihip = floor(hip);
  491 + int iring = (int)floor(sqrt( hip - sqrt(fihip))) + 1;
  492 + int iphi = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
  493 + *theta = acos( -1. + iring*iring / fact2 );
  494 + *phi = (1.*iphi - 0.5) * pi/(2.*iring);
  495 + }
  496 +
  497 + // Return
  498 + return;
  499 +}
  500 +
  501 +
  502 +/***********************************************************************//**
  503 + * @brief Convert pixel index to (theta,phi) angles for nested ordering
  504 + *
  505 + * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
  506 + * @param[out] theta Pointer to result zenith angle in radians.
  507 + * @param[out] phi Pointer to result azimuth angle in radians.
  508 + *
  509 + * @exception GException::out_of_range Pixel index is out of range.
  510 + ***************************************************************************/
  511 +void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi)
  512 +{
  513 + // Check if ipix is in range
  514 + if (ipix < 0 || ipix >= m_num_pixels)
  515 + throw GException::out_of_range(G_PIX2ANG_NEST, ipix, 0, m_num_pixels-1);
  516 +
  517 + // ...
  518 + double fn = 1. * m_nside;
  519 + double fact1 = 1. / (3.*fn*fn);
  520 + double fact2 = 2. / (3.*fn);
  521 + int nl4 = 4 * m_nside;
  522 +
  523 + // Finds the face, and the number in the face
  524 + int npface = m_nside*m_nside;
  525 + int face_num = ipix/npface; // Face number in {0,11}
  526 + int ipf = (int)fmod(ipix,npface); // Pixel number in the face {0,npface-1}
  527 +
  528 + // Finds the x,y on the face from the pixel number
  529 + // (starting from the lowest corner)
  530 + int ip_low = (int)fmod(ipf,1024); // Content of the last 10 bits
  531 + int ip_trunc = ipf/1024; // Truncation of the last 10 bits
  532 + int ip_med = (int)fmod(ip_trunc,1024); // Content of the next 10 bits
  533 + int ip_hi = ip_trunc/1024; // Content of the high weight 10 bits
  534 + int ix = 1024*pix2x[ip_hi] + 32*pix2x[ip_med] + pix2x[ip_low];
  535 + int iy = 1024*pix2y[ip_hi] + 32*pix2y[ip_med] + pix2y[ip_low];
  536 +
  537 + // Transforms this in (horizontal, vertical) coordinates
  538 + int jrt = ix + iy; // 'vertical' in {0,2*(m_nside-1)}
  539 + int jpt = ix - iy; // 'horizontal' in {-m_nside+1,m_nside-1}
  540 +
  541 + // Computes the z coordinate on the sphere
  542 + int jr = jrll[face_num]*m_nside - jrt - 1;
  543 + int nr = m_nside; // Equatorial region (the most frequent)
  544 + double z = (2*m_nside-jr)*fact2;
  545 + int kshift = (int)fmod(jr - m_nside, 2);
  546 +
  547 + // North pole region
  548 + if (jr < m_nside) {
  549 + nr = jr;
  550 + z = 1. - nr*nr*fact1;
  551 + kshift = 0;
  552 + }
  553 +
  554 + // South pole region
  555 + else if (jr > 3*m_nside) {
  556 + nr = nl4 - jr;
  557 + z = - 1. + nr*nr*fact1;
  558 + kshift = 0;
  559 + }
  560 +
  561 + // Computes theta
  562 + *theta = acos(z);
  563 +
  564 + // Computes the phi coordinate on the sphere, in [0,2Pi]
  565 + int jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
  566 + if (jp > nl4) jp = jp - nl4;
  567 + if (jp < 1) jp = jp + nl4;
  568 +
  569 + // Computes Phi
  570 + *phi = (jp - (kshift+1)*0.5) * (pihalf / nr);
  571 +
  572 + // Return
  573 + return;
  574 +}
  575 +
  576 +
335 577 /*==========================================================================
336 578 = =
337 579 = GHealpix friends =
338 580 = =
339 581 ==========================================================================*/
340 582  
  583 +/***********************************************************************//**
  584 + * @brief Output operator
  585 + *
  586 + * @param[in] os Output stream
  587 + * @param[in] column Healpix array to put in output stream
  588 + ***************************************************************************/
  589 +ostream& operator<< (ostream& os, const GHealpix& pixels)
  590 +{
  591 + // Put header in stream
  592 + os << "=== GHealpix ===" << endl;
  593 + os << " Number of base pixel div. .: " << pixels.m_nside << endl;
  594 + os << " Number of pixels ..........: " << pixels.m_num_pixels << endl;
  595 + os << " Pixel vector size .........: " << pixels.m_size_pixels << endl;
  596 + os << " Solid angle ...............: " << scientific << pixels.m_omega
  597 + << fixed << endl;
  598 +
  599 + // Put ordering in stream
  600 + os << " Ordering ..................: ";
  601 + switch (pixels.m_order) {
  602 + case 0:
  603 + os << "Ring" << endl;
  604 + break;
  605 + case 1:
  606 + os << "Nested" << endl;
  607 + break;
  608 + case -1:
  609 + os << "*** Unknown ***" << endl;
  610 + break;
  611 + default:
  612 + os << "*** Invalid ***" << endl;
  613 + break;
  614 + }
  615 +
  616 + // Put coordinate system in stream
  617 + os << " Coordinate system .........: ";
  618 + switch (pixels.m_coordsys) {
  619 + case 0:
  620 + os << "Equatorial (RA,Dec)" << endl;
  621 + break;
  622 + case 1:
  623 + os << "Galactic (l,b)" << endl;
  624 + break;
  625 + case -1:
  626 + os << "*** Unknown ***" << endl;
  627 + break;
  628 + default:
  629 + os << "*** Invalid ***" << endl;
  630 + break;
  631 + }
  632 +
  633 + for (int i = 0; i < pixels.m_num_pixels; ++i)
  634 + os << pixels.m_dir[i] << " ";
  635 +
  636 + // Return output stream
  637 + return os;
  638 +}
  639 +
341 640  
342 641 /*==========================================================================
343 642 = =
... ...
src/GSkyDir.cpp
  View file @fdfb897
... ... @@ -14,6 +14,7 @@
14 14  
15 15 /* __ Includes ___________________________________________________________ */
16 16 #include "GException.hpp"
  17 +#include "GTools.hpp"
17 18 #include "GSkyDir.hpp"
18 19  
19 20 /* __ Namespaces _________________________________________________________ */
... ... @@ -34,9 +35,8 @@
34 35 = =
35 36 ==========================================================================*/
36 37  
37   -/***************************************************************************
38   - * Constructor *
39   - * ----------------------------------------------------------------------- *
  38 +/***********************************************************************//**
  39 + * @brief Constructor
40 40 ***************************************************************************/
41 41 GSkyDir::GSkyDir()
42 42 {
... ... @@ -48,9 +48,10 @@ GSkyDir::GSkyDir()
48 48 }
49 49  
50 50  
51   -/***************************************************************************
52   - * Copy constructor *
53   - * ----------------------------------------------------------------------- *
  51 +/***********************************************************************//**
  52 + * @brief Copy constructor
  53 + *
  54 + * @param[in] dir Sky direction from which class should be instantiated.
54 55 ***************************************************************************/
55 56 GSkyDir::GSkyDir(const GSkyDir& dir)
56 57 {
... ... @@ -65,9 +66,8 @@ GSkyDir::GSkyDir(const GSkyDir&amp; dir)
65 66 }
66 67  
67 68  
68   -/***************************************************************************
69   - * Destructor *
70   - * ----------------------------------------------------------------------- *
  69 +/***********************************************************************//**
  70 + * @brief Destructor
71 71 ***************************************************************************/
72 72 GSkyDir::~GSkyDir()
73 73 {
... ... @@ -85,9 +85,8 @@ GSkyDir::~GSkyDir()
85 85 = =
86 86 ==========================================================================*/
87 87  
88   -/***************************************************************************
89   - * Assignment operator *
90   - * ----------------------------------------------------------------------- *
  88 +/***********************************************************************//**
  89 + * @brief Assignment operator
91 90 ***************************************************************************/
92 91 GSkyDir& GSkyDir::operator= (const GSkyDir& dir)
93 92 {
... ... @@ -116,6 +115,103 @@ GSkyDir&amp; GSkyDir::operator= (const GSkyDir&amp; dir)
116 115 = =
117 116 ==========================================================================*/
118 117  
  118 +/***********************************************************************//**
  119 + * @brief Set equatorial sky direction
  120 + *
  121 + * @param[in] ra Right Ascension in radians.
  122 + * @param[in] dec Declination in radians.
  123 + ***************************************************************************/
  124 +void GSkyDir::radec(const double& ra, const double& dec)
  125 +{
  126 + // Set attributes
  127 + m_has_lb = 0;
  128 + m_has_radec = 1;
  129 +
  130 + // Set direction
  131 + m_ra = ra;
  132 + m_dec = dec;
  133 +
  134 + // Return
  135 + return;
  136 +}
  137 +
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Set galactic sky direction
  141 + *
  142 + * @param[in] l Galactic longitude in radians.
  143 + * @param[in] b Galactic latitude in radians.
  144 + ***************************************************************************/
  145 +void GSkyDir::lb(const double& l, const double& b)
  146 +{
  147 + // Set attributes
  148 + m_has_lb = 1;
  149 + m_has_radec = 0;
  150 +
  151 + // Set direction
  152 + m_l = l;
  153 + m_b = b;
  154 +
  155 + // Return
  156 + return;
  157 +}
  158 +
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Returns galactic longitude in radians
  162 + ***************************************************************************/
  163 +double GSkyDir::l(void)
  164 +{
  165 + // If we have no galactic coordinates then get them now
  166 + if (!m_has_lb && m_has_radec)
  167 + equ2gal();
  168 +
  169 + // Return galactic longitude
  170 + return m_l;
  171 +}
  172 +
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Returns galactic latitude in radians
  176 + ***************************************************************************/
  177 +double GSkyDir::b(void)
  178 +{
  179 + // If we have no galactic coordinates then get them now
  180 + if (!m_has_lb && m_has_radec)
  181 + equ2gal();
  182 +
  183 + // Return galactic latitude
  184 + return m_b;
  185 +}
  186 +
  187 +
  188 +/***********************************************************************//**
  189 + * @brief Returns Right Ascension in radians
  190 + ***************************************************************************/
  191 +double GSkyDir::ra(void)
  192 +{
  193 + // If we have no equatorial coordinates then get them now
  194 + if (!m_has_radec && m_has_lb)
  195 + gal2equ();
  196 +
  197 + // Return Right Ascension
  198 + return m_ra;
  199 +}
  200 +
  201 +
  202 +/***********************************************************************//**
  203 + * @brief Returns Declination in radians
  204 + ***************************************************************************/
  205 +double GSkyDir::dec(void)
  206 +{
  207 + // If we have no equatorial coordinates then get them now
  208 + if (!m_has_radec && m_has_lb)
  209 + gal2equ();
  210 +
  211 + // Return Declination
  212 + return m_dec;
  213 +}
  214 +
119 215  
120 216 /*==========================================================================
121 217 = =
... ... @@ -123,37 +219,46 @@ GSkyDir&amp; GSkyDir::operator= (const GSkyDir&amp; dir)
123 219 = =
124 220 ==========================================================================*/
125 221  
126   -/***************************************************************************
127   - * Initialise class members *
128   - * ----------------------------------------------------------------------- *
  222 +/***********************************************************************//**
  223 + * @brief Initialise class members
129 224 ***************************************************************************/
130 225 void GSkyDir::init_members(void)
131 226 {
132 227 // Initialise members
  228 + m_has_lb = 0;
  229 + m_has_radec = 0;
  230 + m_l = 0.0;
  231 + m_b = 0.0;
  232 + m_ra = 0.0;
  233 + m_dec = 0.0;
133 234  
134 235 // Return
135 236 return;
136 237 }
137 238  
138 239  
139   -/***************************************************************************
140   - * Copy class members *
141   - * ----------------------------------------------------------------------- *
  240 +/***********************************************************************//**
  241 + * @brief Copy class members
  242 + *
  243 + * @param[in] dir Sky direction from which members should be copied
142 244 ***************************************************************************/
143 245 void GSkyDir::copy_members(const GSkyDir& dir)
144 246 {
145 247 // Copy attributes
146   -
147   - // Copy other membres
  248 + m_has_lb = dir.m_has_lb;
  249 + m_has_radec = dir.m_has_radec;
  250 + m_l = dir.m_l;
  251 + m_b = dir.m_b;
  252 + m_ra = dir.m_ra;
  253 + m_dec = dir.m_dec;
148 254  
149 255 // Return
150 256 return;
151 257 }
152 258  
153 259  
154   -/***************************************************************************
155   - * Delete class members *
156   - * ----------------------------------------------------------------------- *
  260 +/***********************************************************************//**
  261 + * @brief Delete class members
157 262 ***************************************************************************/
158 263 void GSkyDir::free_members(void)
159 264 {
... ... @@ -166,12 +271,54 @@ void GSkyDir::free_members(void)
166 271 }
167 272  
168 273  
  274 +/***********************************************************************//**
  275 + * @brief Convert equatorial to galactic coordinates
  276 + *
  277 + * @todo Needs to be implemented
  278 + ***************************************************************************/
  279 +void GSkyDir::equ2gal(void)
  280 +{
  281 + // Return
  282 + return;
  283 +}
  284 +
  285 +
  286 +/***********************************************************************//**
  287 + * @brief Convert galactic to equatorial coordinates
  288 + *
  289 + * @todo Needs to be implemented
  290 + ***************************************************************************/
  291 +void GSkyDir::gal2equ(void)
  292 +{
  293 + // Return
  294 + return;
  295 +}
  296 +
  297 +
169 298 /*==========================================================================
170 299 = =
171 300 = GSkyDir friends =
172 301 = =
173 302 ==========================================================================*/
174 303  
  304 +/***********************************************************************//**
  305 + * @brief Output operator
  306 + *
  307 + * @param[in] os Output stream
  308 + * @param[in] column Sky direction to put in output stream
  309 + ***************************************************************************/
  310 +ostream& operator<< (ostream& os, const GSkyDir& dir)
  311 +{
  312 + // Create coordinate system dependent output
  313 + if (dir.m_has_lb)
  314 + os << "(l,b)=(" << dir.m_l*rad2deg << "," << dir.m_b*rad2deg << ")";
  315 + else if (dir.m_has_radec)
  316 + os << "(RA,Dec)=(" << dir.m_ra*rad2deg << "," << dir.m_dec*rad2deg << ")";
  317 +
  318 + // Return output stream
  319 + return os;
  320 +}
  321 +
175 322  
176 323 /*==========================================================================
177 324 = =
... ...
test/test_python.py
  View file @fdfb897
... ... @@ -9,7 +9,7 @@ import os
9 9 #================#
10 10 def test_fits():
11 11 """
12   - Test GammaLib GFits interface.
  12 + Test GFits interface.
13 13 """
14 14 # Remove test files
15 15 try:
... ... @@ -120,7 +120,7 @@ def test_fits():
120 120 #=================#
121 121 def test_node_array():
122 122 """
123   - Test GammaLib GNodeArray interface.
  123 + Test GNodeArray interface.
124 124 """
125 125 # Set-up vector and data array
126 126 vector = GVector(20)
... ... @@ -153,7 +153,7 @@ def test_node_array():
153 153 #===================#
154 154 def test_lat_response():
155 155 """
156   - Test GammaLib GLATResponse interface.
  156 + Test GLATResponse interface.
157 157 """
158 158 # Remove test file
159 159 try:
... ... @@ -179,14 +179,40 @@ def test_lat_response():
179 179 #======================#
180 180 def test_lat_observation():
181 181 """
182   - Test GammaLib GLATObservation interface.
  182 + Test GLATObservation interface.
183 183 """
184 184 # Allocate LAT Observation
185 185 obs = GLATObservation("data/FT1_253582800.fits.gz", "data/FT2_253582800.fits.gz")
186 186  
187 187 #
188   - print obs.ft1()
189   - print obs.ft2()
  188 + #print obs.ft1()
  189 + #print obs.ft2()
  190 +
  191 +
  192 +#==============#
  193 +# Test Healpix #
  194 +#==============#
  195 +def test_healpix():
  196 + """
  197 + Test GHealpix interface.
  198 + """
  199 + # Open exposure cube
  200 + fits = GFits()
  201 + fits.open("data/expCube_253582800.fits.gz")
  202 +
  203 + # Lood Healpix table
  204 + pixels = GHealpix()
  205 + pixels.load(fits.hdu("EXPOSURE"))
  206 +
  207 + # Dump information
  208 + print pixels
  209 +
  210 + print pixels.pix2ang(0)
  211 + print " ",pixels.pix2ang(0).ra()*180/pi, pixels.pix2ang(0).dec()*180/pi
  212 + print pixels.pix2ang(1)
  213 + print pixels.pix2ang(2)
  214 + print pixels.pix2ang(3)
  215 +
190 216  
191 217 #==========================#
192 218 # Main routine entry point #
... ... @@ -199,4 +225,5 @@ if __name__ == &#39;__main__&#39;:
199 225 test_node_array()
200 226 test_lat_response()
201 227 test_lat_observation()
  228 + test_healpix()
202 229  
203 230 \ No newline at end of file
... ...