source: Sophya/trunk/SophyaLib/Samba/spheregorski.cc@ 228

//
#include "spheregorski.h"
#include "strutil.h"
extern "C" {
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
}
      
extern "C" {
  //void ang2pix_ring_(int&,double&,double&,int&);
  //void pix2ang_ring_(int&, int&, double&, double&);
//void ang2pix_nest_(int&,double&,double&,int&);
  //void pix2ang_nest_(int&, int&, double&, double&);
//void nest2ring_(int&, int&, int&);
//void ring2nest_(int&, int&, int&);
void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,
              float*,float*,float*);
void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,
              double*, double*,double*,double*,double*,float*);
void  input_map_(char*,float*,int&);
void  ecrire_fits_(int&,int&,int&,float*,char*,char*,int&,float&,float&);

}
//*******************************************************************
//++ 
// 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 DOIT OBLIGATOIREMENT ETRE UNE PUISSANCE DE 2 (<= 8192)
//--
//++
//
// Links        Parents
//
//    SphericalMap 
//-- 
//++
//
// Links        Descendants
//
//     
//-- 

/* --Methode-- */
//++
// Titre        Constructeurs
//--
//++

SphereGorski::SphereGorski()

//--
{
  InitNul();
}

/* --Methode-- */
//++
SphereGorski::SphereGorski(char* flnm)

//    Constructeur : charge une image à partir d'un fichier 
//--
{
  InitNul();
  PInPersist s(flnm);
  Read(s);
}

//++
SphereGorski::SphereGorski(int_4 m)

//    Constructeur : m est la variable nside de l'algorithme de Gorski
//    le nombre total de pixels sera Npix =  12*nside**2
//    m DOIT OBLIGATOIREMENT ETRE UNE PUISSANCE DE 2 (<= 8192)
//--
{
  //printf(" initialisation par defaut SphereGorski \n");
  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_4 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);
  if (pix2x_==NULL) {
    pix2x_=new int[1024];
    pix2y_=new int[1024];
    mk_pix2xy(pix2x_,pix2y_); 
  }
  if (x2pix_==NULL) {
    x2pix_=new int[128];
    y2pix_=new int[128];
    mk_xy2pix(x2pix_,y2pix_); 
  }
}






//++
// Titre        Destructeur
//--
//++
SphereGorski::~SphereGorski()

//--
{
  Clear();
}

//++
// Titre        Méthodes
//--




void  SphereGorski::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 er NEST) : tout dependra 
// de pourquoi c'est faire

  // Creation et initialisation du vecteur des contenus des pixels 
  mPix_ = new r_8[nPix_];
  for(int i=0; i<nPix_; i++)  mPix_[i] = 0.;

  // solid angle per pixel   
  omeg_=4*Pi/nPix_;
}

void SphereGorski::InitNul()
//
//    initialise à zéro les variables de classe 
{
  nlmax_=0;
  nmmax_=0;
  iseed_=0;
  fwhm_=0.;
  quadrupole_=0.;
  sym_cut_deg_=0.;
  for (int k=0; k<128;k++) powFile_[k]=' ';
  nSide_=0;
  nPix_ =0;
  mPix_ = NULL;
  pix2x_=NULL;
  pix2y_=NULL;
  x2pix_=NULL;
  y2pix_=NULL;
  pix2xy_=NULL;
}
/* --Methode-- */
void SphereGorski::Clear()

{
  if (mPix_)   delete[] mPix_;
  if (pix2x_)  delete[] pix2x_;
  if (pix2y_)  delete[] pix2y_;
  if (x2pix_)  delete[] x2pix_;
  if (y2pix_)  delete[] y2pix_;
  if (pix2xy_) delete pix2xy_;
}
//++
void SphereGorski::WriteSelf(POutPersist& s) const

//    créer un fichier image
//--
{
  char strg[256];
  if (mInfo_) sprintf(strg, "SphereGorski: NSlices=%6d  NPix=%9d HasInfo",  
                   (int_4)nSide_, (int_4)nPix_);
  else
    sprintf(strg, "SphereGorski: nSide=%6d  nPix=%9d ", 
          (int_4)nSide_, (int_4)nPix_);
  s.PutLine(strg);
  if (mInfo_)  mInfo_->Write(s);
  s.PutI4(nlmax_);
  s.PutI4(nmmax_);
  s.PutI4(iseed_);
  s.PutI4(nSide_);
  s.PutI4(nPix_);
  s.PutR4(fwhm_);
  s.PutR4(quadrupole_);
  s.PutR4(sym_cut_deg_);
  s.PutR8(omeg_);
  s.PutR8s(mPix_, nPix_);
  s.PutLine(powFile_);
  return;
}

/* --Methode-- */
//++
void SphereGorski::ReadSelf(PInPersist& s)

//    relit un fichier d'image
//--
{
  Clear();
  char strg[256];
  s.GetLine(strg, 255);
// Pour savoir s'il y avait un DVList Info associe
  bool hadinfo = false;
  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
  if (hadinfo) {    // Lecture eventuelle du DVList Info
    if (mInfo_ == NULL)  mInfo_ = new DVList;
    mInfo_->Read(s);
  }
  s.GetI4(nlmax_);
  s.GetI4(nmmax_);
  s.GetI4(iseed_);
  s.GetI4(nSide_);
  s.GetI4(nPix_);
  s.GetR4(fwhm_);
  s.GetR4(quadrupole_);
  s.GetR4(sym_cut_deg_);
  mPix_=new r_8[nPix_];
  s.GetR8(omeg_);
  s.GetR8s(mPix_, nPix_);
  s.GetLine(powFile_, 127);

  return;
}


/* --Methode-- */
//++
void SphereGorski::ReadFits(char flnm[])

//    remplit la sphere a partir d'un fichier FITS
//--
{
  strip(flnm,'B',' ');
  if (access(flnm,F_OK) != 0) {perror(flnm); exit(1);}
  if(!nPix_) {
    cout << " ReadFits : SphereGorski non pixelisee " << endl;
    exit(1);
  }
  int npixtot=nPix_;

// quand map et mPix_ auront le meme type, map ne sera plus necessaire

  float* map=new float[12*nSide_*nSide_];
  input_map_(flnm,map,npixtot);
  // Remplissage de la sphère 

  for(int  j=0; j<nPix_; j++ ) mPix_[j]=(r_8)map[j];
  delete [] map;
  

}

/* --Methode-- */
//++
void SphereGorski::WriteFits(char flnm[])

//    ecrit la sphere sur un fichier FITS

//--
{
  strip(flnm,'B',' ');
  if (access(flnm,F_OK) == 0) {
    cout << " SphereGorski::WriteFits : le fichier existe deja" << endl; 
    exit(1);
  }
  //else cout << "un fichier sera cree " << endl;
  if(!nPix_) {
    cout << " WriteFits : SphereGorski non pixelisee " << endl;
    exit(1);
  }


  char infile[128];
  for (int k=0; k< 128; k++) infile[k]=powFile_[k];

// quand map et mPix_ auront le meme type, map ne sera plus necessaire

  float* map=new float[12*nSide_*nSide_];

  for(int  j=0; j<nPix_; j++ ) map[j]=(float)mPix_[j];
  int nlmax=nlmax_;
  int nsmax=nSide_;
  int nmmax=nmmax_;
  int iseed=iseed_;
  float fwhm=fwhm_;
  float quadrupole=quadrupole_;
  ecrire_fits_(nlmax,nsmax,nmmax,map,infile,flnm,iseed,fwhm,quadrupole);
  delete [] map;
  
}

/* --Methode-- */
//++
int_4 SphereGorski::NbPixels() const

//    Retourne le nombre de pixels du découpage
//--
{
  return(nPix_);
}

//++
int_4 SphereGorski::NbThetaSlices() const 

//    Retourne le nombre de tranches en theta sur la sphere
//--
{return int_4(4*nSide_-1);}


//++
void  SphereGorski::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const 

//    Retourne, pour la tranche en theta d'indice 'index' le theta 
//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
//    correspondantes
//--
{

cout << "entree GetThetaSlice, couche no " << index << endl;
  if (index<0 || index > NbThetaSlices()) {
    // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
    cout << " SphereGorski::GetThetaSlice : exceptions  a mettre en place" <<endl;
  THROW(rangeCheckErr);
  }
  int_4 iring=0;
  int lring=0;
  if (index<nSide_-1) {
    iring=2*index*(index+1);
    lring=4*(index+1);
  }else
    if (index<3*nSide_) {
      iring=2*nSide_*(2*index-nSide_+1);
      lring=4*nSide_;
    }else
      {
        int nc=4*nSide_-1-index;
        iring=nPix_-2*nc*(nc+1);
        lring=4*nc;
      }
  phi.Realloc(lring);
  value.Realloc(lring);
  float T=0.;
  float F=0.;
  for (int kk=0; kk<lring;kk++) {
    PixThetaPhi(kk+iring,T,F);
    phi(kk)=F;
    value(kk)=PixVal(kk+iring);
  }
  theta=T;

}



/* --Methode-- */
//++
r_8& SphereGorski::PixVal(int_4 k)

//    Retourne la valeur du contenu du pixel d'indice "RING" k 
//--
{
  if ( (k<0) || (k >= nPix_) ) {
    // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
    cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
  THROW(rangeCheckErr);
  }
  return(mPix_[k]);
}

/* --Methode-- */
//++
r_8 const& SphereGorski::PixVal(int_4 k) const

//    Retourne la valeur du contenu du pixel d'indice "RING" k 
//--
{
  if ( (k<0) || (k >= nPix_) ) {
    //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
    cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
  THROW(rangeCheckErr);
  }
  return(mPix_[k]);
}


//++
r_8& SphereGorski::PixValNest(int_4 k)

//    Retourne la valeur du contenu du pixel d'indice "NESTED" 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 mPix_[nest2ring(nSide_,k)];
}
//++

r_8 const& SphereGorski::PixValNest(int_4 k) const

//    Retourne la valeur du contenu du pixel d'indice "NESTED" 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 mPix_[pix];
}



/* --Methode-- */
//++
int_4 SphereGorski::PixIndexSph(r_4 theta, r_4 phi) const

//    Retourne l'indice "RING" du pixel vers lequel pointe une direction 
//    définie par ses coordonnées sphériques
//--
{
  return ang2pix_ring(nSide_, double(theta), double(phi));
}

//++
int_4 SphereGorski::PixIndexSphNest(r_4 theta, r_4 phi) const

//    Retourne l'indice NESTED" du pixel vers lequel pointe une direction 
//    définie par ses coordonnées sphériques
//--
{
  return ang2pix_nest(nSide_, double(theta), double(phi));
}


/* --Methode-- */
//++
void SphereGorski::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const

//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice "RING" k
//--
{
  double t;
  double p;
  pix2ang_ring(nSide_,k, t,  p);
  teta=(r_4)t;
  phi=(r_4)p;
}

//++
r_8       SphereGorski::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_;
}


//++
void SphereGorski::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const

//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice 
//    NESTED k
//--
{   
  double t;
  double p;
  pix2ang_nest(nSide_, k, t, p);
  teta=(r_4)t;
  phi=(r_4)p;
}

//++
int_4 SphereGorski::NestToRing(int_4 k) 

//    conversion d'index NESTD en un index RING
//
//--
{
  return  nest2ring(nSide_,k);
}
//++
int_4 SphereGorski::RingToNest(int_4 k) 
//
//    conversion d'index RING en un index NESTED
//
//--
{
  return  ring2nest(nSide_,k);
}
//++
void SphereGorski::anharm(int nlmax, float sym_c,float* powspec) 
//
//    analyse en harmoniques spheriques des valeurs des pixels de la 
//    sphere : appel du module anafast (Gorski-Hivon)
//
//    "nlmax"   : multipole maximum, nlmax <= 2*nsmax (cf. Nyquist)
//
//    "sym c"   : coupure symetrique autour de l'equateur (degres)
//
//    "powspec" : tableau resultat (a reserver avant l'appel) de C(l) 
//    (spectre de puissance)
//
//--
//
// Pb a resoudre : dans cette classe les valeurs de pixel sont "double"
// dans anafast le tableau correspondant est "float"
// pour l'instant on duplique les tableaux, il faudra decider quelque chose
//
{ 
  if (nlmax > 2*nSide_) {
    cout << " anharm : nlmax= " << nlmax << 
            " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
    exit(1);
  }
  else {
    nlmax_=nlmax;
    nmmax_=nlmax_;
  }
  sym_cut_deg_=sym_c;
  float* map=new float[nPix_];
  for (int k=0; k<nPix_; k++) map[k]=(float)mPix_[k];
  int nsmax=nSide_;
  int nmmax=nmmax_;
  double sc=(double)sym_cut_deg_;
  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
  if (powspec==NULL) {
    
    cout << 
      " anharm : un tableau de C_l doit etre alloue avant appel " << endl;
    exit(1);
  }
  float* phas_n=new float[2*(nmmax+1)];
  float* phas_s=new float[2*(nmmax+1)];
  float*  dataw =new float[16*nsmax];
  float*  work =new float[16*nsmax];

  anafast_(nsmax,nlmax,nmmax,sc,map,alm_T, powspec,phas_n,phas_s,dataw,work);
  quadrupole_=powspec[2];
  delete [] map;
  delete [] alm_T;
  delete [] phas_n;
  delete [] phas_s;
  delete [] dataw;
  delete [] work;
  
}


//++
void SphereGorski::synharm(int nlmax, int iseed,float fwhm, float* powspec) 
//
//    synthese  des valeurs des pixels de la sphere par l'intermediaire 
//    des coefficients en harmoniques spheriques reconstitues apartir d'un 
//    spectre en puissance : appel du module synfast (Gorski-Hivon)
//
//    powspec est un tableau (a fournir) de C(l) (spectre de puissance)
//    Ce tableau doit contenir les valeur de C(l) par ordre 
//    SEQUENTIEL de l (de l=0 a l=nlmax). IL SERA MODIFIE PAR L'ALGORITHME
//
//    nlmax : multipole maximum (nlmax <= 2*nsmax (cf. Nyquist)
//    iseed : initialisation generation aleatoire (negatif, suggere : -1)
//    fwhm  : largeur totale a mi-hauteur (minutes d'arc, >=0, ex: 5)
//--
// Pb a resoudre : dans cette classe les valeurs de pixel sont "double"
// dans anafast le tableau correspondant est "float"
// pour l'instant on duplique les tableaux, il faudra decider quelque chose

{
  if (nlmax > 2*nSide_) {
    cout << " sphereGorski::synharm: nlmax= " << nlmax << 
            " doit etre <= 2*nsmax (cf. Nyquist), soit : " << 2*nSide_ << endl;
    exit(1);
  }
  else {
    nlmax_=nlmax;
    nmmax_=nlmax_;
    quadrupole_=powspec[2];
  }
  if (powspec==NULL) {
    
    cout << 
      "sphereGorski::synharm : un tableau de C_l doit etre alloue avant appel"
         << endl;
    exit(1);
  }
  iseed_=iseed;
  fwhm_ =fwhm; 
  float* map=new float[nPix_];
  int nsmax=nSide_;
  int nmmax=nmmax_;
  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
  

//    tableaux de travail 
  double* b_north=new double[2*(2*nmmax+1)];
  double* b_south=new double[2*(2*nmmax+1)];
  double* bw=new double[2*4*nsmax];
  double* data=new double[2*4*nsmax];
  double* work=new double[2*4*nsmax];
  float* lread=new float[nlmax+1];
  synfast_(nsmax,nlmax,nmmax,iseed,fwhm, map,alm_T, powspec,
                  b_north,b_south,bw,data,work,lread);
  for (int k=0; k<nPix_; k++) mPix_[k]=(double)map[k];
  delete [] map;
  delete [] alm_T;
  delete [] b_north;
  delete [] b_south;
  delete [] bw; 
  delete [] data;
  delete [] work;
  delete [] lread;

}

int SphereGorski::nest2ring(int nside, int ipnest) const {
  /*
    c=======================================================================
    subroutine nest2ring(nside, ipnest, ipring)
    c=======================================================================
    c     conversion from NESTED to RING pixel number
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)

      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*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + 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}

}

void SphereGorski::mk_pix2xy(int *pix2x,int *pix2y) {
  /*
    c=======================================================================
    subroutine mk_pix2xy
    c=======================================================================
    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
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)

      int kpix, jpix, IX, IY, IP, ID;
      for (int i=0;i<1023;i++ ) pix2x[i]=0;

      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
      }
}

int SphereGorski::ring2nest(int nside, int ipring) const {
  /*
    c=======================================================================
    subroutine ring2nest(nside, ipring, ipnest)
    c=======================================================================
    c     conversion from RING to NESTED pixel number
    c=======================================================================
    */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)

  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 = (int)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 =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)
    + (x2pix_[ix_low]+y2pix_[iy_low]);

  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
}

void SphereGorski::mk_xy2pix(int *x2pix, int *y2pix) {
  /*
    c=======================================================================
    subroutine mk_xy2pix
    c=======================================================================
    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}
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)

  int K,IP,I,J,ID;

  for( int i(0);i<127;i++ ) x2pix[i]=0;
  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;
    }
  }     
  //c      endif
  
}

int SphereGorski::ang2pix_ring(int nside, double theta, double phi) const {
    /*
    c=======================================================================
    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}
  
}

int SphereGorski::ang2pix_nest(int nside, double theta, double phi) const {
  /*
    c=======================================================================
    subroutine ang2pix_nest(nside, theta, phi, ipix)
    c=======================================================================
    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
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)
    
      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 =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)+ (x2pix_[ix_low]+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}
}


void SphereGorski::pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const {
  /*
    c=======================================================================
    c     gives theta and phi corresponding to pixel ipix (RING) 
    c     for a parameter nside
    c=======================================================================
  */
  // 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   = (1.*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  = (int)fmod(ip,nl4) + 1;
    
    fodd  = 0.5 * (1 + fmod((double)(iring+nside),2));//  ! 1 if iring+nside is odd, 1/2 otherwise
    theta = acos( (nl2 - iring) / fact1 );
    phi   = (1.*iphi - fodd) * Pi /(2.*nside);
  }
  else {//! South Polar cap -----------------------------------
    
    ip    = npix - ipix1 + 1;
    hip   = ip/2.;
    fihip = 1.*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   = (1.*iphi - 0.5) * Pi/(2.*iring);
    //    cout << theta << " " << phi << endl;
  }
}

void SphereGorski::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
  /*
    c=======================================================================
    subroutine pix2ang_nest(nside, ipix, theta, phi)
    c=======================================================================
    c     gives theta and phi corresponding to pixel ipix (NESTED) 
    c     for a parameter nside
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)
    
      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*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + 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);

}


void SphereGorski::Pix2XY::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
  /*
    c=======================================================================
    subroutine pix2ang_nest(nside, ipix, theta, phi)
    c=======================================================================
    c     gives theta and phi corresponding to pixel ipix (NESTED) 
    c     for a parameter nside
    c=======================================================================
  */
  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
  //  (16/12/98)
    
      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*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + 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);

}