From fdfb897d4246eadfef1a1ef56d5b4e0fef4daaa8 Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Fri, 14 Mar 2008 07:41:54 +0000 Subject: [PATCH] enhance GHealpix and GSkyDir classes --- include/GException.hpp | 6 ++++++ include/GFitsHDU.hpp | 2 +- include/GHealpix.hpp | 34 +++++++++++++++++++++------------- include/GSkyDir.hpp | 60 +++++++++++++++++++++++++++--------------------------------- pyext/GHealpix.i | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/GLATResponse.i | 2 +- pyext/GSkyDir.i | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/gammalib.i | 2 ++ src/GException.cpp | 13 +++++++++++++ src/GHealpix.cpp | 405 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------- src/GSkyDir.cpp | 193 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------- test/test_python.py | 39 +++++++++++++++++++++++++++++++++------ 12 files changed, 749 insertions(+), 130 deletions(-) create mode 100644 pyext/GHealpix.i create mode 100644 pyext/GSkyDir.i diff --git a/include/GException.hpp b/include/GException.hpp index 3744d6d..a3fc74b 100644 --- a/include/GException.hpp +++ b/include/GException.hpp @@ -233,6 +233,12 @@ public: rsp_invalid_type(string origin, string type); }; + // General Healpix error + class healpix : public GExceptionHandler { + public: + healpix(string origin, string message); + }; + }; #endif /* GEXCEPTION_HPP */ diff --git a/include/GFitsHDU.hpp b/include/GFitsHDU.hpp index 2820477..c0e2f99 100644 --- a/include/GFitsHDU.hpp +++ b/include/GFitsHDU.hpp @@ -85,7 +85,7 @@ private: // Private data area __fitsfile m_fitsfile; //!< FITS file pointer - int m_hdunum; //!< HDU number (starting from 1) + int m_hdunum; //!< HDU number (starting from 0) std::string m_name; //!< HDU name (extname) int m_type; //!< HDU type GFitsHeader* m_header; //!< HDU header diff --git a/include/GHealpix.hpp b/include/GHealpix.hpp index 3660739..f7f32dd 100644 --- a/include/GHealpix.hpp +++ b/include/GHealpix.hpp @@ -21,6 +21,7 @@ /* __ Includes ___________________________________________________________ */ #include "GFitsHDU.hpp" +#include "GSkyDir.hpp" /* __ Namespaces _________________________________________________________ */ @@ -32,6 +33,9 @@ ***************************************************************************/ class GHealpix { + // I/O friends + friend ostream& operator<< (ostream& os, const GHealpix& pixels); + public: // Constructors and destructors GHealpix(); @@ -42,10 +46,11 @@ public: GHealpix& operator= (const GHealpix& pixels); // Methods - void load(const GFitsHDU* hdu); - int nside(void) const; - int num_pixels(void) const; - double omega(void) const; + void load(const GFitsHDU* hdu); + int nside(void) const; + int num_pixels(void) const; + double omega(void) const; + GSkyDir pix2ang(const int& ipix); private: // Private methods @@ -53,16 +58,19 @@ private: void copy_members(const GHealpix& pixels); void free_members(void); GHealpix* clone(void) const; - + void mk_pix2xy(void); + void pix2ang_ring(int ipix, double* theta, double* phi); + void pix2ang_nest(int ipix, double* theta, double* phi); + // Private data area - int m_nside; //!< Number of divisions of each base pixel (1,2,..) - int m_order; //!< Ordering (0=ring, 1=nested) - int m_coordsys; //!< Coordinate system (0=equatorial, 1=galactic) - int m_num_pixels; //!< Number of pixels - int m_size_pixels; //!< Vector size of each pixel - double* m_pixels; //!< Pixel array - double* m_lon; //!< RA or GLON values for each pixel - double* m_lat; //!< DEC or GLAT values for each pixel + int m_nside; //!< Number of divisions of each base pixel (1-8192) + int m_order; //!< Ordering (0=ring, 1=nested, -1=?) + int m_coordsys; //!< Coordinate system (0=equatorial, 1=galactic, -1=?) + int m_num_pixels; //!< Number of pixels + int m_size_pixels; //!< Vector size of each pixel + double m_omega; //!< Solid angle of pixel + double* m_pixels; //!< Pixel array + GSkyDir* m_dir; //!< Sky direction for each pixel }; #endif /* GHEALPIX_HPP */ diff --git a/include/GSkyDir.hpp b/include/GSkyDir.hpp index 4581e04..3275d4c 100644 --- a/include/GSkyDir.hpp +++ b/include/GSkyDir.hpp @@ -24,54 +24,48 @@ /* __ Namespaces _________________________________________________________ */ -/*************************************************************************** - * GSkyDir class definition * - ***************************************************************************/ -/** +/***********************************************************************//** * @class GSkyDir * - * @brief Interface for the sky direction classes. - * - * @author J. Knodlseder - */ + * @brief GSkyDir class interface defintion + ***************************************************************************/ class GSkyDir { + // I/O friends + friend ostream& operator<< (ostream& os, const GSkyDir& dir); + public: - /// Constructor + // Constructors and destructors GSkyDir(); - - /// Copy constructor GSkyDir(const GSkyDir& dir); - - /// Destructor virtual ~GSkyDir(); - /// Assignment operator + // Operators GSkyDir& operator= (const GSkyDir& dir); -protected: - /// Longitude - double m_lon; - - /// Latitude - double m_lat; + // Methods + void radec(const double& ra, const double& dec); + void lb(const double& l, const double& b); + double l(void); + double b(void); + double ra(void); + double dec(void); - /// Protected method that initialises protected members of the class +private: + // Private methods void init_members(void); - - /// Protected method that copies protected members of the class void copy_members(const GSkyDir& dir); - - /// Protected method that frees memory allocated by instances of the class void free_members(void); - -private: - + void equ2gal(void); + void gal2equ(void); + + // Private data area + int m_has_lb; //!< Has galactic coordinates + int m_has_radec; //!< Has equatorial coordinates + double m_l; //!< Galactic longitude in radians + double m_b; //!< Galactic latitude in radians + double m_ra; //!< Right Ascension in radians + double m_dec; //!< Declination in radians }; - -/*************************************************************************** - * Inline methods * - ***************************************************************************/ - #endif /* GSKYDIR_HPP */ diff --git a/pyext/GHealpix.i b/pyext/GHealpix.i new file mode 100644 index 0000000..59682d6 --- /dev/null +++ b/pyext/GHealpix.i @@ -0,0 +1,61 @@ +/*************************************************************************** + * GHealpix.i - Healpix sky representation class SWIG definition * + * ----------------------------------------------------------------------- * + * copyright : (C) 2008 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file GHealpix.i + * @brief GHealpix class SWIG file. + * @author J. Knodlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GHealpix.hpp" +%} + + +/***********************************************************************//** + * @class GHealpix + * + * @brief GHealpix class interface defintion + ***************************************************************************/ +class GHealpix { +public: + // Constructors and destructors + GHealpix(); + GHealpix(const GHealpix& pixels); + virtual ~GHealpix(); + + // Methods + void load(const GFitsHDU* hdu); + int nside(void) const; + int num_pixels(void) const; + double omega(void) const; + GSkyDir pix2ang(const int& ipix); +}; + + +/***********************************************************************//** + * @brief GHealpix class extension + ***************************************************************************/ +%extend GHealpix { + char *__str__() { + static char str_buffer[1001]; + std::ostringstream buffer; + buffer << *self; + std::string str = buffer.str(); + strncpy(str_buffer, (char*)str.c_str(), 1001); + str_buffer[1000] = '\0'; + return str_buffer; + } + GHealpix copy() { + return (*self); + } +}; diff --git a/pyext/GLATResponse.i b/pyext/GLATResponse.i index a5fd7d8..5326539 100644 --- a/pyext/GLATResponse.i +++ b/pyext/GLATResponse.i @@ -1,5 +1,5 @@ /*************************************************************************** - * GLATResponse.hpp - GLAST LAT Response class SWIG definition * + * GLATResponse.i - GLAST LAT Response class SWIG definition * * ----------------------------------------------------------------------- * * copyright : (C) 2008 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * diff --git a/pyext/GSkyDir.i b/pyext/GSkyDir.i new file mode 100644 index 0000000..8ea2dda --- /dev/null +++ b/pyext/GSkyDir.i @@ -0,0 +1,62 @@ +/*************************************************************************** + * GSkyDir.i - Sky direction class SWIG file * + * ----------------------------------------------------------------------- * + * copyright : (C) 2008 by Jurgen Knodlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +/** + * @file GSkyDir.i + * @brief GSkyDir class SWIG file. + * @author J. Knodlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GSkyDir.hpp" +%} + + +/***********************************************************************//** + * @class GSkyDir + * + * @brief GSkyDir class interface defintion + ***************************************************************************/ +class GSkyDir { +public: + // Constructors and destructors + GSkyDir(); + GSkyDir(const GSkyDir& dir); + virtual ~GSkyDir(); + + // Methods + void radec(const double& ra, const double& dec); + void lb(const double& l, const double& b); + double l(void); + double b(void); + double ra(void); + double dec(void); +}; + + +/***********************************************************************//** + * @brief GSkyDir class extension + ***************************************************************************/ +%extend GSkyDir { + char *__str__() { + static char str_buffer[1001]; + std::ostringstream buffer; + buffer << *self; + std::string str = buffer.str(); + strncpy(str_buffer, (char*)str.c_str(), 1001); + str_buffer[1000] = '\0'; + return str_buffer; + } + GSkyDir copy() { + return (*self); + } +}; diff --git a/pyext/gammalib.i b/pyext/gammalib.i index 92c4023..869291f 100755 --- a/pyext/gammalib.i +++ b/pyext/gammalib.i @@ -18,6 +18,8 @@ /* General */ #include "GException.i" #include "GNodeArray.i" +#include "GHealpix.i" +#include "GSkyDir.i" /* Response classes */ #include "GResponse.i" diff --git a/src/GException.cpp b/src/GException.cpp index 2032115..d6b037b 100644 --- a/src/GException.cpp +++ b/src/GException.cpp @@ -519,3 +519,16 @@ GException::rsp_invalid_type::rsp_invalid_type(string origin, m_origin = origin; m_message = "Invalid response type '"+type+"' specified"; } + + +/***********************************************************************//** + * @brief General Healpix error + * + * @param[in] origin Method that throws the error. + * @param[in] message Error message. + ***************************************************************************/ +GException::healpix::healpix(string origin, string message) +{ + m_origin = origin; + m_message = message; +} diff --git a/src/GHealpix.cpp b/src/GHealpix.cpp index 19161c3..823d69e 100644 --- a/src/GHealpix.cpp +++ b/src/GHealpix.cpp @@ -18,6 +18,8 @@ */ /* __ Includes ___________________________________________________________ */ +#include +#include #include "GException.hpp" #include "GTools.hpp" #include "GHealpix.hpp" @@ -25,6 +27,9 @@ /* __ Namespaces _________________________________________________________ */ /* __ Method name definitions ____________________________________________ */ +#define G_LOAD "GHealpix::load(const GFitsHDU*)" +#define G_PIX2ANG_RING "GHealpix::pix2ang_ring(int,double*,double*)" +#define G_PIX2ANG_NEST "GHealpix::pix2ang_nest(int,double*,double*)" /* __ Macros _____________________________________________________________ */ @@ -32,6 +37,13 @@ /* __ Debug definitions __________________________________________________ */ +/* __ Constants __________________________________________________________ */ +const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4}; +const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7}; + +/* __ Static conversion arrays ___________________________________________ */ +static int pix2x[1024]; +static int pix2y[1024]; /*========================================================================== @@ -125,7 +137,10 @@ GHealpix& GHealpix::operator= (const GHealpix& pixels) /***********************************************************************//** * @brief Load Healpix data from table. * - * @param hdu FITS HDU containing the Healpix data + * @param[in] hdu FITS HDU containing the Healpix data + * + * @exception GException::healpix + * Unable to load Healpix data from table. ***************************************************************************/ void GHealpix::load(const GFitsHDU* hdu) { @@ -133,56 +148,66 @@ void GHealpix::load(const GFitsHDU* hdu) free_members(); init_members(); - // Get HDU header - GFitsHeader* hdr = hdu->header(); - if (hdr == NULL) - throw; - // Check if we have a healpix representation - if (hdr->card("PIXTYPE")->string() != "HEALPIX") - throw; + if (hdu->card("PIXTYPE")->string() != "HEALPIX") + throw GException::healpix(G_LOAD, "HDU does not contain Healpix data"); - // Get Healpix resolution and determine number of pixels - m_nside = hdr->card("NSIDE")->integer(); - m_num_pixels = 12*m_nside*m_nside; + // Get Healpix resolution and determine number of pixels and solid angle + m_nside = hdu->card("NSIDE")->integer(); + m_num_pixels = 12 * m_nside * m_nside; + m_omega = fourpi / m_num_pixels; - // Get ordering + // Get ordering scheme from ORDERING keyword + std::string ordering; try { - std::string ordering = hdr->card("ORDERING")->string(); - if (ordering == "RING") - m_order = 0; - else if (ordering == "NESTED") - m_order = 1; - else - throw; + ordering = hdu->card("ORDERING")->string(); } - catch (exception &e) { - m_order = 1; // Default is NESTED + catch (GException::fits_key_not_found &e) { + m_order = -1; // Flag unknown ordering } - - // Get coordinate system + + // Decode ordering scheme string + if (ordering == "RING") + m_order = 0; + else if (ordering == "NESTED") + m_order = 1; + else + throw GException::healpix(G_LOAD, "Invalid ordering"); + + // Get coordinate system. First search for the COORDSYS keyword. + // If not found then search for the HIER_CRD keyword. This has + // been used in older versions of LAT exposure cubes ... + std::string coordsys; try { - std::string coordsys = hdr->card("COORDSYS")->string(); - if (coordsys == "EQUATORIAL") - m_coordsys = 0; - else if (coordsys == "GALACTIC") - m_coordsys = 1; - else - throw; + coordsys = hdu->card("COORDSYS")->string(); } - catch (exception &e) { - m_coordsys = 0; // Default is EQUATORIAL + catch (GException::fits_key_not_found &e) { + try { + coordsys = hdu->card("HIER_CRD")->string(); + } + catch (GException::fits_key_not_found &e) { + m_coordsys = -1; // Flag coordinate system as unknown + } } + // Decode coordinate system string + if (coordsys.find("EQU") == 0) + m_coordsys = 0; + else if (coordsys.find("GAL") == 0) + m_coordsys = 1; + else + throw GException::healpix(G_LOAD, "Invalid coordinate system"); + // Continue only of we have pixels if (m_num_pixels > 0) { // Get first column - GFitsTableCol* col = hdu->column(1); + GFitsTableCol* col = hdu->column(0); // Check column consistency if (col->length() != m_num_pixels) - throw; + throw GException::healpix(G_LOAD, + "NSIDE inconsistent with number of rows"); // Extract vector size of each pixel m_size_pixels = col->number(); @@ -198,7 +223,10 @@ void GHealpix::load(const GFitsHDU* hdu) } } // endif: there were pixels to load - // Calculate longitude and latitude for each pixel + // Compute sky direction for each pixel + m_dir = new GSkyDir[m_num_pixels]; + for (int i = 0; i < m_num_pixels; ++i) + m_dir[i] = pix2ang(i); } // endif: we had pixels @@ -207,7 +235,6 @@ void GHealpix::load(const GFitsHDU* hdu) } - /***********************************************************************//** * @brief Returns number of divisions of the side of each base pixel. ***************************************************************************/ @@ -233,11 +260,49 @@ int GHealpix::num_pixels(void) const ***************************************************************************/ double GHealpix::omega(void) const { - // Set solid angle - double omega = fourpi / m_num_pixels; - // Return solid angle - return omega; + return m_omega; +} + + +/***********************************************************************//** + * @brief Returns sky direction of pixel + * + * @param[in] ipix Pixel number (0,1,...,num_pixels) + ***************************************************************************/ +GSkyDir GHealpix::pix2ang(const int& ipix) +{ + // Declare result + GSkyDir result; + double theta = 0.0; + double phi = 0.0; + + // Perform ordering dependent conversion + switch (m_order) { + case 0: + pix2ang_ring(ipix, &theta, &phi); + break; + case 1: + pix2ang_nest(ipix, &theta, &phi); + break; + default: + break; + } + + // Store coordinate system dependent result + switch (m_coordsys) { + case 0: + result.radec(phi, pihalf-theta); + break; + case 1: + result.lb(phi, pihalf-theta); + break; + default: + break; + } + + // Return result + return result; } @@ -258,9 +323,12 @@ void GHealpix::init_members(void) m_coordsys = 0; m_num_pixels = 0; m_size_pixels = 0; + m_omega = 0.0; m_pixels = NULL; - m_lon = NULL; - m_lat = NULL; + m_dir = NULL; + + // Initiates the array for the pixel number -> (x,y) mapping + mk_pix2xy(); // Return return; @@ -280,6 +348,7 @@ void GHealpix::copy_members(const GHealpix& pixels) m_coordsys = pixels.m_coordsys; m_num_pixels = pixels.m_num_pixels; m_size_pixels = pixels.m_size_pixels; + m_omega = pixels.m_omega; // Copy arrays if (m_num_pixels > 0) { @@ -288,13 +357,10 @@ void GHealpix::copy_members(const GHealpix& pixels) m_pixels = new double[size]; memcpy(m_pixels, pixels.m_pixels, size*sizeof(double)); } - if (pixels.m_lon != NULL) { - m_lon = new double[m_num_pixels]; - memcpy(m_lon, pixels.m_lon, m_num_pixels*sizeof(double)); - } - if (pixels.m_lat != NULL) { - m_lat = new double[m_num_pixels]; - memcpy(m_lat, pixels.m_lat, m_num_pixels*sizeof(double)); + if (pixels.m_dir != NULL) { + m_dir = new GSkyDir[m_num_pixels]; + for (int i = 0; i < m_num_pixels; ++i) + m_dir[i] = pixels.m_dir[i]; } } @@ -310,13 +376,11 @@ void GHealpix::free_members(void) { // Free memory if (m_pixels != NULL) delete [] m_pixels; - if (m_lon != NULL) delete [] m_lon; - if (m_lat != NULL) delete [] m_lat; + if (m_dir != NULL) delete [] m_dir; // Mark memory as free m_pixels = NULL; - m_lon = NULL; - m_lat = NULL; + m_dir = NULL; // Return return; @@ -332,12 +396,247 @@ GHealpix* GHealpix::clone(void) const } +/***********************************************************************//** + * @brief Constructs pixel to (x,y) conversion array + * + * Constructs the array giving x and y in the face from pixel number for the + * nested (quad-cube like) ordering of pixels. The bits corresponding to x + * and y are interleaved in the pixel number one breaks up the pixel number + * by even and odd bits. + ***************************************************************************/ +void GHealpix::mk_pix2xy(void) +{ + // Loop over all pixels of array + for (int kpix = 0; kpix < 1024; ++kpix) { + + // Setup + int jpix = kpix; + int IX = 0; + int IY = 0; + int IP = 1; // Bit position (in x and y) + + // Go through all the bits + while (jpix != 0) { + int ID = (int)fmod(jpix,2); // Bit value (in kpix), goes in ix + jpix = jpix/2; + IX = ID*IP+IX; + ID = (int)fmod(jpix,2); // Bit value (in kpix), goes in iy + jpix = jpix/2; + IY = ID*IP+IY; + IP = 2*IP; // Next bit (in x and y) + } + + // Set array elements + pix2x[kpix] = IX; // in {0,31} + pix2y[kpix] = IY; // in {0,31} + + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Convert pixel index to (theta,phi) angles for ring ordering + * + * @param[in] ipix Pixel index for which (theta,phi) are to be computed. + * @param[out] theta Pointer to result zenith angle in radians. + * @param[out] phi Pointer to result azimuth angle in radians. + * + * @exception GException::out_of_range Pixel index is out of range. + ***************************************************************************/ +void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi) +{ + // Check if ipix is in range + if (ipix < 0 || ipix >= m_num_pixels) + throw GException::out_of_range(G_PIX2ANG_RING, ipix, 0, m_num_pixels-1); + + // Setup + int ipix1 = ipix + 1; // in {1, m_num_pixels} + int nl2 = 2 * m_nside; + int nl4 = 4 * m_nside; + + // Determine points in each polar cap (0 for nside=1) + int ncap = 2 * m_nside * (m_nside-1); + + double fact1 = 1.5 * m_nside; + double fact2 = 3.0 * m_nside * m_nside; + + // Handle North Polar cap + if (ipix1 <= ncap) { + double hip = ipix1/2.; + double fihip = floor(hip); + int iring = (int)floor(sqrt(hip - sqrt(fihip)))+1; + int iphi = ipix1 - 2*iring*(iring - 1); + *theta = acos(1. - iring*iring / fact2); + *phi = (1.*iphi - 0.5) * pi/(2.*iring); + } + + // Handle Equatorial region + else if(ipix1 <= nl2*(5*m_nside+1)) { + int ip = ipix1 - ncap - 1; + int iring = (int)floor(ip/nl4) + m_nside; + int iphi = (int)fmod(ip,nl4) + 1; + double fodd = 0.5 * (1 + fmod((double)(iring+m_nside),2)); + *theta = acos((nl2 - iring) / fact1); + *phi = (1.*iphi - fodd) * pi/(2.*m_nside); + } + + // Handle South Polar cap + else { + int ip = m_num_pixels - ipix1 + 1; + double hip = ip/2.; + double fihip = floor(hip); + int iring = (int)floor(sqrt( hip - sqrt(fihip))) + 1; + int iphi = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1))); + *theta = acos( -1. + iring*iring / fact2 ); + *phi = (1.*iphi - 0.5) * pi/(2.*iring); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Convert pixel index to (theta,phi) angles for nested ordering + * + * @param[in] ipix Pixel index for which (theta,phi) are to be computed. + * @param[out] theta Pointer to result zenith angle in radians. + * @param[out] phi Pointer to result azimuth angle in radians. + * + * @exception GException::out_of_range Pixel index is out of range. + ***************************************************************************/ +void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) +{ + // Check if ipix is in range + if (ipix < 0 || ipix >= m_num_pixels) + throw GException::out_of_range(G_PIX2ANG_NEST, ipix, 0, m_num_pixels-1); + + // ... + double fn = 1. * m_nside; + double fact1 = 1. / (3.*fn*fn); + double fact2 = 2. / (3.*fn); + int nl4 = 4 * m_nside; + + // Finds the face, and the number in the face + int npface = m_nside*m_nside; + int face_num = ipix/npface; // Face number in {0,11} + int ipf = (int)fmod(ipix,npface); // Pixel number in the face {0,npface-1} + + // Finds the x,y on the face from the pixel number + // (starting from the lowest corner) + int ip_low = (int)fmod(ipf,1024); // Content of the last 10 bits + int ip_trunc = ipf/1024; // Truncation of the last 10 bits + int ip_med = (int)fmod(ip_trunc,1024); // Content of the next 10 bits + int ip_hi = ip_trunc/1024; // Content of the high weight 10 bits + int ix = 1024*pix2x[ip_hi] + 32*pix2x[ip_med] + pix2x[ip_low]; + int iy = 1024*pix2y[ip_hi] + 32*pix2y[ip_med] + pix2y[ip_low]; + + // Transforms this in (horizontal, vertical) coordinates + int jrt = ix + iy; // 'vertical' in {0,2*(m_nside-1)} + int jpt = ix - iy; // 'horizontal' in {-m_nside+1,m_nside-1} + + // Computes the z coordinate on the sphere + int jr = jrll[face_num]*m_nside - jrt - 1; + int nr = m_nside; // Equatorial region (the most frequent) + double z = (2*m_nside-jr)*fact2; + int kshift = (int)fmod(jr - m_nside, 2); + + // North pole region + if (jr < m_nside) { + nr = jr; + z = 1. - nr*nr*fact1; + kshift = 0; + } + + // South pole region + else if (jr > 3*m_nside) { + nr = nl4 - jr; + z = - 1. + nr*nr*fact1; + kshift = 0; + } + + // Computes theta + *theta = acos(z); + + // Computes the phi coordinate on the sphere, in [0,2Pi] + int jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2; + if (jp > nl4) jp = jp - nl4; + if (jp < 1) jp = jp + nl4; + + // Computes Phi + *phi = (jp - (kshift+1)*0.5) * (pihalf / nr); + + // Return + return; +} + + /*========================================================================== = = = GHealpix friends = = = ==========================================================================*/ +/***********************************************************************//** + * @brief Output operator + * + * @param[in] os Output stream + * @param[in] column Healpix array to put in output stream + ***************************************************************************/ +ostream& operator<< (ostream& os, const GHealpix& pixels) +{ + // Put header in stream + os << "=== GHealpix ===" << endl; + os << " Number of base pixel div. .: " << pixels.m_nside << endl; + os << " Number of pixels ..........: " << pixels.m_num_pixels << endl; + os << " Pixel vector size .........: " << pixels.m_size_pixels << endl; + os << " Solid angle ...............: " << scientific << pixels.m_omega + << fixed << endl; + + // Put ordering in stream + os << " Ordering ..................: "; + switch (pixels.m_order) { + case 0: + os << "Ring" << endl; + break; + case 1: + os << "Nested" << endl; + break; + case -1: + os << "*** Unknown ***" << endl; + break; + default: + os << "*** Invalid ***" << endl; + break; + } + + // Put coordinate system in stream + os << " Coordinate system .........: "; + switch (pixels.m_coordsys) { + case 0: + os << "Equatorial (RA,Dec)" << endl; + break; + case 1: + os << "Galactic (l,b)" << endl; + break; + case -1: + os << "*** Unknown ***" << endl; + break; + default: + os << "*** Invalid ***" << endl; + break; + } + + for (int i = 0; i < pixels.m_num_pixels; ++i) + os << pixels.m_dir[i] << " "; + + // Return output stream + return os; +} + /*========================================================================== = = diff --git a/src/GSkyDir.cpp b/src/GSkyDir.cpp index b0d0707..13ce874 100644 --- a/src/GSkyDir.cpp +++ b/src/GSkyDir.cpp @@ -14,6 +14,7 @@ /* __ Includes ___________________________________________________________ */ #include "GException.hpp" +#include "GTools.hpp" #include "GSkyDir.hpp" /* __ Namespaces _________________________________________________________ */ @@ -34,9 +35,8 @@ = = ==========================================================================*/ -/*************************************************************************** - * Constructor * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Constructor ***************************************************************************/ GSkyDir::GSkyDir() { @@ -48,9 +48,10 @@ GSkyDir::GSkyDir() } -/*************************************************************************** - * Copy constructor * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] dir Sky direction from which class should be instantiated. ***************************************************************************/ GSkyDir::GSkyDir(const GSkyDir& dir) { @@ -65,9 +66,8 @@ GSkyDir::GSkyDir(const GSkyDir& dir) } -/*************************************************************************** - * Destructor * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Destructor ***************************************************************************/ GSkyDir::~GSkyDir() { @@ -85,9 +85,8 @@ GSkyDir::~GSkyDir() = = ==========================================================================*/ -/*************************************************************************** - * Assignment operator * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Assignment operator ***************************************************************************/ GSkyDir& GSkyDir::operator= (const GSkyDir& dir) { @@ -116,6 +115,103 @@ GSkyDir& GSkyDir::operator= (const GSkyDir& dir) = = ==========================================================================*/ +/***********************************************************************//** + * @brief Set equatorial sky direction + * + * @param[in] ra Right Ascension in radians. + * @param[in] dec Declination in radians. + ***************************************************************************/ +void GSkyDir::radec(const double& ra, const double& dec) +{ + // Set attributes + m_has_lb = 0; + m_has_radec = 1; + + // Set direction + m_ra = ra; + m_dec = dec; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set galactic sky direction + * + * @param[in] l Galactic longitude in radians. + * @param[in] b Galactic latitude in radians. + ***************************************************************************/ +void GSkyDir::lb(const double& l, const double& b) +{ + // Set attributes + m_has_lb = 1; + m_has_radec = 0; + + // Set direction + m_l = l; + m_b = b; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Returns galactic longitude in radians + ***************************************************************************/ +double GSkyDir::l(void) +{ + // If we have no galactic coordinates then get them now + if (!m_has_lb && m_has_radec) + equ2gal(); + + // Return galactic longitude + return m_l; +} + + +/***********************************************************************//** + * @brief Returns galactic latitude in radians + ***************************************************************************/ +double GSkyDir::b(void) +{ + // If we have no galactic coordinates then get them now + if (!m_has_lb && m_has_radec) + equ2gal(); + + // Return galactic latitude + return m_b; +} + + +/***********************************************************************//** + * @brief Returns Right Ascension in radians + ***************************************************************************/ +double GSkyDir::ra(void) +{ + // If we have no equatorial coordinates then get them now + if (!m_has_radec && m_has_lb) + gal2equ(); + + // Return Right Ascension + return m_ra; +} + + +/***********************************************************************//** + * @brief Returns Declination in radians + ***************************************************************************/ +double GSkyDir::dec(void) +{ + // If we have no equatorial coordinates then get them now + if (!m_has_radec && m_has_lb) + gal2equ(); + + // Return Declination + return m_dec; +} + /*========================================================================== = = @@ -123,37 +219,46 @@ GSkyDir& GSkyDir::operator= (const GSkyDir& dir) = = ==========================================================================*/ -/*************************************************************************** - * Initialise class members * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Initialise class members ***************************************************************************/ void GSkyDir::init_members(void) { // Initialise members + m_has_lb = 0; + m_has_radec = 0; + m_l = 0.0; + m_b = 0.0; + m_ra = 0.0; + m_dec = 0.0; // Return return; } -/*************************************************************************** - * Copy class members * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] dir Sky direction from which members should be copied ***************************************************************************/ void GSkyDir::copy_members(const GSkyDir& dir) { // Copy attributes - - // Copy other membres + m_has_lb = dir.m_has_lb; + m_has_radec = dir.m_has_radec; + m_l = dir.m_l; + m_b = dir.m_b; + m_ra = dir.m_ra; + m_dec = dir.m_dec; // Return return; } -/*************************************************************************** - * Delete class members * - * ----------------------------------------------------------------------- * +/***********************************************************************//** + * @brief Delete class members ***************************************************************************/ void GSkyDir::free_members(void) { @@ -166,12 +271,54 @@ void GSkyDir::free_members(void) } +/***********************************************************************//** + * @brief Convert equatorial to galactic coordinates + * + * @todo Needs to be implemented + ***************************************************************************/ +void GSkyDir::equ2gal(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Convert galactic to equatorial coordinates + * + * @todo Needs to be implemented + ***************************************************************************/ +void GSkyDir::gal2equ(void) +{ + // Return + return; +} + + /*========================================================================== = = = GSkyDir friends = = = ==========================================================================*/ +/***********************************************************************//** + * @brief Output operator + * + * @param[in] os Output stream + * @param[in] column Sky direction to put in output stream + ***************************************************************************/ +ostream& operator<< (ostream& os, const GSkyDir& dir) +{ + // Create coordinate system dependent output + if (dir.m_has_lb) + os << "(l,b)=(" << dir.m_l*rad2deg << "," << dir.m_b*rad2deg << ")"; + else if (dir.m_has_radec) + os << "(RA,Dec)=(" << dir.m_ra*rad2deg << "," << dir.m_dec*rad2deg << ")"; + + // Return output stream + return os; +} + /*========================================================================== = = diff --git a/test/test_python.py b/test/test_python.py index 963b068..35ce059 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -9,7 +9,7 @@ import os #================# def test_fits(): """ - Test GammaLib GFits interface. + Test GFits interface. """ # Remove test files try: @@ -120,7 +120,7 @@ def test_fits(): #=================# def test_node_array(): """ - Test GammaLib GNodeArray interface. + Test GNodeArray interface. """ # Set-up vector and data array vector = GVector(20) @@ -153,7 +153,7 @@ def test_node_array(): #===================# def test_lat_response(): """ - Test GammaLib GLATResponse interface. + Test GLATResponse interface. """ # Remove test file try: @@ -179,14 +179,40 @@ def test_lat_response(): #======================# def test_lat_observation(): """ - Test GammaLib GLATObservation interface. + Test GLATObservation interface. """ # Allocate LAT Observation obs = GLATObservation("data/FT1_253582800.fits.gz", "data/FT2_253582800.fits.gz") # - print obs.ft1() - print obs.ft2() + #print obs.ft1() + #print obs.ft2() + + +#==============# +# Test Healpix # +#==============# +def test_healpix(): + """ + Test GHealpix interface. + """ + # Open exposure cube + fits = GFits() + fits.open("data/expCube_253582800.fits.gz") + + # Lood Healpix table + pixels = GHealpix() + pixels.load(fits.hdu("EXPOSURE")) + + # Dump information + print pixels + + print pixels.pix2ang(0) + print " ",pixels.pix2ang(0).ra()*180/pi, pixels.pix2ang(0).dec()*180/pi + print pixels.pix2ang(1) + print pixels.pix2ang(2) + print pixels.pix2ang(3) + #==========================# # Main routine entry point # @@ -199,4 +225,5 @@ if __name__ == '__main__': test_node_array() test_lat_response() test_lat_observation() + test_healpix() \ No newline at end of file -- libgit2 0.21.4