#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_;
}
//++
template<class T>
SphereHEALPix<T>::SphereHEALPix(const SphereHEALPix<T>& s)
  : pixels_(s.pixels_), sliceBeginIndex_(s.sliceBeginIndex_), 
				sliceLenght_(s.sliceLenght_)
//    copy constructor
//--
{
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);

  nSide_= s.nSide_;
  nPix_ = s.nPix_;
  omeg_ = s.omeg_;
}

template<class T>
void SphereHEALPix<T>::CloneOrShare(const SphereHEALPix<T>& a)
{
  pixels_.CloneOrShare(a.pixels_);
  sliceBeginIndex_.CloneOrShare(a.sliceBeginIndex_);
  sliceLenght_.CloneOrShare(a.sliceLenght_);
}

template<class T>
SphereHEALPix<T>& SphereHEALPix<T>::Set(const SphereHEALPix<T>& a)
   {
  if (this != &a) 
    { 
      nSide_= a.nSide_;
     nPix_ = a.nPix_;
     omeg_ = a.omeg_;
     CloneOrShare(a);
     if (mInfo_) delete mInfo_;
     mInfo_ = NULL;
     if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
    }
  return(*this);
   }


//++
// 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