enhance class

knodlseder
1 parent 86674206
Showing 3 changed files with 276 additions and 56 deletions
include/GHealpix.hpp
pyext/GHealpix.i
src/GHealpix.cpp
@@ -49,9 +49,11 @@ public:
     // Methods
     void    read(const GFitsHDU* hdu);
     int     nside(void) const;
+    int     order(void) const;
     int     num_pixels(void) const;
     double  omega(void) const;
     GSkyDir pix2ang(const int& ipix);
+    int     ang2pix(GSkyDir dir) const;
  
 private:
     // Private methods
@@ -59,15 +61,20 @@ private:
     void      copy_members(const GHealpix& pixels);
     void      free_members(void);
     GHealpix* clone(void) const;
+    int       nside2order(int nside);
     void      pix2xy(const int& ipix, int* x, int* y);
+    int       xy2pix(int x, int y) const;
     void      pix2ang_ring(int ipix, double* theta, double* phi);
     void      pix2ang_nest(int ipix, double* theta, double* phi);
+    int       ang2pix_z_phi_ring(double z, double phi) const;
+    int       ang2pix_z_phi_nest(double z, double phi) const;
  
     // Private data area
     int      m_nside;        //!< Number of divisions of each base pixel (1-8192)
     int      m_npface;       //!<
     int      m_ncap;         //!<
-    int      m_order;        //!< Ordering (0=ring, 1=nested, -1=?)
+    int      m_order;        //!< Order
+    int      m_scheme;       //!< Ordering scheme (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
@@ -40,6 +40,7 @@ public:
     int     num_pixels(void) const;
     double  omega(void) const;
     GSkyDir pix2ang(const int& ipix);
+    int     ang2pix(GSkyDir dir) const;
 };
  
  
@@ -19,7 +19,8 @@
  
 /* __ Includes ___________________________________________________________ */
 #include <iostream>
-#include <math.h>
+//#include <math.h>
+#include <cmath>
 #include "GException.hpp"
 #include "GTools.hpp"
 #include "GHealpix.hpp"
@@ -37,15 +38,19 @@
  
 /* __ Debug definitions __________________________________________________ */
  
+/* __ Local prototypes ___________________________________________________ */
+unsigned int isqrt(unsigned int arg);
+
 /* __ 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};
+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};
+const int order_max = 13;
+const int ns_max    = 1 << order_max;
  
 /* __ Static conversion arrays ___________________________________________ */
 static short ctab[0x100];
 static short utab[0x100];
  
-
 /*==========================================================================
  =                                                                         =
  =                     GHealpix constructors/destructors                   =
@@ -178,6 +183,8 @@ void GHealpix::read(const GFitsHDU* hdu)
     m_fact2      = 4.0 / m_num_pixels;
     m_fact1      = 2 * m_nside * m_fact2;
     m_omega      = fourpi / m_num_pixels;
+    m_order      = nside2order(m_nside);
+
  
     // Get ordering scheme from ORDERING keyword
     std::string ordering;
@@ -185,14 +192,14 @@ void GHealpix::read(const GFitsHDU* hdu)
         ordering = hdu->card("ORDERING")->string();
     }
     catch (GException::fits_key_not_found &e) {
-        m_order = -1; // Flag unknown ordering
+        m_scheme = -1; // Flag unknown ordering
     }
  
     // Decode ordering scheme string
     if (ordering == "RING")
-        m_order = 0;
+        m_scheme = 0;
     else if (ordering == "NESTED")
-        m_order = 1;
+        m_scheme = 1;
     else
         throw GException::healpix(G_LOAD, "Invalid ordering");
  
@@ -268,6 +275,16 @@ int GHealpix::nside(void) const
  
  
 /***********************************************************************//**
+ * @brief Returns ordering parameter.
+ ***************************************************************************/
+int GHealpix::order(void) const
+{
+    // Return order
+    return m_order;
+}
+
+
+/***********************************************************************//**
  * @brief Returns number of pixels.
  ***************************************************************************/
 int GHealpix::num_pixels(void) const
@@ -300,7 +317,7 @@ GSkyDir GHealpix::pix2ang(const int&amp; ipix)
     double  phi   = 0.0;
  
     // Perform ordering dependent conversion
-    switch (m_order) {
+    switch (m_scheme) {
     case 0:
         pix2ang_ring(ipix, &theta, &phi);
         break;
@@ -328,6 +345,49 @@ GSkyDir GHealpix::pix2ang(const int&amp; ipix)
 }
  
  
+/***********************************************************************//**
+ * @brief Returns pixel for a given sky direction
+ *
+ * @param[in] ipix Pixel number (0,1,...,num_pixels)
+ ***************************************************************************/
+int GHealpix::ang2pix(GSkyDir dir) const
+{
+    // Declare result
+    int    ipix;
+    double z;
+    double phi;
+
+    // Compute coordinate system dependent (z,phi)
+    switch (m_coordsys) {
+    case 0:
+        z   = cos(pihalf-dir.dec());
+        phi = dir.ra();
+        break;
+    case 1:
+        z   = cos(pihalf-dir.b());
+        phi = dir.l();
+        break;
+    default:
+        break;
+    }
+    
+    // Perform ordering dependent conversion
+    switch (m_scheme) {
+    case 0:
+        ipix = ang2pix_z_phi_ring(z, phi);
+        break;
+    case 1:
+        ipix = ang2pix_z_phi_nest(z, phi);
+        break;
+    default:
+        break;
+    }
+    
+    // Return pixel index
+    return ipix;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                         GHealpix private methods                        =
@@ -344,6 +404,7 @@ void GHealpix::init_members(void)
     m_npface      = 0;
     m_ncap        = 0;
     m_order       = 0;
+    m_scheme      = 0;
     m_coordsys    = 0;
     m_num_pixels  = 0;
     m_size_pixels = 0;
@@ -380,6 +441,7 @@ void GHealpix::copy_members(const GHealpix&amp; pixels)
     m_npface      = pixels.m_npface;
     m_ncap        = pixels.m_ncap;
     m_order       = pixels.m_order;
+    m_scheme      = pixels.m_scheme;
     m_coordsys    = pixels.m_coordsys;
     m_num_pixels  = pixels.m_num_pixels;
     m_size_pixels = pixels.m_size_pixels;
@@ -434,6 +496,32 @@ GHealpix* GHealpix::clone(void) const
  
  
 /***********************************************************************//**
+ * @brief Convert nside to order
+ *
+ * @param[in] nside Number of sides.
+ ***************************************************************************/
+int GHealpix::nside2order(int nside)
+{
+    // Initialise order
+    int order = -1;
+    
+    // Determine order
+    for (int m = 0; m <= order_max; ++m) {
+        int nstest = 1 << m;
+        if (nside == nstest) {
+            order = m;
+            break;
+        }
+        if (nside < nstest)
+            break;
+    }
+    
+    // Return order
+    return order;
+}
+
+
+/***********************************************************************//**
  * @brief Convert pixel index to (x,y) coordinate
  *
  * @param[in] ipix Pixel index for which (x,y) are to be computed.
@@ -456,6 +544,19 @@ void GHealpix::pix2xy(const int&amp; ipix, int* x, int* y)
  
  
 /***********************************************************************//**
+ * @brief Convert (x,y) coordinate to pixel
+ *
+ * @param[in] x x coordinate.
+ * @param[in] y y coordinate.
+ ***************************************************************************/
+int GHealpix::xy2pix(int x, int y) const
+{
+    // Return pixel
+    return utab[x&0xff] | (utab[x>>8]<<16) | (utab[y&0xff]<<1) | (utab[y>>8]<<17);
+}
+
+
+/***********************************************************************//**
  * @brief Convert pixel index to (theta,phi) angles for ring ordering
  *
  * @param[in] ipix Pixel index for which (theta,phi) are to be computed.
@@ -470,48 +571,34 @@ void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi)
     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);
+    if (ipix < m_ncap) {
+        int iring = int(0.5*(1+isqrt(1+2*ipix))); // counted from North pole
+        int iphi  = (ipix+1) - 2*iring*(iring-1);
+        *theta    = acos(1.0 - iring*iring / m_fact2);
+        *phi      = (iphi - 0.5) * pi/(2.0*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);
+    else if (ipix < (m_num_pixels - m_ncap)) {
+        int    ip    = ipix - m_ncap;
+        int    iring = ip/(4*m_nside) + m_nside;   // counted from North pole
+        int    iphi  = ip%(4*m_nside) + 1;
+        double fodd  = ((iring+m_nside)&1) ? 1 : 0.5;
+        int    nl2   = 2*m_nside;
+        *theta       = acos((nl2 - iring) / m_fact1);
+        *phi         = (iphi - fodd) * pi/nl2;
     }
  
-    // Handle South Polar cap
+    // 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);
+        int ip    = m_num_pixels - ipix;
+        int iring = int(0.5*(1+isqrt(2*ip-1)));    // Counted from South pole
+        int iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
+        *theta    = acos(iring*iring / m_fact2 - 1.0);
+        *phi      = (iphi - 0.5) * pi/(2.*iring);
     }
-    
+
     // Return
     return;
 }
@@ -534,7 +621,7 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi)
  
     // Get face number and index in face
     int nl4      = 4 * m_nside;
-    int face_num = ipix/m_npface;            // Face number in {0,11}
+    int face_num = ipix >> (2*m_order);      // Face number in {0,11}
     int ipf      = ipix & (m_npface - 1);
  
     // Get pixel coordinates
@@ -543,7 +630,7 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi)
     pix2xy(ipf, &ix, &iy);
  
     // Computes the z coordinate on the sphere
-    int jr = jrll[face_num]*m_nside - ix - iy - 1;
+    int jr = (jrll[face_num] << m_order) - ix - iy - 1;
  
     // Declare result variables
     int    nr;
@@ -585,6 +672,118 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi)
 }
  
  
+/***********************************************************************//**
+ * @brief Returns pixels which contains angular coordinates (z,phi)
+ *
+ * @param[in] z Cosine of zenith angle - cos(theta).
+ * @param[in] phi Azimuth angle in radians.
+ ***************************************************************************/
+int GHealpix::ang2pix_z_phi_ring(double z, double phi) const
+{
+    // Initialise pixel
+    int ipix = 0;
+    
+    // Setup
+    double za = abs(z);
+    double tt = modulo(phi,twopi) * inv_pihalf; // in [0,4)
+
+    // Equatorial region
+    if (za <= twothird) {
+        double temp1  = m_nside*(0.5+tt);
+        double temp2  = m_nside*z*0.75;
+        int    jp     = int(temp1-temp2);           // index of ascending edge line
+        int    jm     = int(temp1+temp2);           // index of descending edge line
+        int    ir     = m_nside + 1 + jp - jm;      // in {1,2n+1}
+        int    kshift = 1 - (ir & 1);               // kshift=1 if ir even, 0 otherwise
+        int    ip     = (jp+jm-m_nside+kshift+1)/2; // in {0,4n-1}
+        ip            = modulo(ip,4*m_nside);
+        ipix          = m_ncap + (ir-1)*4*m_nside + ip;
+    }
+    
+    // North & South polar caps
+    else {
+        double tp  = tt - int(tt);
+        double tmp = m_nside * sqrt(3*(1-za));
+        int    jp  = int(tp*tmp);       // increasing edge line index
+        int    jm  = int((1.0-tp)*tmp); // decreasing edge line index
+        int    ir  = jp + jm + 1;       // ring number counted from the closest pole
+        int    ip  = int(tt*ir);        // in {0,4*ir-1}
+        ip = modulo(ip,4*ir);
+        if (z>0)
+            ipix = 2*ir*(ir-1) + ip;
+        else
+            ipix = m_num_pixels - 2*ir*(ir+1) + ip;
+    }
+    
+    // Return pixel
+    return ipix;
+}
+    
+    
+/***********************************************************************//**
+ * @brief Returns pixels which contains angular coordinates (z,phi)
+ *
+ * @param[in] z Cosine of zenith angle - cos(theta).
+ * @param[in] phi Azimuth angle in radians.
+ ***************************************************************************/
+int GHealpix::ang2pix_z_phi_nest(double z, double phi) const
+{
+    // Initialise face and pixel numbers
+    int face_num;
+    int ix;
+    int iy;
+    
+    // Setup
+    double za = abs(z);
+    double tt = modulo(phi,twopi) * inv_pihalf; // in [0,4)
+
+    // Equatorial region
+    if (za <= twothird) {
+        double temp1 = ns_max*(0.5+tt);
+        double temp2 = ns_max*z*0.75;
+        int    jp    = int(temp1-temp2); // index of  ascending edge line
+        int    jm    = int(temp1+temp2); // index of descending edge line
+        int    ifp   = jp >> order_max;  // in {0,4}
+        int    ifm   = jm >> order_max;
+        if (ifp == ifm)                  // faces 4 to 7
+            face_num = (ifp==4) ? 4: ifp+4;
+        else if (ifp < ifm)              // (half-)faces 0 to 3
+            face_num = ifp;
+        else                             // (half-)faces 8 to 11
+            face_num = ifm + 8;
+        ix = jm & (ns_max-1);
+        iy = ns_max - (jp & (ns_max-1)) - 1;
+    }
+    
+    // Polar region, za > 2/3
+    else {
+        int    ntt = int(tt);
+        double tp  = tt-ntt;
+        double tmp = ns_max*sqrt(3*(1-za));
+        int    jp  = int(tp*tmp);        // increasing edge line index
+        int    jm  = int((1.0-tp)*tmp);  // decreasing edge line index
+        if (jp >= ns_max) jp = ns_max-1; // for points too close to the boundary
+        if (jm >= ns_max) jm = ns_max-1;
+        if (z >= 0) {
+            face_num = ntt;              // in {0,3}
+            ix       = ns_max - jm - 1;
+            iy       = ns_max - jp - 1;
+        }
+        else {
+            face_num = ntt + 8;          // in {8,11}
+            ix       =  jp;
+            iy       =  jm;
+        }
+    }
+    
+    // Get pixel
+    int ipf = xy2pix(ix, iy);
+    ipf >>= (2*(order_max - m_order));     // in {0, nside**2 - 1}
+    return ipf + (face_num<<(2*m_order));  // in {0, 12*nside**2 - 1}
+}
+
+
+
 /*==========================================================================
  =                                                                         =
  =                             GHealpix friends                            =
@@ -601,15 +800,18 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GHealpix&amp; 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 << " Nside (number of divisions): " << pixels.m_nside << endl;
+    os << " Npface (pixels per face) ..: " << pixels.m_npface << endl;
+    os << " Ncap (number of cap pixels): " << pixels.m_ncap << endl;
+    os << " Npix (number of pixels) ...: " << pixels.m_num_pixels << endl;
+    os << " Order .....................: " << pixels.m_order << endl;
     os << " Pixel vector size .........: " << pixels.m_size_pixels << endl;
     os << " Solid angle ...............: " << scientific << pixels.m_omega 
-       << fixed << endl;
+       << fixed << " sr" << endl;
  
     // Put ordering in stream
     os << " Ordering ..................: ";
-    switch (pixels.m_order) {
+    switch (pixels.m_scheme) {
     case 0:
         os << "Ring" << endl;
         break;
@@ -628,22 +830,19 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GHealpix&amp; pixels)
     os << " Coordinate system .........: ";
     switch (pixels.m_coordsys) {
     case 0:
-        os << "Equatorial (RA,Dec)" << endl;
+        os << "Equatorial (RA,Dec)";
         break;
     case 1:
-        os << "Galactic (l,b)" << endl;
+        os << "Galactic (l,b)";
         break;
     case -1:
-        os << "*** Unknown ***" << endl;
+        os << "*** Unknown ***";
         break;
     default:
-        os << "*** Invalid ***" << endl;
+        os << "*** Invalid ***";
         break;
     }
  
-    for (int i = 0; i < pixels.m_num_pixels; ++i)
-        os << pixels.m_dir[i] << " ";
-
     // Return output stream
     return os;
 }
@@ -654,3 +853,16 @@ ostream&amp; operator&lt;&lt; (ostream&amp; os, const GHealpix&amp; pixels)
  =                     Other functions used by GHealpix                    =
  =                                                                         =
  ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Integer n that fulfills n*n <= arg < (n+1)*(n+1)
+ *
+ * @param[in] arg Argument.
+ *
+ * Returns the integer \a n, which fulfills \a n*n <= arg < (n+1)*(n+1).
+ ***************************************************************************/
+unsigned int isqrt(unsigned int arg)
+{
+    // Return
+    return unsigned(sqrt(arg+0.5));
+}