Changeset 1195 in Sophya for trunk/SophyaLib


Ignore:
Timestamp:
Sep 20, 2000, 5:30:03 PM (25 years ago)
Author:
ansari
Message:

fichier d'utilitaire HEALPix

Location:
trunk/SophyaLib/SkyMap
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/SkyMap/spherehealpix.cc

    r980 r1195  
    77#include "spherehealpix.h"
    88#include "strutil.h"
    9 
     9#include "HEALPixUtils.h"
    1010extern "C"
    1111{
    … …  
    1515}
    1616     
    17 
    18 //*******************************************************************
    19 // Class PIXELS_XY
    20 // Construction des tableaux necessaires a la traduction des indices RING en
    21 // indices NESTED (ou l'inverse)
    22 //*******************************************************************
    23 
    24 PIXELS_XY::PIXELS_XY()
    25 {
    26   pix2x_.ReSize(1024);
    27   pix2x_.Reset();
    28   pix2y_.ReSize(1024);
    29   pix2y_.Reset();
    30   x2pix_.ReSize(128);
    31   x2pix_.Reset();
    32   y2pix_.ReSize(128);
    33   y2pix_.Reset();
    34   mk_pix2xy();
    35   mk_xy2pix();
    36 }
    37 
    38 PIXELS_XY& PIXELS_XY::instance()
    39 {
    40   static PIXELS_XY single;
    41   return (single);
    42 }
    43 
    44 void PIXELS_XY::mk_pix2xy()
    45 {
    46   /*
    47     ==================================================
    48     subroutine mk_pix2xy
    49     ==================================================
    50     c     constructs the array giving x and y in the face from pixel number
    51     c     for the nested (quad-cube like) ordering of pixels
    52     c
    53     c     the bits corresponding to x and y are interleaved in the pixel number
    54     c     one breaks up the pixel number by even and odd bits
    55     ==================================================
    56   */
    57   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    58   //  (16/12/98)
    59 
    60   int kpix, jpix, IX, IY, IP, ID;
    61  
    62   for(kpix = 0; kpix < 1024; kpix++)
    63     {
    64       jpix = kpix;
    65       IX = 0;
    66       IY = 0;
    67       IP = 1 ;//  ! bit position (in x and y)
    68       while( jpix!=0 )
    69         { // ! go through all the bits
    70           ID=jpix%2;//  ! bit value (in kpix), goes in ix
    71           jpix = jpix/2;
    72           IX = ID*IP+IX;
    73      
    74           ID=jpix%2;//  ! bit value (in kpix), goes in iy
    75           jpix = jpix/2;
    76           IY = ID*IP+IY;
    77      
    78           IP = 2*IP;//        ! next bit (in x and y)
    79         }
    80       pix2x_(kpix) = IX;//     ! in 0,31
    81       pix2y_(kpix) = IY;//     ! in 0,31
    82     }
    83 }
    84 
    85 void PIXELS_XY::mk_xy2pix()
    86 {
    87   /*
    88     =================================================
    89     subroutine mk_xy2pix
    90     =================================================
    91     c     sets the array giving the number of the pixel lying in (x,y)
    92     c     x and y are in {1,128}
    93     c     the pixel number is in {0,128**2-1}
    94     c
    95     c     if  i-1 = sum_p=0  b_p * 2^p
    96     c     then ix = sum_p=0  b_p * 4^p
    97     c          iy = 2*ix
    98     c     ix + iy in {0, 128**2 -1}
    99     =================================================
    100   */
    101   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    102   //  (16/12/98)
    103 
    104   int K,IP,I,J,ID;
    105   for(I = 1; I <= 128; I++)
    106     {
    107       J  = I-1;//            !pixel numbers
    108       K  = 0;//
    109       IP = 1;//
    110       truc : if( J==0 )
    111         {
    112           x2pix_(I-1) = K;
    113           y2pix_(I-1) = 2*K;
    114         }
    115       else
    116         {
    117           ID = (int)fmod(J,2);
    118           J  = J/2;
    119           K  = IP*ID+K;
    120           IP = IP*4;
    121           goto truc;
    122         }
    123     }     
    124 }
     17using namespace SOPHYA;
    12518
    12619//*******************************************************************
    … …  
    659552  return  ring2nest(nSide_,k);
    660553}
    661 
    662 
    663 template<class T>
    664 int_4 SphereHEALPix<T>::nest2ring(int_4 nside, int_4 ipnest) const
    665 {
    666   /*
    667     ====================================================
    668     subroutine nest2ring(nside, ipnest, ipring)
    669     ====================================================
    670     c     conversion from NESTED to RING pixel number
    671     ====================================================
    672   */
    673   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    674   //  (16/12/98)
    675  
    676   const PIXELS_XY& PXY= PIXELS_XY::instance();
    677 
    678   int npix, npface, face_num, ncap, n_before;
    679   int ipf, ip_low, ip_trunc, ip_med, ip_hi;
    680   int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
    681   int ns_max=8192;
    682   int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
    683   int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
    684  
    685   if(  nside<1 ||  nside>ns_max ) {
    686     cout << "nside out of range" << endl;
    687     exit(0);
    688   }
    689   npix = 12 *  nside* nside;
    690   if( ipnest<0 || ipnest>npix-1 ) {
    691     cout << "ipnest out of range" << endl;
    692     exit(0);
    693   }
    694 
    695   ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
    696   nl4   = 4* nside;
    697  
    698   //c     finds the face, and the number in the face
    699   npface =  nside* nside;
    700   //cccccc      ip = ipnest - 1         ! in {0,npix-1}
    701  
    702   face_num = ipnest/npface;//  ! face number in {0,11}
    703   ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
    704   //c     finds the x,y on the face (starting from the lowest corner)
    705   //c     from the pixel number
    706   ip_low=ipf%1024;                //   ! content of the last 10 bits
    707   ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
    708   ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
    709   ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
    710  
    711   ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
    712   iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
    713  
    714   //c     transforms this in (horizontal, vertical) coordinates
    715   jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
    716   jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
    717  
    718   //c     computes the z coordinate on the sphere
    719   //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
    720   jr =  jrll[face_num]*nside - jrt - 1;
    721   nr = nside;//                  ! equatorial region (the most frequent)
    722   n_before = ncap + nl4 * (jr - nside);
    723   kshift=(jr - nside)%2;
    724   if( jr<nside ) {//then     ! north pole region
    725     nr = jr;
    726     n_before = 2 * nr * (nr - 1);
    727     kshift = 0;
    728   }
    729   else if( jr>3*nside ) {//then ! south pole region
    730     nr = nl4 - jr;
    731     n_before = npix - 2 * (nr + 1) * nr;
    732     kshift = 0;
    733   }
    734  
    735   //c     computes the phi coordinate on the sphere, in [0,2Pi]
    736   jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
    737  
    738   if( jp>nl4 ) jp = jp - nl4;
    739   if( jp<1 )   jp = jp + nl4;
    740  
    741   int aux=n_before + jp - 1;
    742   return (n_before + jp - 1);// ! in {0, npix-1}
    743 }
    744 
    745 template<class T>
    746 int_4 SphereHEALPix<T>::ring2nest(int_4 nside, int_4 ipring) const
    747 {
    748   /*
    749     ==================================================
    750     subroutine ring2nest(nside, ipring, ipnest)
    751     ==================================================
    752     c     conversion from RING to NESTED pixel number
    753     ==================================================
    754   */
    755   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    756   //  (16/12/98)
    757 
    758   const PIXELS_XY& PXY= PIXELS_XY::instance();
    759 
    760   double fihip, hip;
    761   int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
    762   int     kshift, face_num, nr;
    763   int irn, ire, irm, irs, irt, ifm , ifp;
    764   int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
    765   int ns_max(8192);
    766  
    767   // coordinate of the lowest corner of each face
    768   int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
    769   int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
    770  
    771   if( nside<1 || nside>ns_max ) {
    772     cout << "nside out of range" << endl;
    773     exit(0);
    774   }
    775   npix = 12 * nside*nside;
    776   if( ipring<0 || ipring>npix-1 ) {
    777     cout << "ipring out of range" << endl;
    778     exit(0);
    779   }
    780  
    781   nl2 = 2*nside;
    782   nl4 = 4*nside;
    783   npix = 12*nside*nside;//      ! total number of points
    784   ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
    785   ipring1 = ipring + 1;
    786  
    787   //c     finds the ring number, the position of the ring and the face number
    788   if( ipring1<=ncap ) {//then
    789    
    790     hip   = ipring1/2.;
    791     fihip = floor ( hip );
    792     irn   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
    793     iphi  = ipring1 - 2*irn*(irn - 1);
    794    
    795     kshift = 0;
    796     nr = irn   ;//               ! 1/4 of the number of points on the current ring
    797     face_num = (iphi-1) / irn;// ! in {0,3}
    798   }
    799   else if( ipring1<=nl2*(5*nside+1) ) {//then
    800    
    801     ip    = ipring1 - ncap - 1;
    802     irn   = (int)floor( ip / nl4 ) + nside;//               ! counted from North pole
    803     iphi  = (int)fmod(ip,nl4) + 1;
    804    
    805     kshift  = (int)fmod(irn+nside,2);//  ! 1 if irn+nside is odd, 0 otherwise
    806     nr = nside;
    807     ire =  irn - nside + 1;// ! in {1, 2*nside +1}
    808     irm =  nl2 + 2 - ire;
    809     ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
    810     ifp = (iphi - irm/2 + nside -1) / nside;
    811     if( ifp==ifm ) {//then          ! faces 4 to 7
    812       face_num = (int)fmod(ifp,4) + 4;
    813     }
    814     else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
    815       face_num = ifp;
    816     }
    817     else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
    818       face_num = ifp + 7;
    819     }
    820   }
    821   else {
    822    
    823     ip    = npix - ipring1 + 1;
    824     hip   = ip/2.;
    825     fihip = floor ( hip );
    826     irs   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//  ! counted from South pole
    827     iphi  = 4*irs + 1 - (ip - 2*irs*(irs-1));
    828    
    829     kshift = 0;
    830     nr = irs;
    831     irn   = nl4 - irs;
    832     face_num = (iphi-1) / irs + 8;// ! in {8,11}
    833   }
    834  
    835   //c     finds the (x,y) on the face
    836   irt =   irn  - jrll[face_num]*nside + 1;//       ! in {-nside+1,0}
    837   ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
    838 
    839 
    840   if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
    841  
    842   ix =  (ipt - irt ) / 2;
    843   iy = -(ipt + irt ) / 2;
    844  
    845   ix_low = (int)fmod(ix,128);
    846   ix_hi  = ix/128;
    847   iy_low = (int)fmod(iy,128);
    848   iy_hi  = iy/128;
    849   ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
    850 
    851   return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
    852 }
    853 
    854 template<class T>
    855 int_4 SphereHEALPix<T>::ang2pix_ring(int_4 nside, double theta, double phi) const
    856 {
    857   /*
    858     ==================================================
    859     c     gives the pixel number ipix (RING)
    860     c     corresponding to angles theta and phi
    861     c==================================================
    862   */
    863   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    864   //  (16/12/98)
    865 
    866   int nl2, nl4, ncap, npix, jp, jm, ipix1;
    867   double  z, za, tt, tp, tmp;
    868   int ir, ip, kshift;
    869 
    870   double piover2(Pi/2.);
    871   double twopi(2.*Pi);
    872   double z0(2./3.);
    873   int ns_max(8192);
    874 
    875   if( nside<1 || nside>ns_max ) {
    876     cout << "nside out of range" << endl;
    877     exit(0);
    878   }
    879 
    880   if( theta<0. || theta>Pi) {
    881     cout << "theta out of range" << endl;
    882     exit(0);
    883   }
    884 
    885   z = cos(theta);
    886   za = fabs(z);
    887   if( phi >= twopi)  phi = phi - twopi;
    888   if (phi < 0.)     phi = phi + twopi;
    889   tt = phi / piover2;//  ! in [0,4)
    890 
    891   nl2 = 2*nside;
    892   nl4 = 4*nside;
    893   ncap  = nl2*(nside-1);// ! number of pixels in the north polar cap
    894   npix  = 12*nside*nside;
    895 
    896   if( za <= z0 ) {
    897    
    898     jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of  ascending edge line
    899     jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
    900 
    901     ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
    902     kshift = 0;
    903     if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
    904    
    905     ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
    906     if( ip>nl4 ) ip = ip - nl4;
    907    
    908     ipix1 = ncap + nl4*(ir-1) + ip ;
    909   }
    910   else {
    911 
    912     tp = tt - floor(tt);//      !MOD(tt,1.d0)
    913     tmp = sqrt( 3.*(1. - za) );
    914 
    915     jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
    916     jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
    917 
    918     ir = jp + jm + 1;//        ! ring number counted from the closest pole
    919     ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
    920     if( ip>4*ir ) ip = ip - 4*ir;
    921      
    922     ipix1 = 2*ir*(ir-1) + ip;
    923       if( z<=0. ) {
    924         ipix1 = npix - 2*ir*(ir+1) + ip;
    925       }
    926   }
    927     return (ipix1 - 1);// ! in {0, npix-1}
    928 }
    929 
    930 template<class T>
    931 int_4 SphereHEALPix<T>::ang2pix_nest(int_4 nside, double theta, double phi) const
    932 {
    933   /*
    934     ==================================================
    935     subroutine ang2pix_nest(nside, theta, phi, ipix)
    936     ==================================================
    937     c     gives the pixel number ipix (NESTED)
    938     c     corresponding to angles theta and phi
    939     c
    940     c     the computation is made to the highest resolution available (nside=8192)
    941     c     and then degraded to that required (by integer division)
    942     c     this doesn't cost more, and it makes sure
    943     c     that the treatement of round-off will be consistent
    944     c     for every resolution
    945     ==================================================
    946   */
    947   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    948   //  (16/12/98)
    949  
    950   const PIXELS_XY& PXY= PIXELS_XY::instance();
    951 
    952   double    z, za, z0, tt, tp, tmp;
    953   int face_num,jp,jm;
    954   int ifp, ifm;
    955   int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
    956   double piover2(Pi/2.), twopi(2.*Pi);
    957   int ns_max(8192);
    958  
    959   if( nside<1 || nside>ns_max ) {
    960     cout << "nside out of range" << endl;
    961     exit(0);
    962   }
    963   if( theta<0 || theta>Pi ) {
    964     cout << "theta out of range" << endl;
    965     exit(0);
    966   }
    967   z  = cos(theta);
    968   za = fabs(z);
    969   z0 = 2./3.;
    970   if( phi>=twopi ) phi = phi - twopi;
    971   if( phi<0. )    phi = phi + twopi;
    972   tt = phi / piover2;// ! in [0,4[
    973   if( za<=z0 ) { // then ! equatorial region
    974    
    975     //(the index of edge lines increase when the longitude=phi goes up)
    976     jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
    977     jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
    978    
    979     //c        finds the face
    980     ifp = jp / ns_max;//  ! in {0,4}
    981     ifm = jm / ns_max;
    982     if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then  ! faces 4 to 7
    983     else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
    984     else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
    985    
    986     ix = (int)fmod(jm, ns_max);
    987     iy = ns_max - (int)fmod(jp, ns_max) - 1;
    988   }
    989   else { //! polar region, za > 2/3
    990    
    991     ntt = (int)floor(tt);
    992     if( ntt>=4 ) ntt = 3;
    993     tp = tt - ntt;
    994     tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
    995    
    996     //(the index of edge lines increase when distance from the closest pole goes up)
    997     jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
    998     jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
    999     jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
    1000     jm = (int)min(ns_max-1, jm);
    1001    
    1002     // finds the face and pixel's (x,y)
    1003     if( z>=0 ) {
    1004       face_num = ntt;//  ! in {0,3}
    1005       ix = ns_max - jm - 1;
    1006       iy = ns_max - jp - 1;
    1007     }
    1008     else {
    1009       face_num = ntt + 8;// ! in {8,11}
    1010       ix =  jp;
    1011       iy =  jm;
    1012     }
    1013   }
    1014  
    1015   ix_low = (int)fmod(ix,128);
    1016   ix_hi  =     ix/128;
    1017   iy_low = (int)fmod(iy,128);
    1018   iy_hi  =     iy/128;
    1019   ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
    1020   //  ipf = ipf / pow(ns_max/nside,2.);//  ! in {0, nside**2 - 1}
    1021   //  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
    1022   // $CHECK$  Reza 25/10/99 , pow remplace par *
    1023   ipf = ipf / ((ns_max/nside)*(ns_max/nside));
    1024   return (ipf + face_num*nside*nside);
    1025 }
    1026 
    1027 template<class T>
    1028 void SphereHEALPix<T>::pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const {
    1029   /*
    1030     ===================================================
    1031     c     gives theta and phi corresponding to pixel ipix (RING)
    1032     c     for a parameter nside
    1033     ===================================================
    1034   */
    1035   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    1036   //  (16/12/98)
    1037  
    1038   int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
    1039   double  fact1, fact2, fodd, hip, fihip;
    1040  
    1041   int ns_max(8192);
    1042  
    1043   if( nside<1 || nside>ns_max ) {
    1044     cout << "nside out of range" << endl;
    1045     exit(0);
    1046   }
    1047   npix = 12*nside*nside;      // ! total number of points
    1048   if( ipix<0 || ipix>npix-1 ) {
    1049     cout << "ipix out of range" << endl;
    1050     exit(0);
    1051   }
    1052  
    1053   ipix1 = ipix + 1; // in {1, npix}
    1054   nl2 = 2*nside;
    1055   nl4 = 4*nside;
    1056   ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
    1057   fact1 = 1.5*nside;
    1058   fact2 = 3.0*nside*nside;
    1059  
    1060   if( ipix1 <= ncap ) {  //! North Polar cap -------------
    1061    
    1062     hip   = ipix1/2.;
    1063     fihip = floor(hip);
    1064     iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
    1065     iphi  = ipix1 - 2*iring*(iring - 1);
    1066    
    1067     theta = acos( 1. - iring*iring / fact2 );
    1068     phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
    1069     //    cout << theta << " " << phi << endl;
    1070   }
    1071   else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
    1072    
    1073     ip    = ipix1 - ncap - 1;
    1074     iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
    1075     iphi  = ip%nl4 + 1;
    1076    
    1077     fodd  = 0.5 * (1 + (iring+nside)%2 );//  ! 1 if iring+nside is odd, 1/2 otherwise
    1078     theta = acos( (nl2 - iring) / fact1 );
    1079     phi   = ((double)iphi - fodd) * Pi /(2.*nside);
    1080   }
    1081   else {//! South Polar cap -----------------------------------
    1082    
    1083     ip    = npix - ipix1 + 1;
    1084     hip   = ip/2.;
    1085     fihip = floor(hip);
    1086     iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//     ! counted from South pole
    1087     iphi  = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
    1088    
    1089     theta = acos( -1. + iring*iring / fact2 );
    1090     phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
    1091     //    cout << theta << " " << phi << endl;
    1092   }
    1093 }
    1094 
    1095 template<class T>
    1096 void SphereHEALPix<T>::pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const {
    1097   /*
    1098     ==================================================
    1099     subroutine pix2ang_nest(nside, ipix, theta, phi)
    1100     ==================================================
    1101     c     gives theta and phi corresponding to pixel ipix (NESTED)
    1102     c     for a parameter nside
    1103     ==================================================
    1104   */
    1105   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    1106   //  (16/12/98)
    1107 
    1108   const PIXELS_XY& PXY= PIXELS_XY::instance();
    1109    
    1110   int npix, npface, face_num;
    1111   int ipf, ip_low, ip_trunc, ip_med, ip_hi;
    1112   int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
    1113   double z, fn, fact1, fact2;
    1114   double piover2(Pi/2.);
    1115   int ns_max(8192);
    1116          
    1117   // ! coordinate of the lowest corner of each face
    1118   int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
    1119   int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
    1120  
    1121   if( nside<1 || nside>ns_max ) {
    1122     cout << "nside out of range" << endl;
    1123     exit(0);
    1124   }
    1125   npix = 12 * nside*nside;
    1126   if( ipix<0 || ipix>npix-1 ) {
    1127     cout << "ipix out of range" << endl;
    1128     exit(0);
    1129   }
    1130  
    1131   fn = 1.*nside;
    1132   fact1 = 1./(3.*fn*fn);
    1133   fact2 = 2./(3.*fn);
    1134   nl4   = 4*nside;
    1135 
    1136   //c     finds the face, and the number in the face
    1137   npface = nside*nside;
    1138  
    1139   face_num = ipix/npface;//  ! face number in {0,11}
    1140   ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
    1141  
    1142   //c     finds the x,y on the face (starting from the lowest corner)
    1143   //c     from the pixel number
    1144   ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
    1145   ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
    1146   ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
    1147   ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
    1148  
    1149   ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
    1150   iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
    1151  
    1152   //c     transforms this in (horizontal, vertical) coordinates
    1153   jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
    1154   jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
    1155  
    1156   //c     computes the z coordinate on the sphere
    1157   //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
    1158   jr =  jrll[face_num]*nside - jrt - 1;
    1159   nr = nside;//                  ! equatorial region (the most frequent)
    1160   z  = (2*nside-jr)*fact2;
    1161   kshift = (int)fmod(jr - nside, 2);
    1162   if( jr<nside ) { //then     ! north pole region
    1163     nr = jr;
    1164     z = 1. - nr*nr*fact1;
    1165     kshift = 0;
    1166   }
    1167   else {
    1168     if( jr>3*nside ) {// then ! south pole region
    1169       nr = nl4 - jr;
    1170       z = - 1. + nr*nr*fact1;
    1171       kshift = 0;
    1172     }
    1173   }
    1174   theta = acos(z);
    1175  
    1176   //c     computes the phi coordinate on the sphere, in [0,2Pi]
    1177   //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
    1178   jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
    1179   if( jp>nl4 ) jp = jp - nl4;
    1180   if( jp<1 )   jp = jp + nl4;
    1181   phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
    1182 }
    1183 
    1184554
    1185555

  • trunk/SophyaLib/SkyMap/spherehealpix.h

    r1145 r1195  
    216216inline  SphereHEALPix<T>& operator = (const SphereHEALPix<T>& a)
    217217                                                        {return Set(a);}
    218        
     218
     219       
    219220private :
    220221
    … …  
    223224void SetThetaSlices();
    224225
    225 int_4  nest2ring(int_4 nside,int_4 ipnest) const;
    226 int_4  ring2nest(int_4 nside,int_4 ipring) const;
    227 
    228 int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
    229 int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
    230 void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
    231 void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
     226//int_4  nest2ring(int_4 nside,int_4 ipnest) const;
     227//int_4  ring2nest(int_4 nside,int_4 ipring) const;
     228
     229//int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
     230//int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
     231//void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
     232//void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
    232233inline void setParameters(int_4 nside, int_4 nbpixels, double solangle)
    233234  {
    … …  
    254255
    255256
    256 //////////////////////////////////////////////////////////////////////////
    257 //
    258 // ------------- Classe PIXELS_XY -----------------------
    259 //
    260 class PIXELS_XY
    261 {
    262 
    263 public :
    264 
    265 static PIXELS_XY& instance();
    266 
    267 NDataBlock<int_4> pix2x_;
    268 NDataBlock<int_4> pix2y_;
    269 NDataBlock<int_4> x2pix_;
    270 NDataBlock<int_4> y2pix_;
    271 
    272 private :
    273 
    274 PIXELS_XY();
    275 void mk_pix2xy();
    276 void mk_xy2pix();
    277 };
    278257
    279258} // Fin du namespace
Note: See TracChangeset for help on using the changeset viewer.