Index: trunk/SophyaLib/Samba/spherethetaphi.cc
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 276)
+++ trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 470)
@@ -1,7 +1,29 @@
-//
-// 
 #include "spherethetaphi.h"
 #include "nbmath.h"
+#include <complex>
+#include "piocmplx.h"
 #include <iostream.h>
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template SphereThetaPhi<double>
+#pragma define_template SphereThetaPhi<float>
+#pragma define_template SphereThetaPhi< complex<float> >
+#pragma define_template SphereThetaPhi< complex<double> >
+#pragma define_template FIO_SphereThetaPhi<double>
+#pragma define_template FIO_SphereThetaPhi<float>
+#pragma define_template FIO_SphereThetaPhi< complex<float> >
+#pragma define_template FIO_SphereThetaPhi< complex<double> >
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class SphereThetaPhi<double>;
+template class SphereThetaPhi<float>;
+template class SphereThetaPhi< complex<float> >;
+template class SphereThetaPhi< complex<double> >;
+template class FIO_SphereThetaPhi<double>;
+template class FIO_SphereThetaPhi<float>;
+template class FIO_SphereThetaPhi< complex<float> >;
+template class FIO_SphereThetaPhi< complex<double> >;
+#endif
+
 //***************************************************************
 //++ 
@@ -39,6 +61,6 @@
 //++
 
-
-SphereThetaPhi::SphereThetaPhi()
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi()
 
 //--
@@ -49,5 +71,6 @@
 /* --Methode-- */
 //++
-SphereThetaPhi::SphereThetaPhi(char* flnm)
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(char* flnm)
 
 //    Constructeur : charge une image à partir d'un fichier 
@@ -55,12 +78,14 @@
 {
   InitNul();
-  PInPersist s(flnm);
-  Read(s);
-}
-
-/* --Methode-- */
-
-//++
-SphereThetaPhi::SphereThetaPhi(int_4 m, int_4 pet)
+  cout << " ShereThetaPhi:: constructeur lecture sur fichier A FAIRE " << endl;
+  //PInPersist s(flnm);
+  //Read(s);
+}
+
+/* --Methode-- */
+
+//++
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(int_4 m)
 
 //    Constructeur : m est le nombre de bandes en theta sur un hémisphère 
@@ -71,7 +96,27 @@
 {
   InitNul();
-  Pixelize(m,pet);
-}
-
+  Pixelize(m);
+}
+
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
+{
+  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
+  NTheta_= s.NTheta_;
+  NPix_  = s.NPix_;
+  NPhi_  = new int_4[NTheta_];
+  Theta_ = new r_4[NTheta_+1];
+  TNphi_ = new int_4[NTheta_+1];
+  for(int k = 0; k < NTheta_; k++) 
+    {
+      NPhi_[k] = s.NPhi_[k];
+      Theta_[k]= s.Theta_[k];
+      TNphi_[k]= s.TNphi_[k];
+    }
+  Theta_[NTheta_]= s.Theta_[NTheta_];
+  TNphi_[NTheta_]= s.TNphi_[NTheta_];
+  Omega_ = s.Omega_;
+  pixels_= s.pixels_;
+}
 
 //++
@@ -79,5 +124,6 @@
 //--
 //++
-SphereThetaPhi::~SphereThetaPhi()
+template <class T>
+SphereThetaPhi<T>::~SphereThetaPhi()
 
 //--
@@ -89,147 +135,86 @@
 // Titre	Méthodes
 //--
-void SphereThetaPhi::InitNul()
+template <class T>
+void SphereThetaPhi<T>::InitNul()
 //
 //    initialise à zéro les variables de classe pointeurs
 {
-  mTheta = NULL;
-  mNphi = NULL;
-  mTNphi = NULL;
-  // mPix = NULL;
-  _pixel=NULL;
-  mNTheta = mNPet = 0;
-  mNPix = 0;
-}
-
-/* --Methode-- */
-void SphereThetaPhi::Clear()
-
-{
-  if (mTheta)  delete[] mTheta;
-  if (mNphi)  delete[] mNphi;
-  if (mTNphi) delete[] mTNphi;
-  //if (mPix)   delete[] mPix;
-  if (_pixel) _pixel->Detach();
-  mTheta = NULL;
-  mNphi = NULL;
-  mTNphi = NULL;
-  //mPix = NULL;
-  _pixel=NULL;
-  mNTheta = mNPet = 0;
-  mNPix = 0;
-}
-
-/* --Methode-- */
-//++
-void SphereThetaPhi::WriteSelf(POutPersist& s) const
-
-//    créer un fichier image
-//--
-{
-  char strg[256];
-  if (mInfo_) {sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo", (int_4)mNTheta, (int_4)mNPix);
-  }
-  else { sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d ", 
-               (int_4)mNTheta, (int_4)mNPix);
-  }
-  s.PutLine(strg);
-  if (mInfo_)  s << (*mInfo_);
-  s.PutI4(mNTheta);
-  s.PutI4(mNPet);
-  s.PutI4(mNPix);
-  s.PutR8(mOmeg);
-  s.PutR4s(mTheta, mNTheta+1);
-  s.PutI4s(mNphi, mNTheta);
-  s.PutI4s(mTNphi, mNTheta+1);
-  s.PutR8s(_pixel->Data(), mNPix);
-  //s.PutR8s(mPix, mNPix);
-  return;
-}
-
-/* --Methode-- */
-//++
-void SphereThetaPhi::ReadSelf(PInPersist& s)
-
-//    relit un fichier d'image
+  NTheta_= 0;
+  NPix_  = 0;
+  Theta_ = NULL;
+  NPhi_  = NULL;
+  TNphi_ = NULL;
+  pixels_.Reset();
+}
+
+/* --Methode-- */
+template <class T>
+void SphereThetaPhi<T>::Clear()
+
+{
+  if(Theta_)  delete[] Theta_;
+  if(NPhi_ )  delete[] NPhi_;
+  if(TNphi_)  delete[] TNphi_;
+  InitNul();
+}
+
+//++
+template <class T>
+void SphereThetaPhi<T>::Resize(int_4 m)
+//   re-pixelize the sphere
 //--
 {
   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(mNTheta);
-  s.GetI4(mNPet);
-  s.GetI4(mNPix);
-  s.GetR8(mOmeg);
-  mTheta = new r_4[mNTheta+1];
-  mNphi = new int_4[mNTheta];
-  mTNphi = new int_4[mNTheta+1];
-  //mPix = new r_8[mNPix+1];
-  _pixel=new PDataArray<r_8>(mNPix,true);
-  s.GetR4s(mTheta, mNTheta+1);
-  s.GetI4s(mNphi, mNTheta);
-  s.GetI4s(mTNphi, mNTheta+1);
-  s.GetR8s(_pixel->Data(), mNPix);
-  //s.GetR8s(mPix, mNPix);
-  return;
-}
-
-
-/* --Methode-- */
-//++
-int_4 SphereThetaPhi::NbPixels() const
+  Pixelize(m);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int_4 SphereThetaPhi<T>::NbPixels() const
 
 //    Retourne le nombre de pixels du découpage
 //--
 {
-  return(mNPix);
-}
-
-
-static r_8 dummy_pixel = 0;
-
-/* --Methode-- */
-//++
-r_8& SphereThetaPhi::PixVal(int_4 k)
+  return(NPix_);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+T& SphereThetaPhi<T>::PixVal(int_4 k)
 
 //    Retourne la valeur du contenu du pixel d'indice k
 //--
 {
-  if ( (k<0) || (k >= mNPix) ) {
-    //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
-    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
-  THROW(rangeCheckErr);
-
-    //throw "SphereThetaPhi::PIxVal Pixel index out of range";
-  }
-  return (*(_pixel->Data()+k));
-  //return(mPix[k]);
-}
-//++
-r_8 const& SphereThetaPhi::PixVal(int_4 k) const
+  if((k < 0) || (k >= NPix_)) 
+    {
+      //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
+      THROW(rangeCheckErr);
+    }
+  return pixels_(k);
+}
+
+//++
+template <class T>
+T const& SphereThetaPhi<T>::PixVal(int_4 k) const
 
 //    Retourne la valeur du contenu du pixel d'indice k
 //--
 {
-  if ( (k<0) || (k >= mNPix) ) {
-    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
-  //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
-
-    throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
-  } 
-  return *(_pixel->Data()+k);
-  //return(mPix[k]);
-}
-
-
-/* --Methode-- */
-//++
-int_4 SphereThetaPhi::PixIndexSph(r_4 theta, r_4 phi) const
+  if((k < 0) || (k >= NPix_)) 
+    {
+      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
+      //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+      throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
+    } 
+  return *(pixels_.Data()+k);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int_4 SphereThetaPhi<T>::PixIndexSph(r_4 theta, r_4 phi) const
 
 //    Retourne l'indice du pixel vers lequel pointe une direction définie par 
@@ -246,20 +231,19 @@
   
   // La bande d'indice kt est limitée par les valeurs de theta contenues dans
-  // mTheta[kt] et mTheta[kt+1]
-  for( i=1; i< mNTheta; i++ ) 
-    if( theta < mTheta[i] ) break;
-  
-  dphi= (r_4)DeuxPi/(r_4)mNphi[i-1];
-  
-  if (fgzn)  k= mNPix-mTNphi[i]+(int_4)(phi/dphi);
-  else k= mTNphi[i-1]+(int_4)(phi/dphi); 
-  
+  // Theta_[kt] et Theta_[kt+1]
+  for( i=1; i< NTheta_; i++ ) 
+    if( theta < Theta_[i] ) break;
+  
+  dphi= (r_4)DeuxPi/(r_4)NPhi_[i-1];
+  
+  if (fgzn)  k= NPix_-TNphi_[i]+(int_4)(phi/dphi);
+  else k= TNphi_[i-1]+(int_4)(phi/dphi); 
   return(k);
 }
 
-
-/* --Methode-- */
-//++
-void SphereThetaPhi::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const
 
 //    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
@@ -269,26 +253,25 @@
   bool fgzn = false;
   
-  if( (k < 0) || (k >= mNPix) ) {theta = -99999.; phi = -99999.;  return; }
-  if( k >= mNPix/2)  {fgzn = true;  k = mNPix-1-k; }
+  if( (k < 0) || (k >= NPix_) ) {theta = -99999.; phi = -99999.;  return; }
+  if( k >= NPix_/2)  {fgzn = true;  k = NPix_-1-k; }
 
   // recupère l'indice i de la tranche qui contient le pixel k
-  for( i=0; i< mNTheta; i++ ) 
-    if( k < mTNphi[i+1] ) break;
+  for( i=0; i< NTheta_; i++ ) 
+    if( k < TNphi_[i+1] ) break;
 
   // angle theta
-  theta= 0.5*(mTheta[i]+mTheta[i+1]);
+  theta= 0.5*(Theta_[i]+Theta_[i+1]);
   if (fgzn) theta= Pi-theta;
   
   // angle phi
-  k -= mTNphi[i];
-  phi= (r_4)DeuxPi/(r_4)mNphi[i]*(r_4)(k+.5);
+  k -= TNphi_[i];
+  phi= (r_4)DeuxPi/(r_4)NPhi_[i]*(r_4)(k+.5);
   if (fgzn) phi= DeuxPi-phi;
-
-  return;
-}
-
-
-//++
-r_8         SphereThetaPhi::PixSolAngle(int_4 dummy) const
+}
+
+//++
+template <class T>
+r_8  SphereThetaPhi<T>::PixSolAngle(int_4 dummy) const
+
 //    Pixel Solid angle  (steradians)
 //    All the pixels have the same solid angle. The dummy argument is
@@ -297,12 +280,11 @@
 //--
 {
-  return mOmeg;
-}
-
-
-/* --Methode-- */
-//++
-void SphereThetaPhi::Limits(int_4 k, r_4& tetMin, r_4& tetMax, 
-			   r_4& phiMin, r_4& phiMax) 
+  return Omega_;
+}
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::Limits(int_4 k, r_4& tetMin, r_4& tetMax, r_4& phiMin, r_4& phiMax) 
 
 //    Retourne les valeurs de theta et phi limitant le pixel d'indice k
@@ -313,5 +295,5 @@
   bool fgzn= false;
   
-  if( (k< 0) || (k >= mNPix) ) {
+  if( (k< 0) || (k >= NPix_) ) {
     tetMin= -99999.; 
     phiMin= -99999.;  
@@ -322,29 +304,29 @@
   
   // si on se trouve dans l'hémisphère Sud
-  if( k >= mNPix/2 ) { 
+  if( k >= NPix_/2 ) { 
     fgzn= true; 
-    k= mNPix-1-k; 
+    k= NPix_-1-k; 
   }
   
   // on recupere l'indice i de la tranche qui contient le pixel k
   int i;
-  for( i=0; i< mNTheta; i++ ) 
-    if( k< mTNphi[i+1] ) break;
+  for( i=0; i< NTheta_; i++ ) 
+    if( k< TNphi_[i+1] ) break;
   
   // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= mTheta[i];
-  tetMax= mTheta[i+1];
+  tetMin= Theta_[i];
+  tetMax= Theta_[i+1];
   // valeurs limites de theta dans l'hemisphere Sud
   if (fgzn) {
-    tetMin= Pi-mTheta[i+1];
-    tetMax= Pi-mTheta[i];
+    tetMin= Pi-Theta_[i+1];
+    tetMax= Pi-Theta_[i];
   }
   
   // pixel correspondant dans l'hemisphere Nord
-  if (fgzn) k= mTNphi[i+1]-k+mTNphi[i]-1;
+  if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
   
   // indice j de discretisation ( phi= j*dphi )
-  j= k-mTNphi[i];
-  dphi= (r_4)DeuxPi/(r_4)mNphi[i];
+  j= k-TNphi_[i];
+  dphi= (r_4)DeuxPi/(r_4)NPhi_[i];
   
   // valeurs limites de phi 
@@ -356,5 +338,6 @@
 /* --Methode-- */
 //++
-int_4 SphereThetaPhi::NbThetaSlices() const
+template <class T>
+int_4 SphereThetaPhi<T>::NbThetaSlices() const
 
 //    Retourne le nombre de tranches en theta sur la sphere
@@ -362,5 +345,5 @@
 {
   int_4 nbslices;
-  nbslices= 2*mNTheta;
+  nbslices= 2*NTheta_;
   return(nbslices);
 }
@@ -368,22 +351,21 @@
 /* --Methode-- */
 //++
-int_4 SphereThetaPhi::NPhi(int_4 kt) const
+template <class T>
+int_4 SphereThetaPhi<T>::NPhi(int_4 kt) const
 
 //    Retourne le nombre de pixels en phi de la tranche kt
 //-- 
 {
-  
-  int_4 nbpix;
-  
+  int_4 nbpix;  
   // verification
-  if( (kt< 0) || (kt>= 2*mNTheta) ) return(-1);
+  if( (kt< 0) || (kt>= 2*NTheta_) ) return(-1);
   
   // si on se trouve dans l'hemisphere Sud
-  if( kt >= mNTheta ) {
-    kt= 2*mNTheta-1-kt; 
+  if( kt >= NTheta_ ) {
+    kt= 2*NTheta_-1-kt; 
   }
   
   // nombre de pixels
-  nbpix= mNphi[kt];
+  nbpix= NPhi_[kt];
   return(nbpix);
 }
@@ -392,14 +374,13 @@
 /* --Methode-- */
 //++
-void SphereThetaPhi::Theta(int_4 kt,r_4& tetMin,r_4& tetMax) 
+template <class T>
+void SphereThetaPhi<T>::Theta(int_4 kt,r_4& tetMin,r_4& tetMax) 
 
 //    Retourne les valeurs de theta limitant la tranche kt
 //--
 {
-  
   bool fgzn= false;
-  
   // verification
-  if( (kt< 0) || (kt>= 2*mNTheta) ) {
+  if( (kt< 0) || (kt>= 2*NTheta_) ) {
     tetMin= -99999.; 
     tetMax= -99999.; 
@@ -408,16 +389,16 @@
 
   // si on se trouve dans l'hemisphere Sud
-  if( kt >= mNTheta ) {
+  if( kt >= NTheta_ ) {
     fgzn= true;
-    kt= 2*mNTheta-1-kt; 
+    kt= 2*NTheta_-1-kt; 
   }
   
   // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= mTheta[kt];
-  tetMax= mTheta[kt+1];
+  tetMin= Theta_[kt];
+  tetMax= Theta_[kt+1];
   // valeurs limites de theta dans l'hemisphere Sud
   if (fgzn) {
-    tetMin= Pi-mTheta[kt+1];
-    tetMax= Pi-mTheta[kt];
+    tetMin= Pi-Theta_[kt+1];
+    tetMax= Pi-Theta_[kt];
   }  
 }
@@ -425,14 +406,12 @@
 /* --Methode-- */
 //++
-void SphereThetaPhi::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax) 
+template <class T>
+void SphereThetaPhi<T>::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax) 
 
 //    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
 //--
 {
-  
-  r_4 dphi;
-  
   // verification
-  if( (kt< 0) || (kt>= 2*mNTheta) ) {
+  if( (kt< 0) || (kt>= 2*NTheta_) ) {
     phiMin= -99999.; 
     phiMax= -99999.; 
@@ -441,8 +420,8 @@
   
   // si on se trouve dans l'hemisphere Sud
-  if( kt >= mNTheta ) kt= 2*mNTheta-1-kt;
+  if( kt >= NTheta_ ) kt= 2*NTheta_-1-kt;
   
   // verifie si la tranche kt contient au moins jp pixels
-  if( (jp< 0) || (jp >= mNphi[kt]) ) {
+  if( (jp< 0) || (jp >= NPhi_[kt]) ) {
     phiMin= -88888.; 
     phiMax= -88888.; 
@@ -450,5 +429,5 @@
   }
   
-  dphi= (r_4)DeuxPi/(r_4)mNphi[kt];
+  r_4 dphi= (r_4)DeuxPi/(r_4)NPhi_[kt];
   phiMin= dphi*(r_4)(jp);
   phiMax= dphi*(r_4)(jp+1);
@@ -458,30 +437,32 @@
 /* --Methode-- */
 //++
-int_4 SphereThetaPhi::Index(int_4 kt,int_4 jp) const
+template <class T>
+int_4 SphereThetaPhi<T>::Index(int_4 kt,int_4 jp) const
 
 //    Retourne l'indice du pixel d'indice jp dans la tranche kt 
 //--
 {
-  
   int_4 k;
   bool fgzn= false;
   
   // si on se trouve dans l'hemisphere Sud
-  if( kt >= mNTheta ) {
+  if( kt >= NTheta_ ) {
     fgzn= true; 
-    kt= 2*mNTheta-1-kt;
+    kt= 2*NTheta_-1-kt;
   }
   
   // si la tranche kt contient au moins i pixels
-  if( (jp>=0) && (jp<mNphi[kt]) ) {
-    // dans l'hemisphere Sud
-    if (fgzn) k= mNPix-mTNphi[kt+1]+jp;
-    // dans l'hemisphere Nord
-    else k= mTNphi[kt]+jp;
-  }
-  else{
-    k= 9999;
-    printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
-  }
+  if( (jp>=0) && (jp<NPhi_[kt]) ) 
+    {
+      // dans l'hemisphere Sud
+      if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
+      // dans l'hemisphere Nord
+      else k= TNphi_[kt]+jp;
+    }
+  else
+    {
+      k= 9999;
+      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
+    }
   return(k);
 }
@@ -489,34 +470,34 @@
 /* --Methode-- */
 //++
-void SphereThetaPhi::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 
+template <class T>
+void SphereThetaPhi<T>::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 
 
 //    Retourne les indices kt et jp du pixel d'indice k
 //--
 {
-  
-  bool fgzn= false;
-  
+  bool fgzn= false;  
   // si on se trouve dans l'hemisphere Sud
-  if( k >= mNPix/2 ) { 
+  if( k >= NPix_/2 ) { 
     fgzn= true; 
-    k= mNPix-1-k; 
+    k= NPix_-1-k; 
   }
   
   // on recupere l'indice kt de la tranche qui contient le pixel k
   int i;
-  for( i=0; i< mNTheta; i++ ) 
-    if( k< mTNphi[i+1] ) break;
+  for( i=0; i< NTheta_; i++ ) 
+    if( k< TNphi_[i+1] ) break;
   
   // indice  kt de tranche
-  if (fgzn) kt= 2*mNTheta-1-i;
+  if (fgzn) kt= 2*NTheta_-1-i;
   else kt= i;
   
   // indice jp de pixel
-  if (fgzn) jp= mTNphi[i+1]-k-1;
-  else jp= k-mTNphi[i];
-  
-}
-//++
-void  SphereThetaPhi::Pixelize( int_4 m, int_4 pet) 
+  if (fgzn) jp= TNphi_[i+1]-k-1;
+  else jp= k-TNphi_[i];
+  
+}
+//++
+template <class T>
+void  SphereThetaPhi<T>::Pixelize( int_4 m) 
 
 //    effectue le découpage en pixels (m et pet ont la même signification 
@@ -531,6 +512,5 @@
 {
   int_4 ntotpix,i,j;
-  //r_8 x, omeg, omeg2pi, teti, tetm, tet, nphi, costeti;
-  r_8 omeg2pi;
+  
   // Decodage et controle des arguments d'appel :
   // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
@@ -538,22 +518,18 @@
   if (m > 16384) m = 16384;
   
-  // Au moins 4 et au plus 256 pixels dans la premiere bande decoupee en phi
-  if (pet < 4) pet = 4;
-  if (pet > 256) pet = 256;
-  
   // On memorise les arguments d'appel
-  mNTheta = m;  mNPet = pet;
+  NTheta_ = m;  
   
   // On commence par decouper l'hemisphere z>0.
   // Creation des vecteurs contenant :
   // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
-  mTheta = new r_4[m+1];
+  Theta_= new r_4[m+1];
   
   // Le nombre de pixels en phi de chacune des bandes en theta
-  mNphi = new int_4[m];
+  NPhi_ = new int_4[m];
   
   // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
   // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
-  mTNphi = new int_4[m+1];
+  TNphi_= new int_4[m+1];
   
   // Calcul du nombre total de pixels dans chaque bande pour optimiser
@@ -561,43 +537,44 @@
   
   //calotte polaire
-  mTNphi[0]=0;
-  mNphi[0] = 1;
+  TNphi_[0]= 0;
+  NPhi_[0] = 1;
   
   //bandes jusqu'a l'equateur
-  for(j=1; j < m; j++) {
-    mTNphi[j] = mTNphi[j-1]+mNphi[j-1];
-    mNphi[j]  = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
-  }
+  for(j = 1; j < m; j++) 
+    {
+      TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
+      NPhi_[j] = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
+    }
   
   // Nombre total de pixels sur l'hemisphere
-  ntotpix = mTNphi[m-1]+mNphi[m-1];
-  mTNphi[m] = ntotpix;
+  ntotpix  = TNphi_[m-1]+NPhi_[m-1];
+  TNphi_[m]= ntotpix;
   // et sur la sphere entiere
-  mNPix=2*ntotpix;
+  NPix_= 2*ntotpix;
   
   // Creation et initialisation du vecteur des contenus des pixels 
-  //mPix = new r_8[mNPix];
-  //for(i=0; i<mNPix; i++)  mPix[i] = 0.;
-  _pixel=new PDataArray<r_8>(mNPix,true);
-  for(i=0; i<mNPix; i++)  *(_pixel->Data()+i)=0.;
+  pixels_.ReSize(NPix_);
+  pixels_.Reset();
+
   // Determination des limites des bandes en theta :
   // omeg est l'angle solide couvert par chaque pixel,
-  // une bande donnee kt couvre un angle solide mNphi[kt]*omeg
-  // egal a 2* Pi*(cos mTheta[kt]-cos mTheta[kt+1]). De meme, l'angle solide 
+  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
+  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide 
   //de la
   // calotte allant du pole a la limite haute de la bande kt vaut
-  // 2* Pi*(1.-cos mTheta[kt+1])= mTNphi[kt]*omeg...
-  
-  omeg2pi = 1./(r_8)ntotpix;
-  mOmeg = omeg2pi*DeuxPi;
-  
-  for(j=0; j <= m; j++) {
-    mTheta[j] = acos(1.-(double)mTNphi[j]*omeg2pi);
-  }
-}
-
-
-//++
-void  SphereThetaPhi::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const 
+  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
+  
+  double omeg2pi= 1./(double)ntotpix;
+  Omega_ = omeg2pi*DeuxPi;
+  
+  for(j=0; j <= m; j++) 
+    {
+      Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
+    }
+}
+
+//++
+template <class T>
+void  SphereThetaPhi<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 
@@ -609,27 +586,241 @@
 {
   cout << "entree GetThetaSlice, couche no " << index << endl;
-  if (index<0 || index > NbThetaSlices()) {
-    // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
-    cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
-  THROW(rangeCheckErr);
-  }
-  int_4 iring=Index(index,0);
-  int_4 bid=this->NPhi(index);
-  int lring=bid;
+
+  if(index < 0 || index > NbThetaSlices()) 
+    {
+      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
+      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
+      THROW(rangeCheckErr);
+    }
+
+  int_4 iring= Index(index,0);
+  int_4 bid  = this->NPhi(index);
+  int lring  = bid;
   cout << " iring= " << iring << " lring= " << lring << endl;
-  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-- */
-
+
+  phi.ReSize(lring);
+  value.ReSize(lring);
+  float Te= 0.;
+  float Fi= 0.;
+  for(int kk = 0; kk < lring; kk++) 
+    {
+      PixThetaPhi(kk+iring,Te,Fi);
+      phi(kk)= Fi;
+      value(kk)= PixVal(kk+iring);
+    }
+  theta= Te;
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmNPhi(int_4* array, int_4 m)
+  //remplit  le tableau contenant le nombre de pixels en phi de chacune des bandes en theta 
+  //--
+{
+  NPhi_= new int_4[m];
+  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmTNphi(int_4* array, int_4 m)
+  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes 
+  //--
+{
+  TNphi_= new int_4[m];
+  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmTheta(r_4* array, int_4 m)
+  //remplit  le tableau contenant les valeurs limites de theta
+  //--
+{
+  Theta_= new r_4[m];
+  for(int k = 0; k < m; k++) Theta_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setDataBlock(T* data, int_4 m)
+  // remplit  le vecteur des contenus des pixels 
+{
+  pixels_.FillFrom(m,data);
+} 
+
+template <class T>
+void SphereThetaPhi<T>::print(ostream& os) const
+{
+  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
+  //
+  os << " NTheta_= " << NTheta_ << endl;
+  os << " NPix_    = " << NPix_   << endl;
+  os << " Omega_  =  " << Omega_   << endl;
+
+  os << " contenu de NPhi_ : ";
+  for(int i=0; i < NTheta_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << NPhi_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de Theta_ : ";
+  for(int i=0; i < NTheta_+1; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << Theta_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de TNphi_ : ";
+  for(int i=0; i < NTheta_+1; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << TNphi_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de pixels : ";
+  for(int i=0; i < NPix_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os <<  pixels_(i) <<", ";
+    }
+  os << endl;
+}
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+//////////////////////////////////////////////////////////
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
+{
+  dobj= new SphereThetaPhi<T>;
+  ownobj= true;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
+{
+  dobj= new SphereThetaPhi<T>;
+  ownobj= true;
+  Read(filename);
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
+{
+  dobj= new SphereThetaPhi<T>(obj);
+  ownobj= true;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
+{
+  dobj= obj;
+  ownobj= false;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
+{
+  if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
+{
+  return(dobj);
+}
+
+template <class T>
+void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
+{
+  cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      dobj= new SphereThetaPhi<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 mNTheta;
+  is.GetI4(mNTheta);  
+  dobj->setSizeIndex(mNTheta);
+
+  int_4 mNPix;
+  is.GetI4(mNPix);
+  dobj->setNbPixels(mNPix);
+
+  r_8 mOmeg;
+  is.GetR8(mOmeg);
+  dobj->setPixSolAngle(mOmeg);
+
+  int_4* mNphi= new int_4[mNTheta];
+  is.GetI4s(mNphi, mNTheta);
+  dobj->setmNPhi(mNphi, mNTheta);
+  delete [] mNphi;
+
+  int_4* mTNphi= new int_4[mNTheta+1];
+  is.GetI4s(mTNphi, mNTheta+1);
+  dobj->setmTNphi(mTNphi, mNTheta+1);
+  delete [] mTNphi;
+
+  r_4* mTheta= new r_4[mNTheta+1];
+  is.GetR4s(mTheta, mNTheta+1);
+  dobj->setmTheta(mTheta, mNTheta+1);
+  delete [] mTheta;
+
+  T* mPixels= new T[mNPix];
+  //is.Get(mPixels, mNPix);
+  PIOSReadArray(is, mPixels, mNPix);
+  dobj->setDataBlock(mPixels, mNPix);
+  delete [] mPixels;
+}
+
+template <class T>
+void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
+{
+  cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      cout << " WriteSelf:: dobj= null " << endl;
+      return;
+    }
+
+  char strg[256];
+  int_4 mNTheta= dobj->SizeIndex();
+  int_4 mNPix  = dobj->NbPixels();
+ 
+  if(dobj->ptrInfo()) 
+    {
+      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
+      os.PutLine(strg);
+      os << dobj->Info();
+    }
+  else 
+    { 
+      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
+      os.PutLine(strg);  
+    }
+
+  os.PutI4(mNTheta);  
+  os.PutI4(mNPix);
+  os.PutR8(dobj->PixSolAngle(0));
+  os.PutI4s(dobj->getmNPhi() , mNTheta);
+  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
+  os.PutR4s(dobj->getmTheta(), mNTheta+1);
+  //os.Put((dobj->getDataBlock())->Data(), mNPix);
+  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
+}
