Index: unk/SophyaLib/SkyMap/spheregorski.cc
===================================================================
--- /trunk/SophyaLib/SkyMap/spheregorski.cc	(revision 842)
+++ 	(revision )
@@ -1,1329 +1,0 @@
-#include "machdefs.h"
-#include <math.h>
-#include <complex>
-
-#include "pexceptions.h"
-#include "fiondblock.h"
-#include "spheregorski.h"
-#include "strutil.h"
-
-extern "C" 
-{
-#include <stdio.h>
-#include <stdlib.h>
-#include <unistd.h>
-}
-      
-
-//*******************************************************************
-// Class PIXELS_XY 
-// Construction des tableaux necessaires a la traduction des indices RING en 
-// indices NESTED (ou l'inverse)
-//*******************************************************************
-
-PIXELS_XY::PIXELS_XY()
-{
-  pix2x_.ReSize(1024);
-  pix2x_.Reset();
-  pix2y_.ReSize(1024);
-  pix2y_.Reset();
-  x2pix_.ReSize(128);
-  x2pix_.Reset();
-  y2pix_.ReSize(128);
-  y2pix_.Reset();
-  mk_pix2xy(); 
-  mk_xy2pix(); 
-}
-
-PIXELS_XY& PIXELS_XY::instance()
-{
-  static PIXELS_XY single;
-  return (single);
-}
-
-void PIXELS_XY::mk_pix2xy() 
-{
-  /*
-    ==================================================
-    subroutine mk_pix2xy
-    ==================================================
-    c     constructs the array giving x and y in the face from pixel number
-    c     for the nested (quad-cube like) ordering of pixels
-    c
-    c     the bits corresponding to x and y are interleaved in the pixel number
-    c     one breaks up the pixel number by even and odd bits
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int kpix, jpix, IX, IY, IP, ID;
-  
-  for(kpix = 0; kpix < 1024; kpix++) 
-    {
-      jpix = kpix;
-      IX = 0;
-      IY = 0;
-      IP = 1 ;//  ! bit position (in x and y)
-      while( jpix!=0 )
-	{ // ! go through all the bits
-	  ID=jpix%2;//  ! bit value (in kpix), goes in ix
-	  jpix = jpix/2;
-	  IX = ID*IP+IX;
-      
-	  ID=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)
-	}
-      pix2x_(kpix) = IX;//     ! in 0,31
-      pix2y_(kpix) = IY;//     ! in 0,31
-    }
-}
-
-void PIXELS_XY::mk_xy2pix() 
-{
-  /*
-    =================================================
-    subroutine mk_xy2pix
-    =================================================
-    c     sets the array giving the number of the pixel lying in (x,y)
-    c     x and y are in {1,128}
-    c     the pixel number is in {0,128**2-1}
-    c
-    c     if  i-1 = sum_p=0  b_p * 2^p
-    c     then ix = sum_p=0  b_p * 4^p
-    c          iy = 2*ix
-    c     ix + iy in {0, 128**2 -1}
-    =================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int K,IP,I,J,ID;
-  for(I = 1; I <= 128; I++) 
-    {
-      J  = I-1;//            !pixel numbers
-      K  = 0;//
-      IP = 1;//
-      truc : if( J==0 ) 
-	{
-	  x2pix_(I-1) = K;
-	  y2pix_(I-1) = 2*K;
-	}
-      else 
-	{
-	  ID = (int)fmod(J,2);
-	  J  = J/2;
-	  K  = IP*ID+K;
-	  IP = IP*4;
-	  goto truc;
-	}
-    }     
-}
-
-//*******************************************************************
-//++ 
-// Class	SphereGorski
-//
-// include      spheregorski.h strutil.h
-//
-//    Pixelisation Gorski  
-//
-//
-//|    -----------------------------------------------------------------------
-//|     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
-//|    -----------------------------------------------------------------------
-//
-//    the sphere is split in 12 diamond-faces containing nside**2 pixels each
-//
-//    the numbering of the pixels (in the nested scheme) is similar to
-//    quad-cube
-//    In each face the first pixel is in the lowest corner of the diamond
-//
-//    the faces are                    (x,y) coordinate on each face
-//|        .   .   .   .   <--- North Pole
-//|       / \ / \ / \ / \                          ^        ^     
-//|      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      /      
-//|       \ / \ / \ / \ /                        y   \    /  x   
-//|      4 . 5 . 6 . 7 . 4 <--- equator               \  /                     
-//|       / \ / \ / \ / \                              \/        
-//|      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner  
-//|       \ / \ / \ / \ /                                  
-//|        .   .   .   .   <--- South Pole
-//|
-//    phi:0               2Pi                                 
-//
-//    in the ring scheme pixels are numbered along the parallels
-//    the first parallel is the one closest to the north pole and so on
-//    on each parallel, pixels are numbered starting from the one closest
-//    to phi = 0
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-//++
-//
-// Links	Parents
-//
-//    SphericalMap 
-//-- 
-
-/* --Methode-- */
-//++
-// Titre	Constructors
-//--
-//++
-
-template<class T>
-SphereGorski<T>::SphereGorski()
-
-//--
-{
-  InitNul();
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-SphereGorski<T>::SphereGorski(int_4 m)
-
-//    m is the "nside" of the Gorski algorithm
-//
-//    The total number of pixels will be Npix =  12*nside**2
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-{
-
-  if(m <= 0 || m > 8192) 
-    {
-      cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
-      throw RangeCheckError("SphereGorski<T>::SphereGorski() - Out of bound nside (< 8192)!");
-    }
-  // verifier que m est une puissance de deux 
-  int x= m;
-  while(x%2 == 0) x/=2;
-  if(x != 1) 
-    {  
-      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
-      throw ParmError("SphereGorski<T>::SphereGorski() - nside != 2^n !");
-    }
-  InitNul();
-  Pixelize(m);
-  SetThetaSlices();
-}
-//++
-template<class T>
-SphereGorski<T>::SphereGorski(const SphereGorski<T>& s, bool share)
-  : pixels_(s.pixels_, share)
-//    copy constructor
-//--
-{
-  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
-
-  nSide_= s.nSide_;
-  nPix_ = s.nPix_;
-  omeg_ = s.omeg_;
-}
-
-//++
-// Titre	Destructor
-//--
-//++
-template<class T>
-SphereGorski<T>::~SphereGorski()
-
-//--
-{
-}
-
-//++
-// Titre	Public Methods
-//--
-
-//++
-template<class T>
-void SphereGorski<T>::Resize(int_4 m)
-
-//    m is the "nside" of the Gorski algorithm
-//
-//    The total number of pixels will be Npix =  12*nside**2
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-{
-  if (m<=0 || m> 8192) {
-    cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
-    exit(1);
-  }
-  // verifier que m est une puissance de deux 
-  int x= m;
-  while (x%2==0) x/=2;
-  if(x != 1) 
-    {  
-      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
-      exit(1);
-    }
-  InitNul();
-  Pixelize(m);
-  SetThetaSlices();
-}
-
-template<class T>
-void  SphereGorski<T>::Pixelize( int_4 m) 
-
-//    prépare la pixelisation Gorski (m a la même signification 
-//    que pour le constructeur)
-//
-//    
-{
-  // On memorise les arguments d'appel
-  nSide_= m;  
-
-  // Nombre total de pixels sur la sphere entiere
-  nPix_= 12*nSide_*nSide_;
-
-  // pour le moment les tableaux qui suivent seront ranges dans l'ordre 
-  // de l'indexation GORSKY "RING"
-  // on pourra ulterieurement changer de strategie et tirer profit
-  // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra 
-  // de pourquoi c'est faire
-
-  // Creation et initialisation du vecteur des contenus des pixels 
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
-
-  // solid angle per pixel   
-  omeg_= 4.0*Pi/nPix_;
-}
-
-template<class T>
-void SphereGorski<T>::InitNul()
-//
-//    initialise à zéro les variables de classe 
-{
-  nSide_= 0;
-  nPix_ = 0;
-  omeg_ = 0.;
-//  pixels_.Reset();  -  Il ne faut pas mettre les pixels a zero si share !
-}
-
-/* --Methode-- */
-//++
-template<class T>
-int_4 SphereGorski<T>::NbPixels() const
-
-//    Retourne le nombre de pixels du découpage
-//--
-{  
-  return(nPix_);
-}
-
-//++
-template<class T>
-uint_4 SphereGorski<T>::NbThetaSlices() const 
-
-//    Return number of slices in theta direction on the  sphere
-//--
-{
-  uint_4 nbslices = uint_4(4*nSide_-1);
-  if (nSide_<=0) 
-    {
-      nbslices = 0;
-      throw PException(" sphere not pixelized, NbSlice=0 ");
-    }
-  return nbslices;
-}
-
-//++
-template<class T>
-void SphereGorski<T>::GetThetaSlice(int_4 index, r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 
-
-//    For a theta-slice with index 'index', return : 
-//
-//    the corresponding "theta" 
-//
-//    a vector containing the phi's of the pixels of the slice
-//
-//    a vector containing the corresponding values of pixels 
-//
-//--
-{
-
-  if (index<0 || index >= NbThetaSlices()) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::GetThetaSlice : Pixel index out of range" <<endl;
-      THROW(rangeCheckErr);
-    }
-
-
-  int_4 iring= sliceBeginIndex_(index);
-  int_4 lring  =  sliceLenght_(index);
-
-
-  phi.ReSize(lring);
-  value.ReSize(lring);
-  double TH= 0.;
-  double FI= 0.;
-  for(int_4 kk = 0; kk < lring;kk++) 
-    {
-      PixThetaPhi(kk+iring,TH,FI);
-      phi(kk)= FI;
-      value(kk)= PixVal(kk+iring);
-    }
-  theta= TH;
-}
-//++
-//++
-
-template<class T>
-void SphereGorski<T>::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const 
-
-//    For a theta-slice with index 'sliceIndex', return : 
-//
-//    the corresponding "theta" 
-//    the corresponding "phi" for first pixel of the slice 
-//
-//    a vector containing the indices of the pixels of the slice
-//    (equally distributed in phi)
-//
-//    a vector containing the corresponding values of pixels 
-//
-//--
-{
-
-  if (sliceIndex<0 || sliceIndex >= NbThetaSlices()) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::GetThetaSlice : Pixel index out of range" <<endl;
-      THROW(rangeCheckErr);
-    }
-
-  int_4 iring= sliceBeginIndex_(sliceIndex);
-  int_4 lring  =  sliceLenght_(sliceIndex);
-  pixelIndices.ReSize(lring);
-  value.ReSize(lring);
-  for(int_4 kk = 0; kk < lring;kk++) 
-    {
-      pixelIndices(kk)= kk+iring;
-      value(kk)= PixVal(kk+iring);
-    }
-  PixThetaPhi(iring, theta, phi0);
-}
-//++
-template<class T>
-void SphereGorski<T>::SetThetaSlices()  
-
-//--
-{
-  sliceBeginIndex_.ReSize(4*nSide_-1);
-  sliceLenght_.ReSize(4*nSide_-1);
-  int sliceIndex;
-  for (sliceIndex=0; sliceIndex<  nSide_-1; sliceIndex++)
-    {
-      sliceBeginIndex_(sliceIndex)  = 2*sliceIndex*(sliceIndex+1);
-      sliceLenght_(sliceIndex) = 4*(sliceIndex+1);
-    }
-  for (sliceIndex= nSide_-1; sliceIndex<  3*nSide_; sliceIndex++)
-    {
-      sliceBeginIndex_(sliceIndex)  = 2*nSide_*(2*sliceIndex-nSide_+1);
-      sliceLenght_(sliceIndex) = 4*nSide_;
-    }
-  for (sliceIndex= 3*nSide_; sliceIndex< 4*nSide_-1; sliceIndex++)
-    {
-      int_4 nc= 4*nSide_-1-sliceIndex;
-      sliceBeginIndex_(sliceIndex)  = nPix_-2*nc*(nc+1); 
-      sliceLenght_(sliceIndex) = 4*nc;
-    }
-}
-
-/* --Methode-- */
-//++
-template<class T>
-T& SphereGorski<T>::PixVal(int_4 k)
-
-//    Return value of  pixel with "RING" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::PIxVal : exception  a mettre en place, k= " << k <<endl;
-      THROW(rangeCheckErr);
-    }
-  return pixels_(k);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-T const& SphereGorski<T>::PixVal(int_4 k) const
-
-//    Return value of  pixel with "RING" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::PIxVal : exception  a mettre en place, k= " << k <<endl;
-      THROW(rangeCheckErr);
-    }
-  return *(pixels_.Data()+k);
-}
-
-//++
-template<class T>
-T& SphereGorski<T>::PixValNest(int_4 k)
-
-//    Return value of  pixel with "NESTED" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    { 
-      //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range")); 
-      cout<<" SphereGorski::PIxValNest : exceptions a mettre en place" <<endl;
-      THROW(rangeCheckErr);
-    }
-  return pixels_(nest2ring(nSide_,k));
-}
-//++
-
-template<class T>
-T const& SphereGorski<T>::PixValNest(int_4 k) const
-
-//    Return value of  pixel with "NESTED" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    { 
-      //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range")); 
-      cout<<" SphereGorski::PIxValNest : exceptions a mettre en place" <<endl;
-      THROW(rangeCheckErr);
-  }
-  int_4 pix= nest2ring(nSide_,k);
-  return *(pixels_.Data()+pix);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-bool SphereGorski<T>::ContainsSph(double /*theta*/, double /*phi*/) const
-//--
-{
-return(true);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-int_4 SphereGorski<T>::PixIndexSph(double theta,double phi) const
-
-//    Return "RING" index of the pixel corresponding to 
-//    direction (theta, phi).
-//--
-{
-  return ang2pix_ring(nSide_,theta,phi);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::PixIndexSphNest(double theta,double  phi) const
-
-//    Return "NESTED" index of the pixel corresponding to 
-//    direction (theta, phi).
-//--
-{
-  return ang2pix_nest(nSide_,theta,phi);
-}
-
-
-/* --Methode-- */
-//++
-template<class T>
-void SphereGorski<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
-
-//   Return (theta,phi) coordinates of middle of  pixel with "RING" index k
-//--
-{
-  pix2ang_ring(nSide_,k,theta,phi);
-}
-
-template <class T>
-T SphereGorski<T>::SetPixels(T v)
-{
-pixels_.Reset(v);
-return(v);
-}
-
-//++
-template<class T>
-double SphereGorski<T>::PixSolAngle(int_4 /*dummy*/) const
-//    Pixel Solid angle  (steradians)
-//    All the pixels have the same solid angle. The dummy argument is
-//    for compatibility with eventual pixelizations which would not 
-//    fulfil this requirement.
-//--
-{
-  return omeg_;
-}
-
-//++
-template<class T>
-void SphereGorski<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
-
-//   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
-//--
-{   
-  pix2ang_nest(nSide_,k,theta,phi);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::NestToRing(int_4 k) const
-
-//    translation from NESTED index  into RING index 
-//
-//--
-{
-  return  nest2ring(nSide_,k);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::RingToNest(int_4 k) const
-//
-//    translation from  RING index  into NESTED index 
-//
-//--
-{
-  return  ring2nest(nSide_,k);
-}
-
-
-template<class T>
-int_4 SphereGorski<T>::nest2ring(int_4 nside, int_4 ipnest) const 
-{
-  /*
-    ====================================================
-    subroutine nest2ring(nside, ipnest, ipring)
-    ====================================================
-    c     conversion from NESTED to RING pixel number
-    ====================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  int npix, npface, face_num, ncap, n_before;
-  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
-  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
-  int ns_max=8192;
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
-  
-  if(  nside<1 ||  nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 *  nside* nside;
-  if( ipnest<0 || ipnest>npix-1 ) {
-    cout << "ipnest out of range" << endl;
-    exit(0);
-  }
-
-  ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
-  nl4   = 4* nside;
-  
-  //c     finds the face, and the number in the face
-  npface =  nside* nside;
-  //cccccc      ip = ipnest - 1         ! in {0,npix-1}
-  
-  face_num = ipnest/npface;//  ! face number in {0,11}
-  ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
-  //c     finds the x,y on the face (starting from the lowest corner)
-  //c     from the pixel number
-  ip_low=ipf%1024;                //   ! content of the last 10 bits
-  ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
-  ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
-  ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
-  
-  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
-  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
-  
-  //c     transforms this in (horizontal, vertical) coordinates
-  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
-  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
-  
-  //c     computes the z coordinate on the sphere
-  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
-  jr =  jrll[face_num]*nside - jrt - 1;
-  nr = nside;//                  ! equatorial region (the most frequent)
-  n_before = ncap + nl4 * (jr - nside);
-  kshift=(jr - nside)%2;
-  if( jr<nside ) {//then     ! north pole region
-    nr = jr;
-    n_before = 2 * nr * (nr - 1);
-    kshift = 0;
-  }
-  else if( jr>3*nside ) {//then ! south pole region
-    nr = nl4 - jr;
-    n_before = npix - 2 * (nr + 1) * nr;
-    kshift = 0;
-  }
-  
-  //c     computes the phi coordinate on the sphere, in [0,2Pi]
-  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
-  
-  if( jp>nl4 ) jp = jp - nl4;
-  if( jp<1 )   jp = jp + nl4;
-  
-  int aux=n_before + jp - 1;
-  return (n_before + jp - 1);// ! in {0, npix-1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ring2nest(int_4 nside, int_4 ipring) const 
-{
-  /*
-    ==================================================
-    subroutine ring2nest(nside, ipring, ipnest)
-    ==================================================
-    c     conversion from RING to NESTED pixel number
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  double fihip, hip;
-  int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
-  int     kshift, face_num, nr;
-  int irn, ire, irm, irs, irt, ifm , ifp;
-  int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
-  int ns_max(8192);
-  
-  // coordinate of the lowest corner of each face
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 * nside*nside;
-  if( ipring<0 || ipring>npix-1 ) {
-    cout << "ipring out of range" << endl;
-    exit(0);
-  }
-  
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  npix = 12*nside*nside;//      ! total number of points
-  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
-  ipring1 = ipring + 1;
-  
-  //c     finds the ring number, the position of the ring and the face number
-  if( ipring1<=ncap ) {//then
-    
-    hip   = ipring1/2.;
-    fihip = floor ( hip );
-    irn   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
-    iphi  = ipring1 - 2*irn*(irn - 1);
-    
-    kshift = 0;
-    nr = irn   ;//               ! 1/4 of the number of points on the current ring
-    face_num = (iphi-1) / irn;// ! in {0,3}
-  }
-  else if( ipring1<=nl2*(5*nside+1) ) {//then
-    
-    ip    = ipring1 - ncap - 1;
-    irn   = (int)floor( ip / nl4 ) + nside;//               ! counted from North pole
-    iphi  = (int)fmod(ip,nl4) + 1;
-    
-    kshift  = (int)fmod(irn+nside,2);//  ! 1 if irn+nside is odd, 0 otherwise
-    nr = nside;
-    ire =  irn - nside + 1;// ! in {1, 2*nside +1}
-    irm =  nl2 + 2 - ire;
-    ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
-    ifp = (iphi - irm/2 + nside -1) / nside;
-    if( ifp==ifm ) {//then          ! faces 4 to 7
-      face_num = (int)fmod(ifp,4) + 4;
-    }
-    else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
-      face_num = ifp;
-    }
-    else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
-      face_num = ifp + 7;
-    }
-  }
-  else {
-    
-    ip    = npix - ipring1 + 1;
-    hip   = ip/2.;
-    fihip = floor ( hip );
-    irs   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//  ! counted from South pole
-    iphi  = 4*irs + 1 - (ip - 2*irs*(irs-1));
-    
-    kshift = 0;
-    nr = irs;
-    irn   = nl4 - irs;
-    face_num = (iphi-1) / irs + 8;// ! in {8,11}
-  }
-  
-  //c     finds the (x,y) on the face
-  irt =   irn  - jrll[face_num]*nside + 1;//       ! in {-nside+1,0}
-  ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
-
-
-  if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
-  
-  ix =  (ipt - irt ) / 2;
-  iy = -(ipt + irt ) / 2;
-  
-  ix_low = (int)fmod(ix,128);
-  ix_hi  = ix/128;
-  iy_low = (int)fmod(iy,128);
-  iy_hi  = iy/128;
-  ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
-
-  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ang2pix_ring(int_4 nside, double theta, double phi) const 
-{
-  /*
-    ==================================================
-    c     gives the pixel number ipix (RING) 
-    c     corresponding to angles theta and phi
-    c==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int nl2, nl4, ncap, npix, jp, jm, ipix1;
-  double  z, za, tt, tp, tmp;
-  int ir, ip, kshift;
-
-  double piover2(Pi/2.);
-  double twopi(2.*Pi);
-  double z0(2./3.);
-  int ns_max(8192);
-
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-
-  if( theta<0. || theta>Pi) {
-    cout << "theta out of range" << endl;
-    exit(0);
-  }
-
-  z = cos(theta);
-  za = fabs(z);
-  if( phi >= twopi)  phi = phi - twopi;
-  if (phi < 0.)     phi = phi + twopi;
-  tt = phi / piover2;//  ! in [0,4)
-
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  ncap  = nl2*(nside-1);// ! number of pixels in the north polar cap
-  npix  = 12*nside*nside;
-
-  if( za <= z0 ) {
-    
-    jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of  ascending edge line 
-    jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
-
-    ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
-    kshift = 0;
-    if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
-    
-    ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
-    if( ip>nl4 ) ip = ip - nl4;
-    
-    ipix1 = ncap + nl4*(ir-1) + ip ;
-  }
-  else {
-
-    tp = tt - floor(tt);//      !MOD(tt,1.d0)
-    tmp = sqrt( 3.*(1. - za) );
-
-    jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
-    jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
-
-    ir = jp + jm + 1;//        ! ring number counted from the closest pole
-    ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
-    if( ip>4*ir ) ip = ip - 4*ir;
-      
-    ipix1 = 2*ir*(ir-1) + ip;
-      if( z<=0. ) {
-	ipix1 = npix - 2*ir*(ir+1) + ip;
-      }
-  }
-    return (ipix1 - 1);// ! in {0, npix-1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ang2pix_nest(int_4 nside, double theta, double phi) const 
-{
-  /*
-    ==================================================
-    subroutine ang2pix_nest(nside, theta, phi, ipix)
-    ==================================================
-    c     gives the pixel number ipix (NESTED) 
-    c     corresponding to angles theta and phi
-    c
-    c     the computation is made to the highest resolution available (nside=8192)
-    c     and then degraded to that required (by integer division)
-    c     this doesn't cost more, and it makes sure 
-    c     that the treatement of round-off will be consistent 
-    c     for every resolution
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  double    z, za, z0, tt, tp, tmp;
-  int face_num,jp,jm;
-  int ifp, ifm;
-  int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
-  double piover2(Pi/2.), twopi(2.*Pi);
-  int ns_max(8192);
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  if( theta<0 || theta>Pi ) {
-    cout << "theta out of range" << endl;
-    exit(0);
-  }
-  z  = cos(theta);
-  za = fabs(z);
-  z0 = 2./3.;
-  if( phi>=twopi ) phi = phi - twopi;
-  if( phi<0. )    phi = phi + twopi;
-  tt = phi / piover2;// ! in [0,4[
-  if( za<=z0 ) { // then ! equatorial region
-    
-    //(the index of edge lines increase when the longitude=phi goes up)
-    jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
-    jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
-    
-    //c        finds the face
-    ifp = jp / ns_max;//  ! in {0,4}
-    ifm = jm / ns_max;
-    if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then  ! faces 4 to 7
-    else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
-    else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
-    
-    ix = (int)fmod(jm, ns_max);
-    iy = ns_max - (int)fmod(jp, ns_max) - 1;
-  }
-  else { //! polar region, za > 2/3
-    
-    ntt = (int)floor(tt);
-    if( ntt>=4 ) ntt = 3;
-    tp = tt - ntt;
-    tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
-    
-    //(the index of edge lines increase when distance from the closest pole goes up)
-    jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
-    jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
-    jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
-    jm = (int)min(ns_max-1, jm);
-    
-    // finds the face and pixel's (x,y)
-    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;
-    }
-  }
-  
-  ix_low = (int)fmod(ix,128);
-  ix_hi  =     ix/128;
-  iy_low = (int)fmod(iy,128);
-  iy_hi  =     iy/128;
-  ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
-  //  ipf = ipf / pow(ns_max/nside,2.);//  ! in {0, nside**2 - 1}
-  //  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
-  // $CHECK$  Reza 25/10/99 , pow remplace par *
-  ipf = ipf / ((ns_max/nside)*(ns_max/nside));
-  return (ipf + face_num*nside*nside);
-}
-
-template<class T>
-void SphereGorski<T>::pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const {
-  /*
-    ===================================================
-    c     gives theta and phi corresponding to pixel ipix (RING) 
-    c     for a parameter nside
-    ===================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
-  double  fact1, fact2, fodd, hip, fihip;
-  
-  int ns_max(8192);
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12*nside*nside;      // ! total number of points
-  if( ipix<0 || ipix>npix-1 ) {
-    cout << "ipix out of range" << endl;
-    exit(0);
-  }
-  
-  ipix1 = ipix + 1; // in {1, npix}
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
-  fact1 = 1.5*nside;
-  fact2 = 3.0*nside*nside;
-  
-  if( ipix1 <= ncap ) {  //! North Polar cap -------------
-    
-    hip   = ipix1/2.;
-    fihip = floor(hip);
-    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
-    iphi  = ipix1 - 2*iring*(iring - 1);
-    
-    theta = acos( 1. - iring*iring / fact2 );
-    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
-    //    cout << theta << " " << phi << endl;
-  }
-  else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
-    
-    ip    = ipix1 - ncap - 1;
-    iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
-    iphi  = ip%nl4 + 1;
-    
-    fodd  = 0.5 * (1 + (iring+nside)%2 );//  ! 1 if iring+nside is odd, 1/2 otherwise
-    theta = acos( (nl2 - iring) / fact1 );
-    phi   = ((double)iphi - fodd) * Pi /(2.*nside);
-  }
-  else {//! South Polar cap -----------------------------------
-    
-    ip    = npix - ipix1 + 1;
-    hip   = ip/2.;
-    fihip = floor(hip);
-    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//     ! counted from South pole
-    iphi  = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
-    
-    theta = acos( -1. + iring*iring / fact2 );
-    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
-    //    cout << theta << " " << phi << endl;
-  }
-}
-
-template<class T>
-void SphereGorski<T>::pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const {
-  /*
-    ==================================================
-    subroutine pix2ang_nest(nside, ipix, theta, phi)
-    ==================================================
-    c     gives theta and phi corresponding to pixel ipix (NESTED) 
-    c     for a parameter nside
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-    
-  int npix, npface, face_num;
-  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
-  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
-  double z, fn, fact1, fact2;
-  double piover2(Pi/2.);
-  int ns_max(8192);
-          
-  // ! coordinate of the lowest corner of each face
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 * nside*nside;
-  if( ipix<0 || ipix>npix-1 ) {
-    cout << "ipix out of range" << endl;
-    exit(0);
-  }
-  
-  fn = 1.*nside;
-  fact1 = 1./(3.*fn*fn);
-  fact2 = 2./(3.*fn);
-  nl4   = 4*nside;
-
-  //c     finds the face, and the number in the face
-  npface = nside*nside;
-  
-  face_num = ipix/npface;//  ! face number in {0,11}
-  ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
-  
-  //c     finds the x,y on the face (starting from the lowest corner)
-  //c     from the pixel number
-  ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
-  ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
-  ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
-  ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
-  
-  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
-  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
-  
-  //c     transforms this in (horizontal, vertical) coordinates
-  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
-  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
-  
-  //c     computes the z coordinate on the sphere
-  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
-  jr =  jrll[face_num]*nside - jrt - 1;
-  nr = nside;//                  ! equatorial region (the most frequent)
-  z  = (2*nside-jr)*fact2;
-  kshift = (int)fmod(jr - nside, 2);
-  if( jr<nside ) { //then     ! north pole region
-    nr = jr;
-    z = 1. - nr*nr*fact1;
-    kshift = 0;
-  }
-  else {
-    if( jr>3*nside ) {// then ! south pole region
-      nr = nl4 - jr;
-      z = - 1. + nr*nr*fact1;
-      kshift = 0;
-    }
-  }
-  theta = acos(z);
-  
-  //c     computes the phi coordinate on the sphere, in [0,2Pi]
-  //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
-  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
-  if( jp>nl4 ) jp = jp - nl4;
-  if( jp<1 )   jp = jp + nl4;
-  phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
-}
-
-
-
-template <class T>
-void SphereGorski<T>::print(ostream& os) const
-{
-  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
-  //
-  os << " nSide_ = " << nSide_  << endl;
-  os << " nPix_   = " << nPix_   << endl;
-  os << " omeg_ = " << omeg_   << endl;
-
-  os << " content of pixels : ";
-  for(int i=0; i < nPix_; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os <<  pixels_(i) <<", ";
-    }
-  os << endl;
-
-  os << endl;
-  //const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  //os << endl;  os << " contenu des tableaux conversions "<<endl;
-  //for(int i=0; i < 5; i++)
-  //  {
-  //   os<<PXY.pix2x_(i)<<", "<<PXY.pix2y_(i)<<", "<<PXY.x2pix_(i)<<", "<<PXY.y2pix_(i)<<endl;
-  // }
-  os << endl;
-
-}
-
-//*******************************************************************
-// Class FIO_SphereGorski<T>
-//  Les objets delegues pour la gestion de persistance 
-//*******************************************************************
-
-template <class T>
-FIO_SphereGorski<T>::FIO_SphereGorski()
-{
-  dobj= new SphereGorski<T>;
-  ownobj= true;
-}
-
-template <class T>
-FIO_SphereGorski<T>::FIO_SphereGorski(string const& filename)
-{
-  dobj= new SphereGorski<T>;
-  dobj->pixels_.SetTemp(true);
-  ownobj= true;
-  Read(filename);
-}
-
-template <class T>
-FIO_SphereGorski<T>::FIO_SphereGorski(const SphereGorski<T>& obj)
-{
-  dobj= new SphereGorski<T>(obj, true);
-  dobj->pixels_.SetTemp(true);
-  ownobj= true;
-}
-
-template <class T>
-FIO_SphereGorski<T>::FIO_SphereGorski(SphereGorski<T>* obj)
-{
-  dobj= obj;
-  ownobj= false;
-}
-
-template <class T>
-FIO_SphereGorski<T>::~FIO_SphereGorski()
-{
-  if (ownobj && dobj) delete dobj;
-}
-
-template <class T>
-AnyDataObj* FIO_SphereGorski<T>::DataObj()
-{
-  return(dobj);
-}
-
-
-template <class T>
-void FIO_SphereGorski<T>::SetDataObj(AnyDataObj & o)
-{
-  SphereGorski<T> * po = dynamic_cast< SphereGorski<T> * >(&o);
-  if (po == NULL) return;
-  if (ownobj && dobj) delete dobj;
-  dobj = po; ownobj = false;
-} 
-
-template <class T>
-void FIO_SphereGorski<T>::ReadSelf(PInPersist& is)
-{
-
-  if(dobj == NULL) 
-    {
-      dobj= new SphereGorski<T>;
-      dobj->pixels_.SetTemp(true);
-      ownobj= true;      
-    }
-
-// On lit les 3 premiers uint_8
-uint_8 itab[3];
-is.Get(itab, 3);
-// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Read(is);
-  // Pour savoir s'il y avait un DVList Info associe
-  char strg[256];
-  is.GetLine(strg, 255);
-  bool hadinfo= false;
-  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
-  if(hadinfo) 
-    {    // Lecture eventuelle du DVList Info
-      is >> dobj->Info(); 
-    }
-
-  int_4 nSide;
-  is.GetI4(nSide);
-  int_4 nPix;
-  is.GetI4(nPix);
-  double Omega;
-  is.GetR8(Omega);
-  dobj->setParameters(nSide,nPix, Omega);
-
-// On lit les DataBlocks;
-  FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-  fio_nd.Read(is);
-  FIO_NDataBlock<int> fio_ndi1(&dobj->sliceBeginIndex_);
-  fio_ndi1.Read(is);
-  FIO_NDataBlock<int> fio_ndi2(&dobj->sliceLenght_);
-  fio_ndi2.Read(is);
-}
-
-template <class T>
-void FIO_SphereGorski<T>::WriteSelf(POutPersist& os) const
-{
-  if(dobj == NULL) 
-    {
-      cout << " WriteSelf:: dobj= null " << endl;
-      return;
-    }
-
-//  On ecrit 3 uint_8 
-//  0 : Numero de version,  1 : Size index,  2  reserve a l
-uint_8 itab[3];
-itab[0] = 1;
-itab[1] = dobj->SizeIndex();
-itab[2] = 0;
-os.Put(itab, 3);
-// Let's write the SphereCoordSys object
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Write(os);
-
-  char strg[256];
-  int_4 nSide= dobj->SizeIndex();
-  int_4 nPix = dobj->NbPixels();
- 
-  if(dobj->ptrInfo()) 
-    {
-      sprintf(strg,"SphereGorski: NSide=%6d  NPix=%9d HasInfo",nSide,nPix);
-      os.PutLine(strg);
-      os << dobj->Info();
-    }
-  else 
-    { 
-      sprintf(strg,"SphereGorski: NSide=%6d  NPix=%9d ",nSide,nPix);
-      os.PutLine(strg);  
-    }
-
-  os.PutI4(nSide);
-  os.PutI4(nPix);
-  os.PutR8(dobj->PixSolAngle());
-// On ecrit les dataBlocks 
-  FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-  fio_nd.Write(os);
-  FIO_NDataBlock<int> fio_ndi1(&dobj->sliceBeginIndex_);
-  fio_ndi1.Write(os);
-  FIO_NDataBlock<int> fio_ndi2(&dobj->sliceLenght_);
-  fio_ndi2.Write(os);
-}
-
-#ifdef __CXX_PRAGMA_TEMPLATES__
-#pragma define_template SphereGorski<uint_2>
-#pragma define_template SphereGorski<r_8>
-#pragma define_template SphereGorski<r_4>
-#pragma define_template SphereGorski< complex<r_4> >
-#pragma define_template SphereGorski< complex<r_8> >
-#pragma define_template FIO_SphereGorski<r_8>
-#pragma define_template FIO_SphereGorski<r_4>
-#pragma define_template FIO_SphereGorski< complex<r_4> >
-#pragma define_template FIO_SphereGorski< complex<r_8> >
-#endif
-#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
-template class SphereGorski<uint_2>;
-template class SphereGorski<r_8>;
-template class SphereGorski<r_4>;
-template class SphereGorski< complex<r_4> >;
-template class SphereGorski< complex<r_8> >;
-template class FIO_SphereGorski<r_8>;
-template class FIO_SphereGorski<r_4>;
-template class FIO_SphereGorski< complex<r_4> >;
-template class FIO_SphereGorski< complex<r_8> >;
-#endif
-
Index: unk/SophyaLib/SkyMap/spheregorski.h
===================================================================
--- /trunk/SophyaLib/SkyMap/spheregorski.h	(revision 842)
+++ 	(revision )
@@ -1,275 +1,0 @@
-#ifndef SPHEREGORSKI_SEEN
-#define SPHEREGORSKI_SEEN
-
-#include "sphericalmap.h"
-#include "tvector.h"
-#include "ndatablock.h"
-
-#include "anydataobj.h"
-#include "ppersist.h"
-
-
-// ***************** CLASSE SphereGorski *****************************
-
-//!  class SphereGorski
-/*!
-   Pixelisation Gorski  
-
-
-    -----------------------------------------------------------------------
-     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
-    -----------------------------------------------------------------------
-
-    the sphere is split in 12 diamond-faces containing nside**2 pixels each
-
-    the numbering of the pixels (in the nested scheme) is similar to
-    quad-cube
-    In each face the first pixel is in the lowest corner of the diamond
-
-    the faces are                    (x,y) coordinate on each face
-\verbatim
-        .   .   .   .   <--- North Pole
-       / \ / \ / \ / \                          ^        ^     
-      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      / 
-       \ / \ / \ / \ /                        y   \    /  x  
-      4 . 5 . 6 . 7 . 4 <--- equator               \  /      
-       / \ / \ / \ / \                              \/      
-      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner 
-       \ / \ / \ / \ /      
-        .   .   .   .   <--- South Pole
-\endverbatim
-    phi:0               2Pi                                 
-
-   in the ring scheme pixels are numbered along the parallels
-   the first parallel is the one closest to the north pole and so on  
-   on each parallel, pixels are numbered starting from the one closest
-    to phi = 0
-
-    nside MUST be a power of 2 (<= 8192)
-
-*/
-
-template<class T>
-class FIO_SphereGorski;
-
-template<class T>
-class SphereGorski : public SphericalMap<T>
-{
-
-  friend class FIO_SphereGorski<T>;
-
-public :
-
-SphereGorski();
-/*!    
-  m is the "nside" of the Gorski algorithm
-
-  The total number of pixels will be Npix =  12*nside**2
-
-  nside MUST be a power of 2 (<= 8192)
-*/
-SphereGorski(int_4 m);
-SphereGorski(const SphereGorski<T>& s, bool share=false);
-//!	Destructor
-virtual ~SphereGorski();
-
-/*! Setting blockdata to temporary (see ndatablock documentation) */
-inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};
-// ------------------ Definition of PixelMap abstract methods
-
-/* Nombre de pixels du decoupage */
-/*!    Return number of  pixels of the  splitting */
-virtual int_4 NbPixels() const;
-
-/* Valeur du contenu du pixel d'indice "RING" k  */
-/*!    Return value of  pixel with "RING" index k */
-virtual T& PixVal(int_4 k);
-virtual T const& PixVal(int_4 k) const;
-
-/* Nombre de tranches en theta */ 
-/*!    Return number of slices in theta direction on the  sphere */
-uint_4 NbThetaSlices() const;
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-    a vector containing the phi's of the pixels of the slice
-
-    a vector containing the corresponding values of pixels 
-*/
-void GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const;
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-   the corresponding "phi" for first pixel of the slice 
-
-    a vector containing indices of the pixels of the slice
-
-   (equally distributed in phi)
-
-    a vector containing the corresponding values of pixels 
-*/
- void GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const ;
-
-/* Return true if teta,phi in map  */
-virtual bool ContainsSph(double theta, double phi) const;
-/* Indice "RING" du pixel vers lequel pointe une direction definie par 
-ses  coordonnees spheriques */                                    
-/*!  Return "RING" index of the pixel corresponding to direction (theta, phi).
- */
-virtual int_4 PixIndexSph(double theta,double phi) const;
-
-/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
-virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;
-
-/*! Set all pixels to value v */
-virtual T SetPixels(T v);
-
-/* Pixel Solid angle  (steradians) */
-/*!    Pixel Solid angle  (steradians)
-
-    All the pixels have the same solid angle. The dummy argument is
-    for compatibility with eventual pixelizations which would not 
-   fulfil this requirement.
-*/
-virtual double PixSolAngle(int_4 dummy=0) const;
-
-/*  Acces to the DataBlock  */
-inline       NDataBlock<T>& DataBlock()       {return pixels_;}
-inline const NDataBlock<T>& DataBlock() const {return pixels_;}
-
-// --------------- Specific methods 
-
-/*!    
-  m is the "nside" of the Gorski algorithm
-
-  The total number of pixels will be Npix =  12*nside**2
-
-  nside MUST be a power of 2 (<= 8192)
-*/
-virtual void Resize(int_4 m);
-
-// pour l'instant le tableau est ordonne selon RING, uniquement
-inline virtual char* TypeOfMap() const {return "RING";};
-
-
-/* Valeur du contenu du pixel d'indice "NEST" k */
-/*!    Return value of  pixel with "NESTED" index k */
-virtual T& PixValNest(int_4 k);
-/*!    Return value of  pixel with "NESTED" index k */
-virtual T const& PixValNest(int_4 k) const;
-
-/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
-ses  coordonnees spheriques */                                    
-/*! Return "NESTED" index of the pixel corresponding to direction (theta, phi).
- */
-virtual int_4 PixIndexSphNest(double theta,double phi) const;
-
-/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
-/*!   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
- */
-virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;
-
-/* algorithme de pixelisation */
-void Pixelize(int_4); 
-
-/* convertit index nested en ring  */
-/*!    translation from NESTED index  into RING index */
-int_4 NestToRing(int_4) const;
-
-/* convertit index ring en nested" */
-/*!    translation from  RING index  into NESTED index */
-int_4 RingToNest(int_4) const;
-
-
-/* retourne la valeur du parametre Gorski */
-inline virtual int_4 SizeIndex() const {return(nSide_);}
-
-/* retourne les pointeurs /remplit les tableaux */
-
-
-
-
-/* impression */ 
-void print(ostream& os) const;
-         
-private :
-
-// ------------- méthodes internes ----------------------
-void InitNul();
- void SetThetaSlices(); 
-
-int_4  nest2ring(int_4 nside,int_4 ipnest) const;
-int_4  ring2nest(int_4 nside,int_4 ipring) const;
-
-int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
-int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
-void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
-void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
-inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 
-  {
-    nSide_= nside;
-    nPix_= nbpixels;
-    omeg_= solangle;
-  }
-
-// ------------- variables internes -----------------------
-
-int_4 nSide_;
-int_4 nPix_;
-double omeg_;
-
-NDataBlock<T> pixels_;
-NDataBlock<int> sliceBeginIndex_;
- NDataBlock<int> sliceLenght_;
-
-};
-
-//
-// ------------- Classe pour la gestion de persistance --
-//
-template <class T>
-class FIO_SphereGorski : public PPersist  
-{
-public:
-
-FIO_SphereGorski();
-FIO_SphereGorski(string const & filename); 
-FIO_SphereGorski(const SphereGorski<T>& obj);
-FIO_SphereGorski(SphereGorski<T>* obj);
-virtual ~FIO_SphereGorski();
-virtual AnyDataObj* DataObj();
-virtual void        SetDataObj(AnyDataObj & o);
-inline operator SphereGorski<T>() { return(*dobj); }
-
-protected :
-
-virtual void ReadSelf(PInPersist&);           
-virtual void WriteSelf(POutPersist&) const;  
-SphereGorski<T>* dobj;
-bool ownobj;
-};
-
-//
-// ------------- Classe PIXELS_XY -----------------------
-//
-class PIXELS_XY
-{
-
-public :
-
-static PIXELS_XY& instance();
-
-NDataBlock<int_4> pix2x_;
-NDataBlock<int_4> pix2y_;
-NDataBlock<int_4> x2pix_;
-NDataBlock<int_4> y2pix_;
-
-private :
-
-PIXELS_XY();
-void mk_pix2xy();
-void mk_xy2pix();
-};
-#endif
Index: /trunk/SophyaLib/SkyMap/spherehealpix.cc
===================================================================
--- /trunk/SophyaLib/SkyMap/spherehealpix.cc	(revision 843)
+++ /trunk/SophyaLib/SkyMap/spherehealpix.cc	(revision 843)
@@ -0,0 +1,1167 @@
+#include "machdefs.h"
+#include <math.h>
+#include <complex>
+
+#include "pexceptions.h"
+#include "fiondblock.h"
+#include "spherehealpix.h"
+#include "strutil.h"
+
+extern "C" 
+{
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+}
+      
+
+//*******************************************************************
+// Class PIXELS_XY 
+// Construction des tableaux necessaires a la traduction des indices RING en 
+// indices NESTED (ou l'inverse)
+//*******************************************************************
+
+PIXELS_XY::PIXELS_XY()
+{
+  pix2x_.ReSize(1024);
+  pix2x_.Reset();
+  pix2y_.ReSize(1024);
+  pix2y_.Reset();
+  x2pix_.ReSize(128);
+  x2pix_.Reset();
+  y2pix_.ReSize(128);
+  y2pix_.Reset();
+  mk_pix2xy(); 
+  mk_xy2pix(); 
+}
+
+PIXELS_XY& PIXELS_XY::instance()
+{
+  static PIXELS_XY single;
+  return (single);
+}
+
+void PIXELS_XY::mk_pix2xy() 
+{
+  /*
+    ==================================================
+    subroutine mk_pix2xy
+    ==================================================
+    c     constructs the array giving x and y in the face from pixel number
+    c     for the nested (quad-cube like) ordering of pixels
+    c
+    c     the bits corresponding to x and y are interleaved in the pixel number
+    c     one breaks up the pixel number by even and odd bits
+    ==================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  int kpix, jpix, IX, IY, IP, ID;
+  
+  for(kpix = 0; kpix < 1024; kpix++) 
+    {
+      jpix = kpix;
+      IX = 0;
+      IY = 0;
+      IP = 1 ;//  ! bit position (in x and y)
+      while( jpix!=0 )
+	{ // ! go through all the bits
+	  ID=jpix%2;//  ! bit value (in kpix), goes in ix
+	  jpix = jpix/2;
+	  IX = ID*IP+IX;
+      
+	  ID=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)
+	}
+      pix2x_(kpix) = IX;//     ! in 0,31
+      pix2y_(kpix) = IY;//     ! in 0,31
+    }
+}
+
+void PIXELS_XY::mk_xy2pix() 
+{
+  /*
+    =================================================
+    subroutine mk_xy2pix
+    =================================================
+    c     sets the array giving the number of the pixel lying in (x,y)
+    c     x and y are in {1,128}
+    c     the pixel number is in {0,128**2-1}
+    c
+    c     if  i-1 = sum_p=0  b_p * 2^p
+    c     then ix = sum_p=0  b_p * 4^p
+    c          iy = 2*ix
+    c     ix + iy in {0, 128**2 -1}
+    =================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  int K,IP,I,J,ID;
+  for(I = 1; I <= 128; I++) 
+    {
+      J  = I-1;//            !pixel numbers
+      K  = 0;//
+      IP = 1;//
+      truc : if( J==0 ) 
+	{
+	  x2pix_(I-1) = K;
+	  y2pix_(I-1) = 2*K;
+	}
+      else 
+	{
+	  ID = (int)fmod(J,2);
+	  J  = J/2;
+	  K  = IP*ID+K;
+	  IP = IP*4;
+	  goto truc;
+	}
+    }     
+}
+
+//*******************************************************************
+//++ 
+// Class	SphereHealpix
+//
+// include      SphereHealpix.h strutil.h
+//
+//    Pixelisation Gorski  
+//
+//
+//|    -----------------------------------------------------------------------
+//|     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
+//|    -----------------------------------------------------------------------
+//
+//    the sphere is split in 12 diamond-faces containing nside**2 pixels each
+//
+//    the numbering of the pixels (in the nested scheme) is similar to
+//    quad-cube
+//    In each face the first pixel is in the lowest corner of the diamond
+//
+//    the faces are                    (x,y) coordinate on each face
+//|        .   .   .   .   <--- North Pole
+//|       / \ / \ / \ / \                          ^        ^     
+//|      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      /      
+//|       \ / \ / \ / \ /                        y   \    /  x   
+//|      4 . 5 . 6 . 7 . 4 <--- equator               \  /                     
+//|       / \ / \ / \ / \                              \/        
+//|      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner  
+//|       \ / \ / \ / \ /                                  
+//|        .   .   .   .   <--- South Pole
+//|
+//    phi:0               2Pi                                 
+//
+//    in the ring scheme pixels are numbered along the parallels
+//    the first parallel is the one closest to the north pole and so on
+//    on each parallel, pixels are numbered starting from the one closest
+//    to phi = 0
+//
+//    nside MUST be a power of 2 (<= 8192)
+//--
+//++
+//
+// Links	Parents
+//
+//    SphericalMap 
+//-- 
+
+/* --Methode-- */
+//++
+// Titre	Constructors
+//--
+//++
+
+template<class T>
+SphereHealpix<T>::SphereHealpix()
+
+//--
+{
+  InitNul();
+  pixels_.Reset();
+}
+
+//++
+template<class T>
+SphereHealpix<T>::SphereHealpix(int_4 m)
+
+//    m is the "nside" of the Gorski algorithm
+//
+//    The total number of pixels will be Npix =  12*nside**2
+//
+//    nside MUST be a power of 2 (<= 8192)
+//--
+{
+
+  if(m <= 0 || m > 8192) 
+    {
+      cout << "SphereHealpix : m hors bornes [0,8192], m= " << m << endl;
+      throw RangeCheckError("SphereHealpix<T>::SphereHealpix() - Out of bound nside (< 8192)!");
+    }
+  // verifier que m est une puissance de deux 
+  int x= m;
+  while(x%2 == 0) x/=2;
+  if(x != 1) 
+    {  
+      cout<<"SphereHealpix: m doit etre une puissance de deux, m= "<<m<<endl;
+      throw ParmError("SphereHealpix<T>::SphereHealpix() - nside != 2^n !");
+    }
+  InitNul();
+  Pixelize(m);
+  SetThetaSlices();
+}
+//++
+template<class T>
+SphereHealpix<T>::SphereHealpix(const SphereHealpix<T>& s, bool share)
+  : pixels_(s.pixels_, share), sliceBeginIndex_(s.sliceBeginIndex_, share), 
+				sliceLenght_(s.sliceLenght_, share)
+//    copy constructor
+//--
+{
+  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
+
+  nSide_= s.nSide_;
+  nPix_ = s.nPix_;
+  omeg_ = s.omeg_;
+}
+
+//++
+// Titre	Destructor
+//--
+//++
+template<class T>
+SphereHealpix<T>::~SphereHealpix()
+
+//--
+{
+}
+
+//++
+// Titre	Public Methods
+//--
+
+//++
+template<class T>
+void SphereHealpix<T>::Resize(int_4 m)
+
+//    m is the "nside" of the Gorski algorithm
+//
+//    The total number of pixels will be Npix =  12*nside**2
+//
+//    nside MUST be a power of 2 (<= 8192)
+//--
+{
+  if (m<=0 || m> 8192) {
+    cout << "SphereHealpix : m hors bornes [0,8192], m= " << m << endl;
+    exit(1);
+  }
+  // verifier que m est une puissance de deux 
+  int x= m;
+  while (x%2==0) x/=2;
+  if(x != 1) 
+    {  
+      cout<<"SphereHealpix: m doit etre une puissance de deux, m= "<<m<<endl;
+      exit(1);
+    }
+  InitNul();
+  Pixelize(m);
+  SetThetaSlices();
+}
+
+template<class T>
+void  SphereHealpix<T>::Pixelize( int_4 m) 
+
+//    prépare la pixelisation Gorski (m a la même signification 
+//    que pour le constructeur)
+//
+//    
+{
+  // On memorise les arguments d'appel
+  nSide_= m;  
+
+  // Nombre total de pixels sur la sphere entiere
+  nPix_= 12*nSide_*nSide_;
+
+  // pour le moment les tableaux qui suivent seront ranges dans l'ordre 
+  // de l'indexation GORSKY "RING"
+  // on pourra ulterieurement changer de strategie et tirer profit
+  // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra 
+  // de pourquoi c'est faire
+
+  // Creation et initialisation du vecteur des contenus des pixels 
+  pixels_.ReSize(nPix_);
+  pixels_.Reset();
+
+  // solid angle per pixel   
+  omeg_= 4.0*Pi/nPix_;
+}
+
+template<class T>
+void SphereHealpix<T>::InitNul()
+//
+//    initialise à zéro les variables de classe 
+{
+  nSide_= 0;
+  nPix_ = 0;
+  omeg_ = 0.;
+//  pixels_.Reset();  -  Il ne faut pas mettre les pixels a zero si share !
+}
+
+/* --Methode-- */
+//++
+template<class T>
+int_4 SphereHealpix<T>::NbPixels() const
+
+//    Retourne le nombre de pixels du découpage
+//--
+{  
+  return(nPix_);
+}
+
+//++
+template<class T>
+uint_4 SphereHealpix<T>::NbThetaSlices() const 
+
+//    Return number of slices in theta direction on the  sphere
+//--
+{
+  uint_4 nbslices = uint_4(4*nSide_-1);
+  if (nSide_<=0) 
+    {
+      nbslices = 0;
+      throw PException(" sphere not pixelized, NbSlice=0 ");
+    }
+  return nbslices;
+}
+
+//++
+template<class T>
+void SphereHealpix<T>::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 
+
+//    For a theta-slice with index 'index', return : 
+//
+//    the corresponding "theta" 
+//
+//    a vector containing the phi's of the pixels of the slice
+//
+//    a vector containing the corresponding values of pixels 
+//
+//--
+{
+
+  if (index<0 || index >= NbThetaSlices()) 
+    {
+      // THROW(out_of_range("SphereHealpix::PIxVal Pixel index out of range")); 
+      cout << " SphereHealpix::GetThetaSlice : Pixel index out of range" <<endl;
+      throw RangeCheckError(" SphereHealpix::GetThetaSlice : Pixel index out of range");
+    }
+
+
+  int_4 iring= sliceBeginIndex_(index);
+  int_4 lring  =  sliceLenght_(index);
+
+   phi.ReSize(lring);
+   value.ReSize(lring);
+
+  double TH= 0.;
+  double FI= 0.;
+  for(int_4 kk = 0; kk < lring;kk++) 
+    {
+      PixThetaPhi(kk+iring,TH,FI);
+      phi(kk)= FI;
+      value(kk)= PixVal(kk+iring);
+    }
+  theta= TH;
+}
+//++
+//++
+
+template<class T>
+void SphereHealpix<T>::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const 
+
+//    For a theta-slice with index 'sliceIndex', return : 
+//
+//    the corresponding "theta" 
+//    the corresponding "phi" for first pixel of the slice 
+//
+//    a vector containing the indices of the pixels of the slice
+//    (equally distributed in phi)
+//
+//    a vector containing the corresponding values of pixels 
+//
+//--
+{
+
+  if (sliceIndex<0 || sliceIndex >= NbThetaSlices()) 
+    {
+      // THROW(out_of_range("SphereHealpix::PIxVal Pixel index out of range")); 
+      cout << " SphereHealpix::GetThetaSlice : Pixel index out of range" <<endl;
+      throw RangeCheckError(" SphereHealpix::GetThetaSlice : Pixel index out of range");
+    }
+  int_4 iring= sliceBeginIndex_(sliceIndex);
+  int_4 lring  =  sliceLenght_(sliceIndex);
+  pixelIndices.ReSize(lring);
+  value.ReSize(lring);
+
+  for(int_4 kk = 0; kk < lring;kk++) 
+    {
+      pixelIndices(kk)= kk+iring;
+      value(kk)= PixVal(kk+iring);
+    }
+  PixThetaPhi(iring, theta, phi0);
+}
+//++
+template<class T>
+void SphereHealpix<T>::SetThetaSlices()  
+
+//--
+{
+  sliceBeginIndex_.ReSize(4*nSide_-1);
+  sliceLenght_.ReSize(4*nSide_-1);
+  for (int sliceIndex=0; sliceIndex<  nSide_-1; sliceIndex++)
+    {
+      sliceBeginIndex_(sliceIndex)  = 2*sliceIndex*(sliceIndex+1);
+      sliceLenght_(sliceIndex) = 4*(sliceIndex+1);
+    }
+  for (int sliceIndex= nSide_-1; sliceIndex<  3*nSide_; sliceIndex++)
+    {
+      sliceBeginIndex_(sliceIndex)  = 2*nSide_*(2*sliceIndex-nSide_+1);
+      sliceLenght_(sliceIndex) = 4*nSide_;
+    }
+  for (int sliceIndex= 3*nSide_; sliceIndex< 4*nSide_-1; sliceIndex++)
+    {
+      int_4 nc= 4*nSide_-1-sliceIndex;
+      sliceBeginIndex_(sliceIndex)  = nPix_-2*nc*(nc+1); 
+      sliceLenght_(sliceIndex) = 4*nc;
+    }
+}
+
+/* --Methode-- */
+//++
+template<class T>
+T& SphereHealpix<T>::PixVal(int_4 k)
+
+//    Return value of  pixel with "RING" index k 
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    {
+      throw RangeCheckError("SphereHealpix::PIxVal Pixel index out of range");
+    }
+  return pixels_(k);
+}
+
+/* --Methode-- */
+//++
+template<class T>
+T const& SphereHealpix<T>::PixVal(int_4 k) const
+
+//    Return value of  pixel with "RING" index k 
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    {
+      throw RangeCheckError("SphereHealpix::PIxVal Pixel index out of range");
+    }
+  return *(pixels_.Data()+k);
+}
+
+//++
+template<class T>
+T& SphereHealpix<T>::PixValNest(int_4 k)
+
+//    Return value of  pixel with "NESTED" index k 
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    { 
+      throw RangeCheckError("SphereHealpix::PIxValNest Pixel index out of range");
+    }
+  return pixels_(nest2ring(nSide_,k));
+}
+//++
+
+template<class T>
+T const& SphereHealpix<T>::PixValNest(int_4 k) const
+
+//    Return value of  pixel with "NESTED" index k 
+//--
+{
+  if((k < 0) || (k >= nPix_)) 
+    { 
+      throw RangeCheckError("SphereHealpix::PIxValNest Pixel index out of range");
+  }
+  int_4 pix= nest2ring(nSide_,k);
+  return *(pixels_.Data()+pix);
+}
+
+/* --Methode-- */
+//++
+template<class T>
+bool SphereHealpix<T>::ContainsSph(double /*theta*/, double /*phi*/) const
+//--
+{
+return(true);
+}
+
+/* --Methode-- */
+//++
+template<class T>
+int_4 SphereHealpix<T>::PixIndexSph(double theta,double phi) const
+
+//    Return "RING" index of the pixel corresponding to 
+//    direction (theta, phi).
+//--
+{
+  return ang2pix_ring(nSide_,theta,phi);
+}
+
+//++
+template<class T>
+int_4 SphereHealpix<T>::PixIndexSphNest(double theta,double  phi) const
+
+//    Return "NESTED" index of the pixel corresponding to 
+//    direction (theta, phi).
+//--
+{
+  return ang2pix_nest(nSide_,theta,phi);
+}
+
+
+/* --Methode-- */
+//++
+template<class T>
+void SphereHealpix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
+
+//   Return (theta,phi) coordinates of middle of  pixel with "RING" index k
+//--
+{
+  pix2ang_ring(nSide_,k,theta,phi);
+}
+
+template <class T>
+T SphereHealpix<T>::SetPixels(T v)
+{
+pixels_.Reset(v);
+return(v);
+}
+
+//++
+template<class T>
+double SphereHealpix<T>::PixSolAngle(int_4 /*dummy*/) const
+//    Pixel Solid angle  (steradians)
+//    All the pixels have the same solid angle. The dummy argument is
+//    for compatibility with eventual pixelizations which would not 
+//    fulfil this requirement.
+//--
+{
+  return omeg_;
+}
+
+//++
+template<class T>
+void SphereHealpix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
+
+//   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
+//--
+{   
+  pix2ang_nest(nSide_,k,theta,phi);
+}
+
+//++
+template<class T>
+int_4 SphereHealpix<T>::NestToRing(int_4 k) const
+
+//    translation from NESTED index  into RING index 
+//
+//--
+{
+  return  nest2ring(nSide_,k);
+}
+
+//++
+template<class T>
+int_4 SphereHealpix<T>::RingToNest(int_4 k) const
+//
+//    translation from  RING index  into NESTED index 
+//
+//--
+{
+  return  ring2nest(nSide_,k);
+}
+
+
+template<class T>
+int_4 SphereHealpix<T>::nest2ring(int_4 nside, int_4 ipnest) const 
+{
+  /*
+    ====================================================
+    subroutine nest2ring(nside, ipnest, ipring)
+    ====================================================
+    c     conversion from NESTED to RING pixel number
+    ====================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+  
+  const PIXELS_XY& PXY= PIXELS_XY::instance();
+
+  int npix, npface, face_num, ncap, n_before;
+  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
+  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
+  int ns_max=8192;
+  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
+  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
+  
+  if(  nside<1 ||  nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12 *  nside* nside;
+  if( ipnest<0 || ipnest>npix-1 ) {
+    cout << "ipnest out of range" << endl;
+    exit(0);
+  }
+
+  ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
+  nl4   = 4* nside;
+  
+  //c     finds the face, and the number in the face
+  npface =  nside* nside;
+  //cccccc      ip = ipnest - 1         ! in {0,npix-1}
+  
+  face_num = ipnest/npface;//  ! face number in {0,11}
+  ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
+  //c     finds the x,y on the face (starting from the lowest corner)
+  //c     from the pixel number
+  ip_low=ipf%1024;                //   ! content of the last 10 bits
+  ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
+  ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
+  ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
+  
+  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
+  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
+  
+  //c     transforms this in (horizontal, vertical) coordinates
+  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
+  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
+  
+  //c     computes the z coordinate on the sphere
+  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
+  jr =  jrll[face_num]*nside - jrt - 1;
+  nr = nside;//                  ! equatorial region (the most frequent)
+  n_before = ncap + nl4 * (jr - nside);
+  kshift=(jr - nside)%2;
+  if( jr<nside ) {//then     ! north pole region
+    nr = jr;
+    n_before = 2 * nr * (nr - 1);
+    kshift = 0;
+  }
+  else if( jr>3*nside ) {//then ! south pole region
+    nr = nl4 - jr;
+    n_before = npix - 2 * (nr + 1) * nr;
+    kshift = 0;
+  }
+  
+  //c     computes the phi coordinate on the sphere, in [0,2Pi]
+  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
+  
+  if( jp>nl4 ) jp = jp - nl4;
+  if( jp<1 )   jp = jp + nl4;
+  
+  int aux=n_before + jp - 1;
+  return (n_before + jp - 1);// ! in {0, npix-1}
+}
+
+template<class T>
+int_4 SphereHealpix<T>::ring2nest(int_4 nside, int_4 ipring) const 
+{
+  /*
+    ==================================================
+    subroutine ring2nest(nside, ipring, ipnest)
+    ==================================================
+    c     conversion from RING to NESTED pixel number
+    ==================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  const PIXELS_XY& PXY= PIXELS_XY::instance();
+
+  double fihip, hip;
+  int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
+  int     kshift, face_num, nr;
+  int irn, ire, irm, irs, irt, ifm , ifp;
+  int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
+  int ns_max(8192);
+  
+  // coordinate of the lowest corner of each face
+  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
+  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12 * nside*nside;
+  if( ipring<0 || ipring>npix-1 ) {
+    cout << "ipring out of range" << endl;
+    exit(0);
+  }
+  
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  npix = 12*nside*nside;//      ! total number of points
+  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
+  ipring1 = ipring + 1;
+  
+  //c     finds the ring number, the position of the ring and the face number
+  if( ipring1<=ncap ) {//then
+    
+    hip   = ipring1/2.;
+    fihip = floor ( hip );
+    irn   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
+    iphi  = ipring1 - 2*irn*(irn - 1);
+    
+    kshift = 0;
+    nr = irn   ;//               ! 1/4 of the number of points on the current ring
+    face_num = (iphi-1) / irn;// ! in {0,3}
+  }
+  else if( ipring1<=nl2*(5*nside+1) ) {//then
+    
+    ip    = ipring1 - ncap - 1;
+    irn   = (int)floor( ip / nl4 ) + nside;//               ! counted from North pole
+    iphi  = (int)fmod(ip,nl4) + 1;
+    
+    kshift  = (int)fmod(irn+nside,2);//  ! 1 if irn+nside is odd, 0 otherwise
+    nr = nside;
+    ire =  irn - nside + 1;// ! in {1, 2*nside +1}
+    irm =  nl2 + 2 - ire;
+    ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
+    ifp = (iphi - irm/2 + nside -1) / nside;
+    if( ifp==ifm ) {//then          ! faces 4 to 7
+      face_num = (int)fmod(ifp,4) + 4;
+    }
+    else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
+      face_num = ifp;
+    }
+    else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
+      face_num = ifp + 7;
+    }
+  }
+  else {
+    
+    ip    = npix - ipring1 + 1;
+    hip   = ip/2.;
+    fihip = floor ( hip );
+    irs   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//  ! counted from South pole
+    iphi  = 4*irs + 1 - (ip - 2*irs*(irs-1));
+    
+    kshift = 0;
+    nr = irs;
+    irn   = nl4 - irs;
+    face_num = (iphi-1) / irs + 8;// ! in {8,11}
+  }
+  
+  //c     finds the (x,y) on the face
+  irt =   irn  - jrll[face_num]*nside + 1;//       ! in {-nside+1,0}
+  ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
+
+
+  if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
+  
+  ix =  (ipt - irt ) / 2;
+  iy = -(ipt + irt ) / 2;
+  
+  ix_low = (int)fmod(ix,128);
+  ix_hi  = ix/128;
+  iy_low = (int)fmod(iy,128);
+  iy_hi  = iy/128;
+  ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
+
+  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
+}
+
+template<class T>
+int_4 SphereHealpix<T>::ang2pix_ring(int_4 nside, double theta, double phi) const 
+{
+  /*
+    ==================================================
+    c     gives the pixel number ipix (RING) 
+    c     corresponding to angles theta and phi
+    c==================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  int nl2, nl4, ncap, npix, jp, jm, ipix1;
+  double  z, za, tt, tp, tmp;
+  int ir, ip, kshift;
+
+  double piover2(Pi/2.);
+  double twopi(2.*Pi);
+  double z0(2./3.);
+  int ns_max(8192);
+
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+
+  if( theta<0. || theta>Pi) {
+    cout << "theta out of range" << endl;
+    exit(0);
+  }
+
+  z = cos(theta);
+  za = fabs(z);
+  if( phi >= twopi)  phi = phi - twopi;
+  if (phi < 0.)     phi = phi + twopi;
+  tt = phi / piover2;//  ! in [0,4)
+
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  ncap  = nl2*(nside-1);// ! number of pixels in the north polar cap
+  npix  = 12*nside*nside;
+
+  if( za <= z0 ) {
+    
+    jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of  ascending edge line 
+    jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
+
+    ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
+    kshift = 0;
+    if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
+    
+    ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
+    if( ip>nl4 ) ip = ip - nl4;
+    
+    ipix1 = ncap + nl4*(ir-1) + ip ;
+  }
+  else {
+
+    tp = tt - floor(tt);//      !MOD(tt,1.d0)
+    tmp = sqrt( 3.*(1. - za) );
+
+    jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
+    jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
+
+    ir = jp + jm + 1;//        ! ring number counted from the closest pole
+    ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
+    if( ip>4*ir ) ip = ip - 4*ir;
+      
+    ipix1 = 2*ir*(ir-1) + ip;
+      if( z<=0. ) {
+	ipix1 = npix - 2*ir*(ir+1) + ip;
+      }
+  }
+    return (ipix1 - 1);// ! in {0, npix-1}
+}
+
+template<class T>
+int_4 SphereHealpix<T>::ang2pix_nest(int_4 nside, double theta, double phi) const 
+{
+  /*
+    ==================================================
+    subroutine ang2pix_nest(nside, theta, phi, ipix)
+    ==================================================
+    c     gives the pixel number ipix (NESTED) 
+    c     corresponding to angles theta and phi
+    c
+    c     the computation is made to the highest resolution available (nside=8192)
+    c     and then degraded to that required (by integer division)
+    c     this doesn't cost more, and it makes sure 
+    c     that the treatement of round-off will be consistent 
+    c     for every resolution
+    ==================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+  
+  const PIXELS_XY& PXY= PIXELS_XY::instance();
+
+  double    z, za, z0, tt, tp, tmp;
+  int face_num,jp,jm;
+  int ifp, ifm;
+  int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
+  double piover2(Pi/2.), twopi(2.*Pi);
+  int ns_max(8192);
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  if( theta<0 || theta>Pi ) {
+    cout << "theta out of range" << endl;
+    exit(0);
+  }
+  z  = cos(theta);
+  za = fabs(z);
+  z0 = 2./3.;
+  if( phi>=twopi ) phi = phi - twopi;
+  if( phi<0. )    phi = phi + twopi;
+  tt = phi / piover2;// ! in [0,4[
+  if( za<=z0 ) { // then ! equatorial region
+    
+    //(the index of edge lines increase when the longitude=phi goes up)
+    jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
+    jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
+    
+    //c        finds the face
+    ifp = jp / ns_max;//  ! in {0,4}
+    ifm = jm / ns_max;
+    if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then  ! faces 4 to 7
+    else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
+    else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
+    
+    ix = (int)fmod(jm, ns_max);
+    iy = ns_max - (int)fmod(jp, ns_max) - 1;
+  }
+  else { //! polar region, za > 2/3
+    
+    ntt = (int)floor(tt);
+    if( ntt>=4 ) ntt = 3;
+    tp = tt - ntt;
+    tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
+    
+    //(the index of edge lines increase when distance from the closest pole goes up)
+    jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
+    jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
+    jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
+    jm = (int)min(ns_max-1, jm);
+    
+    // finds the face and pixel's (x,y)
+    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;
+    }
+  }
+  
+  ix_low = (int)fmod(ix,128);
+  ix_hi  =     ix/128;
+  iy_low = (int)fmod(iy,128);
+  iy_hi  =     iy/128;
+  ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
+  //  ipf = ipf / pow(ns_max/nside,2.);//  ! in {0, nside**2 - 1}
+  //  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
+  // $CHECK$  Reza 25/10/99 , pow remplace par *
+  ipf = ipf / ((ns_max/nside)*(ns_max/nside));
+  return (ipf + face_num*nside*nside);
+}
+
+template<class T>
+void SphereHealpix<T>::pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const {
+  /*
+    ===================================================
+    c     gives theta and phi corresponding to pixel ipix (RING) 
+    c     for a parameter nside
+    ===================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+  
+  int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
+  double  fact1, fact2, fodd, hip, fihip;
+  
+  int ns_max(8192);
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12*nside*nside;      // ! total number of points
+  if( ipix<0 || ipix>npix-1 ) {
+    cout << "ipix out of range" << endl;
+    exit(0);
+  }
+  
+  ipix1 = ipix + 1; // in {1, npix}
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
+  fact1 = 1.5*nside;
+  fact2 = 3.0*nside*nside;
+  
+  if( ipix1 <= ncap ) {  //! North Polar cap -------------
+    
+    hip   = ipix1/2.;
+    fihip = floor(hip);
+    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
+    iphi  = ipix1 - 2*iring*(iring - 1);
+    
+    theta = acos( 1. - iring*iring / fact2 );
+    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
+    //    cout << theta << " " << phi << endl;
+  }
+  else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
+    
+    ip    = ipix1 - ncap - 1;
+    iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
+    iphi  = ip%nl4 + 1;
+    
+    fodd  = 0.5 * (1 + (iring+nside)%2 );//  ! 1 if iring+nside is odd, 1/2 otherwise
+    theta = acos( (nl2 - iring) / fact1 );
+    phi   = ((double)iphi - fodd) * Pi /(2.*nside);
+  }
+  else {//! South Polar cap -----------------------------------
+    
+    ip    = npix - ipix1 + 1;
+    hip   = ip/2.;
+    fihip = floor(hip);
+    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//     ! counted from South pole
+    iphi  = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
+    
+    theta = acos( -1. + iring*iring / fact2 );
+    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
+    //    cout << theta << " " << phi << endl;
+  }
+}
+
+template<class T>
+void SphereHealpix<T>::pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const {
+  /*
+    ==================================================
+    subroutine pix2ang_nest(nside, ipix, theta, phi)
+    ==================================================
+    c     gives theta and phi corresponding to pixel ipix (NESTED) 
+    c     for a parameter nside
+    ==================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  const PIXELS_XY& PXY= PIXELS_XY::instance();
+    
+  int npix, npface, face_num;
+  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
+  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
+  double z, fn, fact1, fact2;
+  double piover2(Pi/2.);
+  int ns_max(8192);
+          
+  // ! coordinate of the lowest corner of each face
+  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
+  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12 * nside*nside;
+  if( ipix<0 || ipix>npix-1 ) {
+    cout << "ipix out of range" << endl;
+    exit(0);
+  }
+  
+  fn = 1.*nside;
+  fact1 = 1./(3.*fn*fn);
+  fact2 = 2./(3.*fn);
+  nl4   = 4*nside;
+
+  //c     finds the face, and the number in the face
+  npface = nside*nside;
+  
+  face_num = ipix/npface;//  ! face number in {0,11}
+  ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
+  
+  //c     finds the x,y on the face (starting from the lowest corner)
+  //c     from the pixel number
+  ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
+  ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
+  ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
+  ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
+  
+  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
+  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
+  
+  //c     transforms this in (horizontal, vertical) coordinates
+  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
+  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
+  
+  //c     computes the z coordinate on the sphere
+  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
+  jr =  jrll[face_num]*nside - jrt - 1;
+  nr = nside;//                  ! equatorial region (the most frequent)
+  z  = (2*nside-jr)*fact2;
+  kshift = (int)fmod(jr - nside, 2);
+  if( jr<nside ) { //then     ! north pole region
+    nr = jr;
+    z = 1. - nr*nr*fact1;
+    kshift = 0;
+  }
+  else {
+    if( jr>3*nside ) {// then ! south pole region
+      nr = nl4 - jr;
+      z = - 1. + nr*nr*fact1;
+      kshift = 0;
+    }
+  }
+  theta = acos(z);
+  
+  //c     computes the phi coordinate on the sphere, in [0,2Pi]
+  //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
+  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
+  if( jp>nl4 ) jp = jp - nl4;
+  if( jp<1 )   jp = jp + nl4;
+  phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
+}
+
+
+
+template <class T>
+void SphereHealpix<T>::print(ostream& os) const
+{
+  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
+  //
+  os << " nSide_ = " << nSide_  << endl;
+  os << " nPix_   = " << nPix_   << endl;
+  os << " omeg_ = " << omeg_   << endl;
+
+  os << " content of pixels : ";
+  for(int i=0; i < nPix_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os <<  pixels_(i) <<", ";
+    }
+  os << endl;
+
+  os << endl;
+  //const PIXELS_XY& PXY= PIXELS_XY::instance();
+
+  //os << endl;  os << " contenu des tableaux conversions "<<endl;
+  //for(int i=0; i < 5; i++)
+  //  {
+  //   os<<PXY.pix2x_(i)<<", "<<PXY.pix2y_(i)<<", "<<PXY.x2pix_(i)<<", "<<PXY.y2pix_(i)<<endl;
+  // }
+  os << endl;
+
+}
+
+
+
+//*******************************************************************
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template SphereHealpix<uint_2>
+#pragma define_template SphereHealpix<r_8>
+#pragma define_template SphereHealpix<r_4>
+#pragma define_template SphereHealpix< complex<r_4> >
+#pragma define_template SphereHealpix< complex<r_8> >
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class SphereHealpix<uint_2>;
+template class SphereHealpix<r_8>;
+template class SphereHealpix<r_4>;
+template class SphereHealpix< complex<r_4> >;
+template class SphereHealpix< complex<r_8> >;
+#endif
+
Index: /trunk/SophyaLib/SkyMap/spherehealpix.h
===================================================================
--- /trunk/SophyaLib/SkyMap/spherehealpix.h	(revision 843)
+++ /trunk/SophyaLib/SkyMap/spherehealpix.h	(revision 843)
@@ -0,0 +1,264 @@
+#ifndef SPHEREHEALPIX_SEEN
+#define SPHEREHEALPIX_SEEN
+
+#include "sphericalmap.h"
+#include "tvector.h"
+#include "ndatablock.h"
+
+#include "anydataobj.h"
+#include "ppersist.h"
+
+
+namespace SOPHYA {
+
+
+// ***************** CLASSE SphereHealpix *****************************
+
+//!  class SphereHealpix
+/*!
+   Pixelisation Gorski  
+
+
+    -----------------------------------------------------------------------
+     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
+    -----------------------------------------------------------------------
+
+    the sphere is split in 12 diamond-faces containing nside**2 pixels each
+
+    the numbering of the pixels (in the nested scheme) is similar to
+    quad-cube
+    In each face the first pixel is in the lowest corner of the diamond
+
+    the faces are                    (x,y) coordinate on each face
+\verbatim
+        .   .   .   .   <--- North Pole
+       / \ / \ / \ / \                          ^        ^     
+      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      / 
+       \ / \ / \ / \ /                        y   \    /  x  
+      4 . 5 . 6 . 7 . 4 <--- equator               \  /      
+       / \ / \ / \ / \                              \/      
+      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner 
+       \ / \ / \ / \ /      
+        .   .   .   .   <--- South Pole
+\endverbatim
+    phi:0               2Pi                                 
+
+   in the ring scheme pixels are numbered along the parallels
+   the first parallel is the one closest to the north pole and so on  
+   on each parallel, pixels are numbered starting from the one closest
+    to phi = 0
+
+    nside MUST be a power of 2 (<= 8192)
+
+*/
+
+
+template<class T>
+class FIO_SphereHealpix;
+
+template<class T>
+class FITS_SphereHealpix;
+
+template<class T>
+class SphereHealpix : public SphericalMap<T>
+{
+
+
+  friend class FIO_SphereHealpix<T>;
+  friend class FITS_SphereHealpix<T>;
+
+public :
+
+SphereHealpix();
+/*!    
+  m is the "nside" of the Gorski algorithm
+
+  The total number of pixels will be Npix =  12*nside**2
+
+  nside MUST be a power of 2 (<= 8192)
+*/
+SphereHealpix(int_4 m);
+SphereHealpix(const SphereHealpix<T>& s, bool share=false);
+//!	Destructor
+virtual ~SphereHealpix();
+
+/*! Setting blockdata to temporary (see ndatablock documentation) */
+inline virtual void SetTemp(bool temp=false) const 
+  {
+    pixels_.SetTemp(temp);
+    sliceBeginIndex_.SetTemp(temp);
+    sliceLenght_.SetTemp(temp);
+  };
+// ------------------ Definition of PixelMap abstract methods
+
+/* Nombre de pixels du decoupage */
+/*!    Return number of  pixels of the  splitting */
+virtual int_4 NbPixels() const;
+
+/* Valeur du contenu du pixel d'indice "RING" k  */
+/*!    Return value of  pixel with "RING" index k */
+virtual T& PixVal(int_4 k);
+virtual T const& PixVal(int_4 k) const;
+
+/* Nombre de tranches en theta */ 
+/*!    Return number of slices in theta direction on the  sphere */
+uint_4 NbThetaSlices() const;
+/*!   For a theta-slice with index 'index', return : 
+
+   the corresponding "theta" 
+
+    a vector containing the phi's of the pixels of the slice
+
+    a vector containing the corresponding values of pixels 
+*/
+virtual void GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const;
+/*!   For a theta-slice with index 'index', return : 
+
+   the corresponding "theta" 
+
+   the corresponding "phi" for first pixel of the slice 
+
+    a vector containing indices of the pixels of the slice
+
+   (equally distributed in phi)
+
+    a vector containing the corresponding values of pixels 
+*/
+virtual void GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const ;
+
+/* Return true if teta,phi in map  */
+virtual bool ContainsSph(double theta, double phi) const;
+/* Indice "RING" du pixel vers lequel pointe une direction definie par 
+ses  coordonnees spheriques */                                    
+/*!  Return "RING" index of the pixel corresponding to direction (theta, phi).
+ */
+virtual int_4 PixIndexSph(double theta,double phi) const;
+
+/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
+virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;
+
+/*! Set all pixels to value v */
+virtual T SetPixels(T v);
+
+/* Pixel Solid angle  (steradians) */
+/*!    Pixel Solid angle  (steradians)
+
+    All the pixels have the same solid angle. The dummy argument is
+    for compatibility with eventual pixelizations which would not 
+   fulfil this requirement.
+*/
+virtual double PixSolAngle(int_4 dummy=0) const;
+
+/*  Acces to the DataBlock  */
+inline       NDataBlock<T>& DataBlock()       {return pixels_;}
+inline const NDataBlock<T>& DataBlock() const {return pixels_;}
+
+// --------------- Specific methods 
+
+/*!    
+  m is the "nside" of the Gorski algorithm
+
+  The total number of pixels will be Npix =  12*nside**2
+
+  nside MUST be a power of 2 (<= 8192)
+*/
+virtual void Resize(int_4 m);
+
+// pour l'instant le tableau est ordonne selon RING, uniquement
+inline virtual char* TypeOfMap() const {return "RING";};
+
+
+/* Valeur du contenu du pixel d'indice "NEST" k */
+/*!    Return value of  pixel with "NESTED" index k */
+virtual T& PixValNest(int_4 k);
+/*!    Return value of  pixel with "NESTED" index k */
+virtual T const& PixValNest(int_4 k) const;
+
+/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
+ses  coordonnees spheriques */                                    
+/*! Return "NESTED" index of the pixel corresponding to direction (theta, phi).
+ */
+virtual int_4 PixIndexSphNest(double theta,double phi) const;
+
+/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
+/*!   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
+ */
+virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;
+
+/* algorithme de pixelisation */
+void Pixelize(int_4); 
+
+/* convertit index nested en ring  */
+/*!    translation from NESTED index  into RING index */
+int_4 NestToRing(int_4) const;
+
+/* convertit index ring en nested" */
+/*!    translation from  RING index  into NESTED index */
+int_4 RingToNest(int_4) const;
+
+
+/* retourne la valeur du parametre Gorski */
+inline virtual int_4 SizeIndex() const {return(nSide_);}
+
+/* impression */ 
+void print(ostream& os) const;
+         
+private :
+
+// ------------- méthodes internes ----------------------
+void InitNul();
+void SetThetaSlices(); 
+
+int_4  nest2ring(int_4 nside,int_4 ipnest) const;
+int_4  ring2nest(int_4 nside,int_4 ipring) const;
+
+int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
+int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
+void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
+void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
+inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 
+  {
+    nSide_= nside;
+    nPix_= nbpixels;
+    omeg_= solangle;
+  }
+
+// ------------- variables internes -----------------------
+
+int_4 nSide_;
+int_4 nPix_;
+double omeg_;
+
+NDataBlock<T> pixels_;
+NDataBlock<int> sliceBeginIndex_;
+NDataBlock<int> sliceLenght_;
+
+};
+
+
+//////////////////////////////////////////////////////////////////////////
+//
+// ------------- Classe PIXELS_XY -----------------------
+//
+class PIXELS_XY
+{
+
+public :
+
+static PIXELS_XY& instance();
+
+NDataBlock<int_4> pix2x_;
+NDataBlock<int_4> pix2y_;
+NDataBlock<int_4> x2pix_;
+NDataBlock<int_4> y2pix_;
+
+private :
+
+PIXELS_XY();
+void mk_pix2xy();
+void mk_xy2pix();
+};
+
+} // Fin du namespace
+
+#endif
