Index: trunk/SophyaLib/Samba/spheregorski.cc
===================================================================
--- trunk/SophyaLib/Samba/spheregorski.cc	(revision 276)
+++ trunk/SophyaLib/Samba/spheregorski.cc	(revision 470)
@@ -2,5 +2,9 @@
 #include "spheregorski.h"
 #include "strutil.h"
-extern "C" {
+#include <complex>
+#include "piocmplx.h"
+
+extern "C" 
+{
 #include <stdio.h>
 #include <stdlib.h>
@@ -8,19 +12,119 @@
 }
       
-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&);
-
-}
+extern "C" 
+{
+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*);
+}
+
+//*******************************************************************
+// Class PIXELS_XY 
+//  Construction des tableaux necessaires a la traduction des indices RING en 
+//  indices NESTED (ou l'inverse)
+//*******************************************************************
+
+PIXELS_XY::PIXELS_XY()
+{
+  cout << " appel du constructeur PIXELS_XY " <<endl;
+  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;
+	}
+    }     
+}
+
 //*******************************************************************
 //++ 
@@ -81,25 +185,16 @@
 //++
 
-SphereGorski::SphereGorski()
-
-//--
-{
+template<class T>
+SphereGorski<T>::SphereGorski()
+
+//--
+{
+  cout<<" appel du constructeur SphereGorski ()" <<endl;
   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)
+//++
+template<class T>
+SphereGorski<T>::SphereGorski(int_4 m)
 
 //    Constructeur : m est la variable nside de l'algorithme de Gorski
@@ -108,55 +203,87 @@
 //--
 {
-  //printf(" initialisation par defaut SphereGorski \n");
+  cout<<" appel du constructeur SphereGorski (m)" <<endl;
+
+  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);
+}
+
+template<class T>
+SphereGorski<T>::SphereGorski(const SphereGorski<T>& s)
+{
+  cout << " constructeur de recopie " << endl;
+  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
+
+  nSide_= s.nSide_;
+  nPix_ = s.nPix_;
+  omeg_ = s.omeg_;
+
+  pixels_= s.pixels_;
+
+  nlmax_= s.nlmax_;
+  nmmax_= s.nmmax_;
+  iseed_= s.iseed_;
+  fwhm_ = s.fwhm_;
+  quadrupole_ = s.quadrupole_;
+  sym_cut_deg_= s.sym_cut_deg_;
+  strcpy(powFile_,s.powFile_);
+}
+
+//++
+// Titre	Destructeur
+//--
+//++
+template<class T>
+SphereGorski<T>::~SphereGorski()
+
+//--
+{
+  InitNul();
+}
+
+//++
+// Titre	Méthodes
+//--
+
+//++
+template<class T>
+void SphereGorski<T>::Resize(int_4 m)
+
+//    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)
+//--
+{
   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 
+  // 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);
-  }
-
+  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) 
+}
+
+template<class T>
+void  SphereGorski<T>::Pixelize( int_4 m) 
 
 //    prépare la pixelisation Gorski (m a la même signification 
@@ -166,205 +293,67 @@
 //--
 {
-  
   // 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
+  // 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 
-  mPix_ = new r_8[nPix_];
-  for(int i=0; i<nPix_; i++)  mPix_[i] = 0.;
+  pixels_.ReSize(nPix_);
+  pixels_.Reset();
 
   // solid angle per pixel   
-  omeg_=4*Pi/nPix_;
-}
-
-void SphereGorski::InitNul()
+  omeg_= 4*Pi/nPix_;
+}
+
+template<class T>
+void SphereGorski<T>::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;
-}
+  nSide_= 0;
+  nPix_ = 0;
+  omeg_ = 0.;
+  pixels_.Reset();
+
+  nlmax_= 0;
+  nmmax_= 0;
+  iseed_= 0;
+  fwhm_ = 0.;
+  quadrupole_ = 0.;
+  sym_cut_deg_= 0.;
+  for(int k = 0; k < 128; k++) powFile_[k]=' ';
+}
+
 /* --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_)  s << (*mInfo_);
-  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;
-    s >> (*mInfo_);
-  }
-  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
+//++
+template<class T>
+int_4 SphereGorski<T>::NbPixels() const
 
 //    Retourne le nombre de pixels du découpage
 //--
-{
+{  
   return(nPix_);
 }
 
 //++
-int_4 SphereGorski::NbThetaSlices() const 
+template<class T>
+int_4 SphereGorski<T>::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 
+{
+  return int_4(4*nSide_-1);
+}
+
+//++
+template<class T>
+void  SphereGorski<T>::GetThetaSlice(int_4 index, r_4& theta, TVector<float>& phi, TVector<T>& value) const 
 
 //    Retourne, pour la tranche en theta d'indice 'index' le theta 
@@ -374,106 +363,116 @@
 //--
 {
-
-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;
-
-}
-
-
+  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.ReSize(lring);
+  value.ReSize(lring);
+  float TH=0.;
+  float F =0.;
+  for(int kk = 0; kk < lring;kk++) 
+    {
+      PixThetaPhi(kk+iring,TH,F);
+      phi(kk)= F;
+      value(kk)= PixVal(kk+iring);
+    }
+  theta= TH;
+}
 
 /* --Methode-- */
 //++
-r_8& SphereGorski::PixVal(int_4 k)
+template<class T>
+T& SphereGorski<T>::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]);
+  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 pixels_(k);
 }
 
 /* --Methode-- */
 //++
-r_8 const& SphereGorski::PixVal(int_4 k) const
+template<class T>
+T const& SphereGorski<T>::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)
+  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 *(pixels_.Data()+k);
+}
+
+//++
+template<class T>
+T& SphereGorski<T>::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
+  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
 
 //    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];
-}
-
-
+  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-- */
 //++
-int_4 SphereGorski::PixIndexSph(r_4 theta, r_4 phi) const
+template<class T>
+int_4 SphereGorski<T>::PixIndexSph(r_4 theta, r_4 phi) const
 
 //    Retourne l'indice "RING" du pixel vers lequel pointe une direction 
@@ -481,9 +480,10 @@
 //--
 {
-  return ang2pix_ring(nSide_, double(theta), double(phi));
-}
-
-//++
-int_4 SphereGorski::PixIndexSphNest(r_4 theta, r_4 phi) const
+  return ang2pix_ring(nSide_,double(theta),double(phi));
+}
+
+//++
+template<class T>
+int_4 SphereGorski<T>::PixIndexSphNest(r_4 theta, r_4 phi) const
 
 //    Retourne l'indice NESTED" du pixel vers lequel pointe une direction 
@@ -491,5 +491,5 @@
 //--
 {
-  return ang2pix_nest(nSide_, double(theta), double(phi));
+  return ang2pix_nest(nSide_,double(theta),double(phi));
 }
 
@@ -497,5 +497,6 @@
 /* --Methode-- */
 //++
-void SphereGorski::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const
+template<class T>
+void SphereGorski<T>::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
@@ -504,11 +505,12 @@
   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
+  pix2ang_ring(nSide_,k, t, p);
+  teta= (r_4)t;
+  phi = (r_4)p;
+}
+
+//++
+template<class T>
+r_8 SphereGorski<T>::PixSolAngle(int_4 dummy) const
 //    Pixel Solid angle  (steradians)
 //    All the pixels have the same solid angle. The dummy argument is
@@ -520,7 +522,7 @@
 }
 
-
-//++
-void SphereGorski::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const
+//++
+template<class T>
+void SphereGorski<T>::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const
 
 //    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice 
@@ -531,10 +533,11 @@
   double p;
   pix2ang_nest(nSide_, k, t, p);
-  teta=(r_4)t;
-  phi=(r_4)p;
-}
-
-//++
-int_4 SphereGorski::NestToRing(int_4 k) 
+  teta= (r_4)t;
+  phi = (r_4)p;
+}
+
+//++
+template<class T>
+int_4 SphereGorski<T>::NestToRing(int_4 k) const
 
 //    conversion d'index NESTD en un index RING
@@ -544,6 +547,8 @@
   return  nest2ring(nSide_,k);
 }
-//++
-int_4 SphereGorski::RingToNest(int_4 k) 
+
+//++
+template<class T>
+int_4 SphereGorski<T>::RingToNest(int_4 k) const
 //
 //    conversion d'index RING en un index NESTED
@@ -553,6 +558,9 @@
   return  ring2nest(nSide_,k);
 }
-//++
-void SphereGorski::anharm(int nlmax, float sym_c,float* powspec) 
+
+/*
+//++
+template<class T> 
+void SphereGorski<T>::anharm(int nlmax, float sym_c,float* powspec) 
 //
 //    analyse en harmoniques spheriques des valeurs des pixels de la 
@@ -573,8 +581,8 @@
 //
 { 
-  if (nlmax > 2*nSide_) {
-    cout << " anharm : nlmax= " << nlmax << 
-            " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
-    exit(1);
+    if (nlmax > 2*nSide_) {
+          cout << " anharm : nlmax= " << nlmax << 
+                          " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
+         exit(1);
   }
   else {
@@ -584,5 +592,5 @@
   sym_cut_deg_=sym_c;
   float* map=new float[nPix_];
-  for (int k=0; k<nPix_; k++) map[k]=(float)mPix_[k];
+  for (int k=0; k<nPix_; k++) map[k]=(float)pixels_(k);
   int nsmax=nSide_;
   int nmmax=nmmax_;
@@ -607,11 +615,12 @@
   delete [] phas_s;
   delete [] dataw;
-  delete [] work;
-  
-}
-
-
-//++
-void SphereGorski::synharm(int nlmax, int iseed,float fwhm, float* powspec) 
+  delete [] work;  
+}
+*/
+
+/*
+//++
+template<class T>
+void SphereGorski<T>::synharm(int nlmax, int iseed,float fwhm, float* powspec) 
 //
 //    synthese  des valeurs des pixels de la sphere par l'intermediaire 
@@ -654,8 +663,7 @@
   int nsmax=nSide_;
   int nmmax=nmmax_;
-  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
-  
-
-//    tableaux de travail 
+  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)];
@@ -666,5 +674,5 @@
   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];
+  for (int k=0; k<nPix_; k++) pixels_(k) = (T)map[k];
   delete [] map;
   delete [] alm_T;
@@ -675,137 +683,102 @@
   delete [] work;
   delete [] lread;
-
-}
-
-int SphereGorski::nest2ring(int nside, int ipnest) const {
+}
+*/
+
+template<class T>
+int SphereGorski<T>::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)
+  
+  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 SphereGorski<T>::ring2nest(int nside, int 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;
@@ -898,58 +871,18 @@
   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=(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}
 }
 
-void SphereGorski::mk_xy2pix(int *x2pix, int *y2pix) {
+template<class T>
+int SphereGorski<T>::ang2pix_ring(int nside, double theta, double phi) const 
+{
   /*
-    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=======================================================================
-  */
+    c==================================================
+    */
   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
   //  (16/12/98)
@@ -1017,12 +950,13 @@
   }
     return (ipix1 - 1);// ! in {0, npix-1}
-  
-}
-
-int SphereGorski::ang2pix_nest(int nside, double theta, double phi) const {
+}
+
+template<class T>
+int SphereGorski<T>::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
@@ -1033,89 +967,91 @@
     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 {
+  
+  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}
+}
+
+template<class T>
+void SphereGorski<T>::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)
@@ -1178,183 +1114,338 @@
 }
 
-void SphereGorski::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
+template<class T>
+void SphereGorski<T>::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;
+
+  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);
+}
+
+// retourne le nom du fichier qui contient le spectre de puissance
+template<class T>
+void SphereGorski<T>::powfile(char filename[]) const 
+{
+  bool status = false;
+  for (int k=0; k< 128; k++)
+    {
+      if( powFile_[k] != ' ' )
+	{
+	  status = true;
+	  break;
 	}
-      }
-      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);
-
-}
-
-
-
+    }
+  if ( status ) 
+    {
+      strcpy(filename,powFile_);
+    }
+  else
+    {
+      strcpy(filename,"no file");
+    }
+}
+
+template<class T>
+void SphereGorski<T>::getParafast(int_4& nlmax,int_4& nmmax,int_4& iseed,float& fwhm,float& quadr,float& cut) const
+{
+  nlmax= nlmax_;
+  nmmax= nmmax_;
+  iseed= iseed_;
+  fwhm = fwhm_;
+  quadr= quadrupole_;
+  cut  = sym_cut_deg_;
+}
+
+template<class T>
+void SphereGorski<T>::setParafast(int_4 nlmax,int_4 nmmax,int_4 iseed,float fwhm,float quadr,float cut,char* filename)
+{
+  nlmax_= nlmax;
+  nmmax_= nmmax;
+  iseed_= iseed;
+  fwhm_ = fwhm;
+  quadrupole_ = quadr;
+  sym_cut_deg_= cut; 
+  strcpy(powFile_,filename);
+}
+
+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 << " contenu de 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;
+
+  os << " les parametres des modules anafast et synfast " <<endl;
+  os << " nlmax, nmmax & iseed= " <<nlmax_<<", "<<nmmax_<<", "<<iseed_<<endl;
+  os << " fwhm, quadr & cut = "<<fwhm_<<", "<<quadrupole_<<", "<<sym_cut_deg_<<endl;
+  os << " powfile= " << powFile_<<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>;
+  ownobj= true;
+  Read(filename);
+}
+
+template <class T>
+FIO_SphereGorski<T>::FIO_SphereGorski(const SphereGorski<T>& obj)
+{
+  dobj= new SphereGorski<T>(obj);
+  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>::ReadSelf(PInPersist& is)
+{
+  cout << " FIO_SphereGorski:: ReadSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      dobj= new SphereGorski<T>;
+    }
+
+  // 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);
+  dobj->setSizeIndex(nSide);
+
+  int_4 nPix;
+  is.GetI4(nPix);
+  dobj->setNbPixels(nPix);
+
+  r_8 Omega;
+  is.GetR8(Omega);
+  dobj->setPixSolAngle(Omega);
+
+  T* pixels= new T[nPix];
+  PIOSReadArray(is, pixels, nPix);
+  dobj->setDataBlock(pixels, nPix);
+  delete [] pixels;
+
+  int_4 nlmax,nmmax,iseed;
+  is.GetI4(nlmax);
+  is.GetI4(nmmax);
+  is.GetI4(iseed);
+  
+  float fwhm,quadr,cut;
+  is.GetR4(fwhm);
+  is.GetR4(quadr);
+  is.GetR4(cut);
+
+  char powfl[128];
+  is.GetLine(powfl, 127);
+
+  dobj->setParafast(nlmax,nmmax,iseed,fwhm,quadr,cut,powfl);
+}
+
+template <class T>
+void FIO_SphereGorski<T>::WriteSelf(POutPersist& os) const
+{
+  cout << " FIO_SphereGorski:: WriteSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      cout << " WriteSelf:: dobj= null " << endl;
+      return;
+    }
+
+  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(0));
+
+  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
+
+  int_4 nlmax,nmmax,iseed;
+  float fwhm,quadr,cut;
+  dobj->getParafast(nlmax,nmmax,iseed,fwhm,quadr,cut);
+  os.PutI4(nlmax);
+  os.PutI4(nmmax);
+  os.PutI4(iseed);
+  os.PutR4(fwhm);
+  os.PutR4(quadr);
+  os.PutR4(cut);
+
+  char powfl[128];
+  dobj->powfile(powfl);
+  os.PutLine(powfl);
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template SphereGorski<double>
+#pragma define_template SphereGorski<float>
+#pragma define_template SphereGorski< complex<float> >
+#pragma define_template SphereGorski< complex<double> >
+#pragma define_template FIO_SphereGorski<double>
+#pragma define_template FIO_SphereGorski<float>
+#pragma define_template FIO_SphereGorski< complex<float> >
+#pragma define_template FIO_SphereGorski< complex<double> >
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class SphereGorski<double>;
+template class SphereGorski<float>;
+template class SphereGorski< complex<float> >;
+template class SphereGorski< complex<double> >;
+template class FIO_SphereGorski<double>;
+template class FIO_SphereGorski<float>;
+template class FIO_SphereGorski< complex<float> >;
+template class FIO_SphereGorski< complex<double> >;
+#endif
+
