healpix_base.cc

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Copyright (C) 2003, 2004, 2005, 2006 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include "healpix_base.h"
00033 #include "cxxutils.h"
00034 #include "pointing.h"
00035 #include "arr.h"
00036 #include "geom_utils.h"
00037 
00038 using namespace std;
00039 
00040 short Healpix_Base::ctab[];
00041 short Healpix_Base::utab[];
00042 
00043 const nside_dummy SET_NSIDE=nside_dummy();
00044 
00045 Healpix_Base::Tablefiller::Tablefiller()
00046   {
00047   for (int m=0; m<0x100; ++m)
00048     {
00049     ctab[m] =
00050          (m&0x1 )       | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
00051       | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4);
00052     utab[m] =
00053          (m&0x1 )       | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
00054       | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7);
00055     }
00056   }
00057 
00058 Healpix_Base::Tablefiller Healpix_Base::Filler;
00059 
00060 const int Healpix_Base::jrll[] = { 2,2,2,2,3,3,3,3,4,4,4,4 };
00061 const int Healpix_Base::jpll[] = { 1,3,5,7,0,2,4,6,1,3,5,7 };
00062 
00063 int Healpix_Base::npix2nside (int npix)
00064   {
00065   int res=isqrt(npix/12);
00066   planck_assert (npix==res*res*12, "npix2nside: invalid argument");
00067   return res;
00068   }
00069 
00070 int Healpix_Base::ring_above (double z) const
00071   {
00072   double az=abs(z);
00073   if (az>twothird) // polar caps
00074     {
00075     int iring = int(nside_*sqrt(3*(1-az)));
00076     return (z>0) ? iring : 4*nside_-iring-1;
00077     }
00078   else // ----- equatorial region ---------
00079     return int(nside_*(2-1.5*z));
00080   }
00081 
00082 void Healpix_Base::in_ring(int iz, double phi0, double dphi,
00083   vector<int> &listir) const
00084   {
00085   int nr, ir, ipix1;
00086   double shift=0.5;
00087 
00088   if (iz<nside_) // north pole
00089     {
00090     ir = iz;
00091     nr = ir*4;
00092     ipix1 = 2*ir*(ir-1);        //    lowest pixel number in the ring
00093     }
00094   else if (iz>(3*nside_)) // south pole
00095     {
00096     ir = 4*nside_ - iz;
00097     nr = ir*4;
00098     ipix1 = npix_ - 2*ir*(ir+1); // lowest pixel number in the ring
00099     }
00100   else // equatorial region
00101     {
00102     ir = iz - nside_ + 1;           //    within {1, 2*nside + 1}
00103     nr = nside_*4;
00104     if ((ir&1)==0) shift = 0;
00105     ipix1 = ncap_ + (ir-1)*nr; // lowest pixel number in the ring
00106     }
00107 
00108   int ipix2 = ipix1 + nr - 1;       //    highest pixel number in the ring
00109 
00110    // ----------- constructs the pixel list --------------
00111   if (dphi > (pi-1e-7))
00112     for (int i=ipix1; i<=ipix2; ++i) listir.push_back(i);
00113   else
00114     {
00115     int ip_lo = ifloor<int>(nr*inv_twopi*(phi0-dphi) - shift)+1;
00116     int ip_hi = ifloor<int>(nr*inv_twopi*(phi0+dphi) - shift);
00117     int pixnum = ip_lo+ipix1;
00118     if (pixnum<ipix1) pixnum += nr;
00119     for (int i=ip_lo; i<=ip_hi; ++i, ++pixnum)
00120       {
00121       if (pixnum>ipix2) pixnum -= nr;
00122       listir.push_back(pixnum);
00123       }
00124     }
00125   }
00126 
00127 void Healpix_Base::nest2xyf (int pix, int &ix, int &iy, int &face_num) const
00128   {
00129   face_num = pix>>(2*order_);
00130   pix &= (npface_-1);
00131   int raw = (pix&0x5555) | ((pix&0x55550000)>>15);
00132   ix = ctab[raw&0xff] | (ctab[raw>>8]<<4);
00133   pix >>= 1;
00134   raw = (pix&0x5555) | ((pix&0x55550000)>>15);
00135   iy = ctab[raw&0xff] | (ctab[raw>>8]<<4);
00136   }
00137 
00138 int Healpix_Base::xyf2nest (int ix, int iy, int face_num) const
00139   {
00140   return (face_num<<(2*order_)) +
00141       (utab[ix&0xff] | (utab[ix>>8]<<16)
00142     | (utab[iy&0xff]<<1) | (utab[iy>>8]<<17));
00143   }
00144 
00145 void Healpix_Base::ring2xyf (int pix, int &ix, int &iy, int &face_num) const
00146   {
00147   int iring, iphi, kshift, nr;
00148 
00149   int nl2 = 2*nside_;
00150 
00151   if (pix<ncap_) // North Polar cap
00152     {
00153     iring = int(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00154     iphi  = (pix+1) - 2*iring*(iring-1);
00155     kshift = 0;
00156     nr = iring;
00157     face_num=0;
00158     int tmp = iphi-1;
00159     if (tmp>=(2*iring))
00160       {
00161       face_num=2;
00162       tmp-=2*iring;
00163       }
00164     if (tmp>=iring) ++face_num;
00165     }
00166   else if (pix<(npix_-ncap_)) // Equatorial region
00167     {
00168     int ip = pix - ncap_;
00169     if (order_>=0)
00170       {
00171       iring = (ip>>(order_+2)) + nside_; // counted from North pole
00172       iphi  = (ip&(4*nside_-1)) + 1;
00173       }
00174     else
00175       {
00176       iring = (ip/(4*nside_)) + nside_; // counted from North pole
00177       iphi  = (ip%(4*nside_)) + 1;
00178       }
00179     kshift = (iring+nside_)&1;
00180     nr = nside_;
00181     unsigned int ire = iring-nside_+1;
00182     unsigned int irm = nl2+2-ire;
00183     int ifm, ifp;
00184     if (order_>=0)
00185       {
00186       ifm = (iphi - ire/2 + nside_ -1) >> order_;
00187       ifp = (iphi - irm/2 + nside_ -1) >> order_;
00188       }
00189     else
00190       {
00191       ifm = (iphi - ire/2 + nside_ -1) / nside_;
00192       ifp = (iphi - irm/2 + nside_ -1) / nside_;
00193       }
00194     if (ifp == ifm) // faces 4 to 7
00195       face_num = (ifp==4) ? 4 : ifp+4;
00196     else if (ifp<ifm) // (half-)faces 0 to 3
00197       face_num = ifp;
00198     else // (half-)faces 8 to 11
00199       face_num = ifm + 8;
00200     }
00201   else // South Polar cap
00202     {
00203     int ip = npix_ - pix;
00204     iring = int(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00205     iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
00206     kshift = 0;
00207     nr = iring;
00208     iring = 2*nl2-iring;
00209     face_num=8;
00210     int tmp = iphi-1;
00211     if (tmp>=(2*nr))
00212       {
00213       face_num=10;
00214       tmp-=2*nr;
00215       }
00216     if (tmp>=nr) ++face_num;
00217     }
00218 
00219   int irt = iring - (jrll[face_num]*nside_) + 1;
00220   int ipt = 2*iphi- jpll[face_num]*nr - kshift -1;
00221   if (ipt>=nl2) ipt-=8*nside_;
00222 
00223   ix =  (ipt-irt) >>1;
00224   iy =(-(ipt+irt))>>1;
00225   }
00226 
00227 int Healpix_Base::xyf2ring (int ix, int iy, int face_num) const
00228   {
00229   int nl4 = 4*nside_;
00230   int jr = (jrll[face_num]*nside_) - ix - iy  - 1;
00231 
00232   int nr, kshift, n_before;
00233   if (jr<nside_)
00234     {
00235     nr = jr;
00236     n_before = 2*nr*(nr-1);
00237     kshift = 0;
00238     }
00239   else if (jr > 3*nside_)
00240     {
00241     nr = nl4-jr;
00242     n_before = npix_ - 2*(nr+1)*nr;
00243     kshift = 0;
00244     }
00245   else
00246     {
00247     nr = nside_;
00248     n_before = ncap_ + (jr-nside_)*nl4;
00249     kshift = (jr-nside_)&1;
00250     }
00251 
00252   int jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
00253   if (jp>nl4)
00254     jp-=nl4;
00255   else
00256     if (jp<1) jp+=nl4;
00257 
00258   return n_before + jp - 1;
00259   }
00260 
00261 double Healpix_Base::ring2z (int ring) const
00262   {
00263   if (ring<nside_)
00264     return 1 - ring*ring*fact2_;
00265   if (ring <=3*nside_)
00266     return (2*nside_-ring)*fact1_;
00267   ring=4*nside_ - ring;
00268   return ring*ring*fact2_ - 1;
00269   }
00270 
00271 int Healpix_Base::pix2ring (int pix) const
00272   {
00273   if (scheme_==RING)
00274     {
00275     if (pix<ncap_) // North Polar cap
00276       return int(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00277     else if (pix<(npix_-ncap_)) // Equatorial region
00278       {
00279       int ip  = pix - ncap_;
00280       return ip/(4*nside_) + nside_; // counted from North pole
00281       }
00282     else // South Polar cap
00283       {
00284       int ip = npix_ - pix;
00285       return 4*nside_ - int(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00286       }
00287     }
00288   else
00289     {
00290     int face_num, ix, iy;
00291     nest2xyf(pix,ix,iy,face_num);
00292     return (jrll[face_num]<<order_) - ix - iy - 1;
00293     }
00294   }
00295 
00296 int Healpix_Base::nest2ring (int pix) const
00297   {
00298   planck_assert(order_>=0, "nest2ring: need hierarchical map");
00299   int ix, iy, face_num;
00300   nest2xyf (pix, ix, iy, face_num);
00301   return xyf2ring (ix, iy, face_num);
00302   }
00303 
00304 int Healpix_Base::ring2nest (int pix) const
00305   {
00306   planck_assert(order_>=0, "ring2nest: need hierarchical map");
00307   int ix, iy, face_num;
00308   ring2xyf (pix, ix, iy, face_num);
00309   return xyf2nest (ix, iy, face_num);
00310   }
00311 
00312 int Healpix_Base::nest2peano (int pix) const
00313   {
00314   static const unsigned char subpix[8][4] = {
00315     { 0, 1, 3, 2 }, { 3, 0, 2, 1 }, { 2, 3, 1, 0 }, { 1, 2, 0, 3 },
00316     { 0, 3, 1, 2 }, { 1, 0, 2, 3 }, { 2, 1, 3, 0 }, { 3, 2, 0, 1 } };
00317   const unsigned char subpath[8][4] = {
00318     { 4, 0, 6, 0 }, { 7, 5, 1, 1 }, { 2, 4, 2, 6 }, { 3, 3, 7, 5 },
00319     { 0, 2, 4, 4 }, { 5, 1, 5, 3 }, { 6, 6, 0, 2 }, { 1, 7, 3, 7 } };
00320   static const unsigned char face2path[12] = {
00321     2, 5, 2, 5, 3, 6, 3, 6, 2, 3, 2, 3 };
00322   static const unsigned char face2peanoface[12] = {
00323     0, 5, 6, 11, 10, 1, 4, 7, 2, 3, 8, 9 };
00324 
00325   int face = pix>>(2*order_);
00326   unsigned char path = face2path[face];
00327   int result = 0;
00328 
00329   for (int shift=2*order_-2; shift>=0; shift-=2)
00330     {
00331     unsigned char spix = (pix>>shift) & 0x3;
00332     result <<= 2;
00333     result |= subpix[path][spix];
00334     path=subpath[path][spix];
00335     }
00336 
00337   return result + (int(face2peanoface[face])<<(2*order_));
00338   }
00339 
00340 int Healpix_Base::peano2nest (int pix) const
00341   {
00342   static const unsigned char subpix[8][4] = {
00343     { 0, 1, 3, 2 }, { 1, 3, 2, 0 }, { 3, 2, 0, 1 }, { 2, 0, 1, 3 },
00344     { 0, 2, 3, 1 }, { 1, 0, 2, 3 }, { 3, 1, 0, 2 }, { 2, 3, 1, 0 } };
00345   static const unsigned char subpath[8][4] = {
00346     { 4, 0, 0, 6 }, { 5, 1, 1, 7 }, { 6, 2, 2, 4 }, { 7, 3, 3, 5 },
00347     { 0, 4, 4, 2 }, { 1, 5, 5, 3 }, { 2, 6, 6, 0 }, { 3, 7, 7, 1 } };
00348   static const unsigned char face2path[12] = {
00349     2, 6, 2, 3, 3, 5, 2, 6, 2, 3, 3, 5 };
00350   static const unsigned char peanoface2face[12] = {
00351     0, 5, 8, 9, 6, 1, 2, 7, 10, 11, 4, 3 };
00352 
00353   int face = pix>>(2*order_);
00354   unsigned char path = face2path[face];
00355   int result = 0;
00356 
00357   for (int shift=2*order_-2; shift>=0; shift-=2)
00358     {
00359     unsigned char spix = (pix>>shift) & 0x3;
00360     result <<= 2;
00361     result |= subpix[path][spix];
00362     path=subpath[path][spix];
00363     }
00364 
00365   return result + (int(peanoface2face[face])<<(2*order_));
00366   }
00367 
00368 int Healpix_Base::ang2pix_z_phi (double z, double phi) const
00369   {
00370   double za = abs(z);
00371   double tt = fmodulo(phi,twopi) * inv_halfpi; // in [0,4)
00372 
00373   if (scheme_==RING)
00374     {
00375     if (za<=twothird) // Equatorial region
00376       {
00377       double temp1 = nside_*(0.5+tt);
00378       double temp2 = nside_*z*0.75;
00379       int jp = int(temp1-temp2); // index of  ascending edge line
00380       int jm = int(temp1+temp2); // index of descending edge line
00381 
00382       // ring number counted from z=2/3
00383       int ir = nside_ + 1 + jp - jm; // in {1,2n+1}
00384       int kshift = 1-(ir&1); // kshift=1 if ir even, 0 otherwise
00385 
00386       int ip = (jp+jm-nside_+kshift+1)/2; // in {0,4n-1}
00387       ip = imodulo(ip,4*nside_);
00388 
00389       return ncap_ + (ir-1)*4*nside_ + ip;
00390       }
00391     else  // North & South polar caps
00392       {
00393       double tp = tt-int(tt);
00394       double tmp = nside_*sqrt(3*(1-za));
00395 
00396       int jp = int(tp*tmp); // increasing edge line index
00397       int jm = int((1.0-tp)*tmp); // decreasing edge line index
00398 
00399       int ir = jp+jm+1; // ring number counted from the closest pole
00400       int ip = int(tt*ir); // in {0,4*ir-1}
00401       ip = imodulo(ip,4*ir);
00402 
00403       if (z>0)
00404         return 2*ir*(ir-1) + ip;
00405       else
00406         return npix_ - 2*ir*(ir+1) + ip;
00407       }
00408     }
00409   else // scheme_ == NEST
00410     {
00411     int face_num, ix, iy;
00412 
00413     if (za<=twothird) // Equatorial region
00414       {
00415       double temp1 = nside_*(0.5+tt);
00416       double temp2 = nside_*(z*0.75);
00417       int jp = int(temp1-temp2); // index of  ascending edge line
00418       int jm = int(temp1+temp2); // index of descending edge line
00419       int ifp = jp >> order_;  // in {0,4}
00420       int ifm = jm >> order_;
00421       if (ifp == ifm)           // faces 4 to 7
00422         face_num = (ifp==4) ? 4: ifp+4;
00423       else if (ifp < ifm)       // (half-)faces 0 to 3
00424         face_num = ifp;
00425       else                      // (half-)faces 8 to 11
00426         face_num = ifm + 8;
00427 
00428       ix = jm & (nside_-1);
00429       iy = nside_ - (jp & (nside_-1)) - 1;
00430       }
00431     else // polar region, za > 2/3
00432       {
00433       int ntt = int(tt);
00434       if (ntt>=4) ntt=3;
00435       double tp = tt-ntt;
00436       double tmp = nside_*sqrt(3*(1-za));
00437 
00438       int jp = int(tp*tmp); // increasing edge line index
00439       int jm = int((1.0-tp)*tmp); // decreasing edge line index
00440       if (jp>=nside_) jp = nside_-1; // for points too close to the boundary
00441       if (jm>=nside_) jm = nside_-1;
00442       if (z >= 0)
00443         {
00444         face_num = ntt;  // in {0,3}
00445         ix = nside_ - jm - 1;
00446         iy = nside_ - jp - 1;
00447         }
00448       else
00449         {
00450         face_num = ntt + 8; // in {8,11}
00451         ix =  jp;
00452         iy =  jm;
00453         }
00454       }
00455 
00456     return xyf2nest(ix,iy,face_num);
00457     }
00458   }
00459 
00460 void Healpix_Base::pix2ang_z_phi (int pix, double &z, double &phi) const
00461   {
00462   if (scheme_==RING)
00463     {
00464     if (pix<ncap_) // North Polar cap
00465       {
00466       int iring = int(0.5*(1+isqrt(1+2*pix))); //counted from North pole
00467       int iphi  = (pix+1) - 2*iring*(iring-1);
00468 
00469       z = 1.0 - (iring*iring)*fact2_;
00470       phi = (iphi-0.5) * halfpi/iring;
00471       }
00472     else if (pix<(npix_-ncap_)) // Equatorial region
00473       {
00474       int ip  = pix - ncap_;
00475       int iring = ip/(4*nside_) + nside_; // counted from North pole
00476       int iphi  = ip%(4*nside_) + 1;
00477       // 1 if iring+nside is odd, 1/2 otherwise
00478       double fodd = ((iring+nside_)&1) ? 1 : 0.5;
00479 
00480       int nl2 = 2*nside_;
00481       z = (nl2-iring)*fact1_;
00482       phi = (iphi-fodd) * pi/nl2;
00483       }
00484     else // South Polar cap
00485       {
00486       int ip = npix_ - pix;
00487       int iring = int(0.5*(1+isqrt(2*ip-1))); //counted from South pole
00488       int iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
00489 
00490       z = -1.0 + (iring*iring)*fact2_;
00491       phi = (iphi-0.5) * halfpi/iring;
00492       }
00493     }
00494   else
00495     {
00496     int nl4 = nside_*4;
00497 
00498     int face_num, ix, iy;
00499     nest2xyf(pix,ix,iy,face_num);
00500 
00501     int jr = (jrll[face_num]<<order_) - ix - iy - 1;
00502 
00503     int nr, kshift;
00504     if (jr<nside_)
00505       {
00506       nr = jr;
00507       z = 1 - nr*nr*fact2_;
00508       kshift = 0;
00509       }
00510     else if (jr > 3*nside_)
00511       {
00512       nr = nl4-jr;
00513       z = nr*nr*fact2_ - 1;
00514       kshift = 0;
00515       }
00516     else
00517       {
00518       nr = nside_;
00519       z = (2*nside_-jr)*fact1_;
00520       kshift = (jr-nside_)&1;
00521       }
00522 
00523     int jp = (jpll[face_num]*nr + ix -iy + 1 + kshift) / 2;
00524     if (jp>nl4) jp-=nl4;
00525     if (jp<1) jp+=nl4;
00526 
00527     phi = (jp-(kshift+1)*0.5)*(halfpi/nr);
00528     }
00529   }
00530 
00531 void Healpix_Base::query_disc (const pointing &ptg, double radius,
00532   vector<int>& listpix) const
00533   {
00534   listpix.clear();
00535 
00536   double dth1 = fact2_;
00537   double dth2 = fact1_;
00538   double cosang = cos(radius);
00539 
00540   double z0 = cos(ptg.theta);
00541   double xa = 1./sqrt((1-z0)*(1+z0));
00542 
00543   double rlat1  = ptg.theta - radius;
00544   double zmax = cos(rlat1);
00545   int irmin = ring_above (zmax)+1;
00546 
00547   if (rlat1<=0) // north pole in the disc
00548     for (int m=1; m<irmin; ++m) // rings completely in the disc
00549       in_ring (m, 0, pi, listpix);
00550 
00551   double rlat2  = ptg.theta + radius;
00552   double zmin = cos(rlat2);
00553   int irmax = ring_above (zmin);
00554 
00555 // ------------- loop on ring number ---------------------
00556   for (int iz=irmin; iz<=irmax; ++iz) // rings partially in the disc
00557     {
00558     double z;
00559     if (iz<nside_) // north polar cap
00560       z = 1.0 - iz*iz*dth1;
00561     else if (iz <= (3*nside_)) // tropical band + equat.
00562       z = (2*nside_-iz) * dth2;
00563     else
00564       z = -1.0 + (4*nside_-iz)*(4*nside_-iz)*dth1;
00565 
00566 // --------- phi range in the disc for each z ---------
00567     double x = (cosang-z*z0)*xa;
00568     double ysq = 1-z*z-x*x;
00569     planck_assert(ysq>=0, "error in query_disc()");
00570     double dphi=atan2(sqrt(ysq),x);
00571     in_ring (iz, ptg.phi, dphi, listpix);
00572     }
00573 
00574   if (rlat2>=pi) // south pole in the disc
00575     for (int m=irmax+1; m<(4*nside_); ++m)  // rings completely in the disc
00576       in_ring (m, 0, pi, listpix);
00577 
00578   if (scheme_==NEST)
00579     for (unsigned int m=0; m<listpix.size(); ++m)
00580       listpix[m] = ring2nest(listpix[m]);
00581   }
00582 
00583 void Healpix_Base::get_ring_info (int ring, int &startpix, int &ringpix,
00584   double &costheta, double &sintheta, bool &shifted) const
00585   {
00586   planck_assert(scheme_==RING,"map must be in RING scheme");
00587   int northring = (ring>2*nside_) ? 4*nside_-ring : ring;
00588   if (northring < nside_)
00589     {
00590     double tmp = northring*northring*fact2_;
00591     costheta = 1 - tmp;
00592     sintheta = sqrt(tmp*(2-tmp));
00593     ringpix = 4*northring;
00594     shifted = true;
00595     startpix = 2*northring*(northring-1);
00596     }
00597   else
00598     {
00599     costheta = (2*nside_-northring)*fact1_;
00600     sintheta = sqrt((1+costheta)*(1-costheta));
00601     ringpix = 4*nside_;
00602     shifted = ((northring-nside_) & 1) == 0;
00603     startpix = ncap_ + (northring-nside_)*ringpix;
00604     }
00605   if (northring != ring) // southern hemisphere
00606     {
00607     costheta = -costheta;
00608     startpix = npix_ - startpix - ringpix;
00609     }
00610   }
00611 
00612 void Healpix_Base::neighbors (int pix, fix_arr<int,8> &result) const
00613   {
00614   static const int xoffset[] = { -1,-1, 0, 1, 1, 1, 0,-1 };
00615   static const int yoffset[] = {  0, 1, 1, 1, 0,-1,-1,-1 };
00616   static const int facearray[][12] =
00617         { {  8, 9,10,11,-1,-1,-1,-1,10,11, 8, 9 },   // S
00618           {  5, 6, 7, 4, 8, 9,10,11, 9,10,11, 8 },   // SE
00619           { -1,-1,-1,-1, 5, 6, 7, 4,-1,-1,-1,-1 },   // E
00620           {  4, 5, 6, 7,11, 8, 9,10,11, 8, 9,10 },   // SW
00621           {  0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11 },   // center
00622           {  1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 },   // NE
00623           { -1,-1,-1,-1, 7, 4, 5, 6,-1,-1,-1,-1 },   // W
00624           {  3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 },   // NW
00625           {  2, 3, 0, 1,-1,-1,-1,-1, 0, 1, 2, 3 } }; // N
00626   static const int swaparray[][12] =
00627         { {  0,0,0,0,0,0,0,0,3,3,3,3 },   // S
00628           {  0,0,0,0,0,0,0,0,6,6,6,6 },   // SE
00629           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // E
00630           {  0,0,0,0,0,0,0,0,5,5,5,5 },   // SW
00631           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // center
00632           {  5,5,5,5,0,0,0,0,0,0,0,0 },   // NE
00633           {  0,0,0,0,0,0,0,0,0,0,0,0 },   // W
00634           {  6,6,6,6,0,0,0,0,0,0,0,0 },   // NW
00635           {  3,3,3,3,0,0,0,0,0,0,0,0 } }; // N
00636 
00637   int ix, iy, face_num;
00638   (scheme_==RING) ?
00639     ring2xyf(pix,ix,iy,face_num) : nest2xyf(pix,ix,iy,face_num);
00640 
00641   const int nsm1 = nside_-1;
00642   if ((ix>0)&&(ix<nsm1)&&(iy>0)&&(iy<nsm1))
00643     {
00644     if (scheme_==RING)
00645       for (int m=0; m<8; ++m)
00646         result[m] = xyf2ring(ix+xoffset[m],iy+yoffset[m],face_num);
00647     else
00648       for (int m=0; m<8; ++m)
00649         result[m] = xyf2nest(ix+xoffset[m],iy+yoffset[m],face_num);
00650     }
00651   else
00652     {
00653     for (int i=0; i<8; ++i)
00654       {
00655       int x=ix+xoffset[i];
00656       int y=iy+yoffset[i];
00657       int nbnum=4;
00658       if (x<0)
00659         { x+=nside_; nbnum-=1; }
00660       else if (x>=nside_)
00661         { x-=nside_; nbnum+=1; }
00662       if (y<0)
00663         { y+=nside_; nbnum-=3; }
00664       else if (y>=nside_)
00665         { y-=nside_; nbnum+=3; }
00666 
00667       int f = facearray[nbnum][face_num];
00668       if (f>=0)
00669         {
00670         if (swaparray[nbnum][face_num]&1) x=nside_-x-1;
00671         if (swaparray[nbnum][face_num]&2) y=nside_-y-1;
00672         if (swaparray[nbnum][face_num]&4) std::swap(x,y);
00673         result[i] = (scheme_==RING) ? xyf2ring(x,y,f) : xyf2nest(x,y,f);
00674         }
00675       else
00676         result[i] = -1;
00677       }
00678     }
00679   }
00680 
00681 void Healpix_Base::get_ring_info2 (int ring, int &startpix, int &ringpix,
00682   double &theta, bool &shifted) const
00683   {
00684   int northring = (ring>2*nside_) ? 4*nside_-ring : ring;
00685   if (northring < nside_)
00686     {
00687     double tmp = northring*northring*fact2_;
00688     double costheta = 1 - tmp;
00689     double sintheta = sqrt(tmp*(2-tmp));
00690     theta = atan2(sintheta,costheta);
00691     ringpix = 4*northring;
00692     shifted = true;
00693     startpix = 2*northring*(northring-1);
00694     }
00695   else
00696     {
00697     theta = acos((2*nside_-northring)*fact1_);
00698     ringpix = 4*nside_;
00699     shifted = ((northring-nside_) & 1) == 0;
00700     startpix = ncap_ + (northring-nside_)*ringpix;
00701     }
00702   if (northring != ring) // southern hemisphere
00703     {
00704     theta = pi-theta;
00705     startpix = npix_ - startpix - ringpix;
00706     }
00707   }
00708 
00709 void Healpix_Base::get_interpol (const pointing &ptg, fix_arr<int,4> &pix,
00710   fix_arr<double,4> &wgt) const
00711   {
00712   double z = cos (ptg.theta);
00713   int ir1 = ring_above(z);
00714   int ir2 = ir1+1;
00715   double theta1, theta2, w1, tmp, dphi;
00716   int sp,nr;
00717   bool shift;
00718   int i1,i2;
00719   if (ir1>0)
00720     {
00721     get_ring_info2 (ir1, sp, nr, theta1, shift);
00722     dphi = twopi/nr;
00723     tmp = (ptg.phi/dphi - .5*shift);
00724     i1 = (tmp<0) ? int(tmp)-1 : int(tmp);
00725     w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
00726     i2 = i1+1;
00727     if (i1<0) i1 +=nr;
00728     if (i2>=nr) i2 -=nr;
00729     pix[0] = sp+i1; pix[1] = sp+i2;
00730     wgt[0] = 1-w1; wgt[1] = w1;
00731     }
00732   if (ir2<(4*nside_))
00733     {
00734     get_ring_info2 (ir2, sp, nr, theta2, shift);
00735     dphi = twopi/nr;
00736     tmp = (ptg.phi/dphi - .5*shift);
00737     i1 = (tmp<0) ? int(tmp)-1 : int(tmp);
00738     w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
00739     i2 = i1+1;
00740     if (i1<0) i1 +=nr;
00741     if (i2>=nr) i2 -=nr;
00742     pix[2] = sp+i1; pix[3] = sp+i2;
00743     wgt[2] = 1-w1; wgt[3] = w1;
00744     }
00745 
00746   if (ir1==0)
00747     {
00748     double wtheta = ptg.theta/theta2;
00749     wgt[2] *= wtheta; wgt[3] *= wtheta;
00750     double fac = (1-wtheta)*0.25;
00751     wgt[0] = fac; wgt[1] = fac; wgt[2] += fac; wgt[3] +=fac;
00752     pix[0] = (pix[2]+2)%4;
00753     pix[1] = (pix[3]+2)%4;
00754     }
00755   else if (ir2==4*nside_)
00756     {
00757     double wtheta = (ptg.theta-theta1)/(pi-theta1);
00758     wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
00759     double fac = wtheta*0.25;
00760     wgt[0] += fac; wgt[1] += fac; wgt[2] = fac; wgt[3] =fac;
00761     pix[2] = ((pix[0]+2)&3)+npix_-4;
00762     pix[3] = ((pix[1]+2)&3)+npix_-4;
00763     }
00764   else
00765     {
00766     double wtheta = (ptg.theta-theta1)/(theta2-theta1);
00767     wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
00768     wgt[2] *= wtheta; wgt[3] *= wtheta;
00769     }
00770 
00771   if (scheme_==NEST)
00772     for (int m=0; m<pix.size(); ++m)
00773       pix[m] = ring2nest(pix[m]);
00774   }
00775 
00776 double Healpix_Base::max_pixrad() const
00777   {
00778   vec3 va,vb;
00779   va.set_z_phi (2./3., pi/(4*nside_));
00780   double t1 = 1.-1./nside_;
00781   t1*=t1;
00782   vb.set_z_phi (1-t1/3, 0);
00783   return v_angle(va,vb);
00784   }