Changeset 470 in Sophya

Timestamp:

Oct 15, 1999, 5:43:30 PM (26 years ago)

Author:

ansari

Message:

versions templatees, NdataBlocks etc. 15-OCT-99-GLM

Location:

trunk/SophyaLib/Samba

Files:

• circle.cc (modified) (3 diffs)
• circle.h (modified) (3 diffs)
• exclude (modified) (1 diff)
• localmap.cc (modified) (5 diffs)
• localmap.h (modified) (2 diffs)
• mlobe.cc (modified) (1 diff)
• mlobe.h (modified) (1 diff)
• pixelmap.h (modified) (1 diff)
• spheregorski.cc (modified) (24 diffs)
• spheregorski.h (modified) (1 diff)
• spherethetaphi.cc (modified) (26 diffs)
• spherethetaphi.h (modified) (1 diff)
• sphericalmap.h (modified) (2 diffs)

trunk/SophyaLib/Samba/circle.cc

@@ -120 +120 @@
 }
 
-UnitVector Circle::ConvCircleSphere(double psi) const 
+UnitVector Circle::ConvToSphere(double psi) const 
 {
   psi=mod(psi,pi2);
@@ -157 +157 @@
 {
   psi=mod(psi,pi2);
-  return ConvCircleSphere(psi).EPhi();
+  return ConvToSphere(psi).EPhi();
 }
 
@@ -163 +163 @@
 {
   psi=mod(psi,pi2);
-  return ConvCircleSphere(psi).ETheta();
+  return ConvToSphere(psi).ETheta();
 }
 

trunk/SophyaLib/Samba/circle.h

@@ -7 +7 @@
 #include "unitvector.h"
 #include "utilgeom.h"
+#include "geometry.h"
 
-class Circle 
+class Circle : public Geometry
 {
   
@@ -27 +28 @@
   void SetApertureAngle(const Circle&);
   
-  // psi contient les 4 valeurs des angles d'intersection. -1 si les cercles ne se coupent pas
+  // psi contient les 4 valeurs des angles d intersection. -1 si les cercles ne se coupent pas
   // voir la numerotation dans le .cc
   bool Intersection(const Circle&, double* psi) const;
 
   // donne le UnitVector correspondant a une position donnee sur le cercle
-  UnitVector ConvCircleSphere(double psi) const;
+  UnitVector ConvToSphere(double psi) const;
 
   // donne le UnitVector correspondant la tangente au cercle en une position donnee sur le cercle
@@ -40 +41 @@
   UnitVector EPhi(double psi) const;
 
-  // donne l'autre vecteur tangent (orthogonal a EPhi)
+  // donne l autre vecteur tangent (orthogonal a EPhi)
   UnitVector ETheta(double psi) const;
 
-   // donne l'angle de separation dans [0,2Pi] en une position donnee sur le cercle et EPhi
+   // donne l angle de separation dans [0,2Pi] en une position donnee sur le cercle et EPhi
   double SepAngleTanEPhi02PI(double psi) const;
 

trunk/SophyaLib/Samba/exclude

@@ +1 @@
+lambuilder.cc

trunk/SophyaLib/Samba/localmap.cc

@@ -1 +1 @@
 #include "localmap.h"
+#include "nbmath.h"
+#include <complex>
+#include "piocmplx.h"
+
+#include <string.h>
 #include <iostream.h>
-#include "nbmath.h"
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template LocalMap<double>
+#pragma define_template LocalMap<float>
+#pragma define_template LocalMap< complex<double> >
+#pragma define_template LocalMap< complex<float> >
+#pragma define_template FIO_LocalMap<double>
+#pragma define_template FIO_LocalMap<float>
+#pragma define_template FIO_LocalMap< complex<float> >
+#pragma define_template FIO_LocalMap< complex<double> >
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class LocalMap<double>;
+template class LocalMap<float>;
+template class LocalMap< complex<double> >;
+template class LocalMap< complex<float> >;
+template class FIO_LocalMap<double>;
+template class FIO_LocalMap<float>;
+template class FIO_LocalMap< complex<float> >;
+template class FIO_LocalMap< complex<double> >;
+#endif
+
 //*****************************************************************************
 //++ 
@@ -52 +78 @@
 //--
 //++
-
-LocalMap::LocalMap()
-
-//--
-
-{InitNul();}
-
-//++
-
-LocalMap::LocalMap(int_4 nx, int_4 ny) : mSzX_(nx), mSzY_(ny)
-
+template<class T>
+LocalMap<T>::LocalMap()
+//
+// Constructeur par defaut 
 //--
 {
   InitNul();
-  mNPix_=nx*ny;
-  mPix_ = new r_8[mNPix_];
-}
-
-
-//++
-// Titre        Destructeur
-//--
-//++
-LocalMap::~LocalMap()
-
-//--
-{
-  Clear();
-}
-
-
-
-void        LocalMap::WriteSelf(POutPersist& s) const
-{
-  char strg[256];
-  if (mInfo_) {sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d HasInfo",  mSzX_, mSzY_);} else  {
-    sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d ",  mSzX_, mSzY_);
-  }         
-  s.PutLine(strg);
-  if (mInfo_)  s << (*mInfo_);
-  s.PutI4(mSzX_);
-  s.PutI4(mSzY_);
-  s.PutI4(mNPix_);
-  s.PutI4(x0_);
-  s.PutI4(y0_);
-  s.PutI4(originFlag_);
-  s.PutI4(SzFlag_);
-  s.PutR4(cos_angle_);
-  s.PutR4(sin_angle_);
-  s.PutR8(theta0_);
-  s.PutR8( phi0_);
-  s.PutR8( tgAngleX_);
-  s.PutR8(tgAngleY_);
-  s.PutR8s(mPix_, mNPix_);
-}
-void        LocalMap::ReadSelf(PInPersist& s)
-{
-  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(mSzX_);
-  s.GetI4(mSzY_);
-  s.GetI4(mNPix_);
-  s.GetI4(x0_);
-  s.GetI4(y0_);
-  s.GetI4(originFlag_);
-  s.GetI4(SzFlag_);
-  s.GetR4(cos_angle_);
-  s.GetR4(sin_angle_);
-  s.GetR8(theta0_);
-  s.GetR8( phi0_);
-  s.GetR8( tgAngleX_);
-  s.GetR8(tgAngleY_);
-  mPix_ = new r_8[mNPix_+1];
-  s.GetR8s(mPix_, mNPix_);
+}
+
+//++
+template<class T>
+LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
+//
+// Constructeur
+//--
+{
+  InitNul();
+  nPix_= nx*ny;
+  pixels_.ReSize(nPix_);
+  pixels_.Reset();
+}
+
+//++
+template<class T>
+LocalMap<T>::LocalMap(const LocalMap<T>& lm)
+//
+// Constructeur de recopie
+//--
+{
+  cout<<" LocalMap:: Appel du constructeur de recopie " << endl;
+
+  if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);
+  nSzX_= lm.nSzX_;
+  nSzY_= lm.nSzY_;
+  nPix_= lm.nPix_;
+  originFlag_= lm.originFlag_;
+  extensFlag_= lm.extensFlag_;
+  x0_= lm.x0_;
+  y0_= lm.y0_;
+  theta0_= lm.theta0_;
+  phi0_= lm.phi0_;
+  angle_= lm.angle_;
+  cos_angle_= lm.cos_angle_;
+  sin_angle_= lm.sin_angle_;
+  angleX_= lm.angleX_;
+  angleY_= lm.angleY_;
+  tgAngleX_= lm.tgAngleX_;
+  tgAngleY_= lm.tgAngleY_;
+  pixels_= lm.pixels_;
+}
+
+//++
+template<class T>
+LocalMap<T>::~LocalMap()
+//
+// Destructeur
+//--
+{
+  InitNul();
+}
+
+//++
+template<class T>
+void LocalMap<T>::InitNul()
+//
+// Initialise à zéro certaines variables internes
+//--
+{ 
+  originFlag_= false;
+  extensFlag_= false;
+  cos_angle_= 1.0;
+  sin_angle_= 0.0;
+  pixels_.Reset();
+}
+
+//++
+template<class T>
+void LocalMap<T>::ReSize(int_4 nx, int_4 ny) 
+//
+// Redimensionne l'espace de stokage des pixels
+//--
+{
+  InitNul();
+  nSzX_ = nx;
+  nSzY_ = ny;
+  nPix_= nx*ny;
+  pixels_.ReSize(nPix_);
+  pixels_.Reset();
 }
 
@@ -142 +174 @@
 
 //++
-int_4 LocalMap::NbPixels() const
-
+template<class T>
+int_4 LocalMap<T>::NbPixels() const
+//
 //    Retourne le nombre de pixels du découpage
 //--
 {
-  return(mNPix_);
-}
-
-//++
-r_8& LocalMap::PixVal(int_4 k)
-
+  return(nPix_);
+}
+
+//++
+template<class T>
+T& LocalMap<T>::PixVal(int_4 k)
+//
 //    Retourne la valeur du contenu du pixel d'indice k
 //--
 {
-  if ( (k<0) || (k >= mNPix_) ) {
-    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-    //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-  THROW(rangeCheckErr);
-
-    //throw "LocalMap::PIxVal Pixel index out of range";
-  }
-  //return (*(pixel_->Data()+k));
-  return(mPix_[k]);
-}
-
-//++
-
-r_8 const& LocalMap::PixVal(int_4 k) const
-
+  if((k < 0) || (k >= nPix_)) 
+    {
+      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+      //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      THROW(rangeCheckErr);
+      //throw "LocalMap::PIxVal Pixel index out of range";
+    }
+  return(pixels_(k));
+}
+
+//++
+
+template<class T>
+T const& LocalMap<T>::PixVal(int_4 k) const
+//
 //    Retourne la valeur du contenu du pixel d'indice k
 //--
 {
-  if ( (k<0) || (k >= mNPix_) ) {
-    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
-  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-
-    throw "LocalMap::PIxVal Pixel index out of range"; 
-  } 
-  //return *(pixel_->Data()+k);
-  return(mPix_[k]);
-}
-
-//++
-int_4 LocalMap::PixIndexSph(float theta, float phi) const
-
+  if((k < 0) || (k >= nPix_)) 
+    {
+      cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      
+      throw "LocalMap::PIxVal Pixel index out of range"; 
+    } 
+  return *(pixels_.Data()+k);
+}
+
+//++
+template<class T>
+int_4 LocalMap<T>::PixIndexSph(float theta, float phi) const
+//
 //    Retourne l'indice du pixel vers lequel pointe une direction définie par 
 //    ses coordonnées sphériques
@@ -192 +227 @@
 {
   int i,j;
-  if (!(originFlag_==1) || !(SzFlag_==1) ) {
-    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
-    exit(0);
-  }
-  // theta et phi en coordonnees relatives (on se ramene a une situation 
-  // par rapport au plan de reference)
+  if(!(originFlag_) || !(extensFlag_)) 
+    {
+      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
+      exit(0);
+    }
+
+  // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)
   float theta_aux=theta;
   float phi_aux=phi;
   UserToReference(theta_aux, phi_aux);
+
   // coordonnees dans le plan local en unites de pixels
   float x,y;
-  AngleProjToPix(theta_aux, phi_aux, x,y);
-
-  float xmin=-x0_-0.5;
-  float xmax=xmin+mSzX_;
-  if( (x >  xmax) || (x < xmin ) ) return(-1);
-  float xcurrent=xmin;
-  for( i=0; i < mSzX_; i++ ) {
-    xcurrent += 1.;
-    if( x < xcurrent ) break;
-  }
-  float ymin=-y0_-0.5;
-  float ymax=ymin+mSzY_;
-  if( (y >  ymax) || (y < ymin ) ) return(-1);
-  float ycurrent=ymin;
-  for( j=0; j < mSzY_; j++ ) {
-    ycurrent += 1.;
-    if( y < ycurrent ) break;
-  }
-
-
-  return (j*mSzX_+i);
-}
-
-//++
-void  LocalMap::PixThetaPhi(int_4 k, float& theta, float& phi) const
-
+  AngleProjToPix(theta_aux, phi_aux, x, y);
+
+  float xmin= -x0_-0.5;
+  float xmax= xmin+nSzX_;
+  if((x >  xmax) || (x < xmin)) return(-1);
+  float xcurrent= xmin;
+  for(i = 0; i < nSzX_; i++ ) 
+    {
+      xcurrent += 1.;
+      if( x < xcurrent ) break;
+    }
+  float ymin= -y0_-0.5;
+  float ymax= ymin+nSzY_;
+  if((y >  ymax) || (y < ymin)) return(-1);
+  float ycurrent= ymin;
+  for(j = 0; j < nSzY_; j++ ) 
+    {
+      ycurrent += 1.;
+      if( y < ycurrent ) break;
+    }
+  return (j*nSzX_+i);
+}
+
+//++
+template<class T>
+void LocalMap<T>::PixThetaPhi(int_4 k, float& theta, float& phi) const
+//
 //    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
 //--
 {
-  if (!(originFlag_==1) || !(SzFlag_==1) ) {
-    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
-    exit(0);
-  }
+  if(!(originFlag_) || !(extensFlag_)) 
+    {
+      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
+      exit(0);
+    }
+
   int i,j;
   Getij(k,i,j);
-  float X=float(i-x0_);
-  float Y=float(j-y0_);
+
+  float X= float(i-x0_);
+  float Y= float(j-y0_);
   // situation de ce pixel dans le plan de reference
-  float x=X*cos_angle_-Y*sin_angle_;
-  float y=X*sin_angle_+Y* cos_angle_;
+  float x= X*cos_angle_-Y*sin_angle_;
+  float y= X*sin_angle_+Y* cos_angle_;
   // projection sur la sphere
-  PixProjToAngle(x, y,theta, phi);
+  PixProjToAngle(x, y, theta, phi);
   // retour au plan utilisateur
   ReferenceToUser(theta, phi);
@@ -250 +291 @@
  
 //++
-r_8    LocalMap::PixSolAngle(int_4 k) const
-//    Pixel Solid angle  (steradians)
-//    All the pixels have not necessarly the same size in (theta, phi)
-//    because of the projection scheme which is not yet fixed. 
-//    
+template<class T>
+r_8 LocalMap<T>::PixSolAngle(int_4 k) const
+//
+// Pixel Solid angle  (steradians)
+// All the pixels have not necessarly the same size in (theta, phi)
+// because of the projection scheme which is not yet fixed.   
 //--
 {
   int i,j;
   Getij(k,i,j);
-  float X=float(i-x0_);
-  float Y=float(j-y0_);
-  float XR=X+float(i)*0.5;
-  float XL=X-float(i)*0.5;
-  float YU=Y+float(j)*0.5;
-  float YL=Y-float(j)*0.5;
+  float X= float(i-x0_);
+  float Y= float(j-y0_);
+  float XR= X+float(i)*0.5;
+  float XL= X-float(i)*0.5;
+  float YU= Y+float(j)*0.5;
+  float YL= Y-float(j)*0.5;
   // situation  dans le plan de reference
-  float x0=XL*cos_angle_-YL*sin_angle_;
-  float y0=XL*sin_angle_+YL* cos_angle_;
-  float xa=XR*cos_angle_-YL*sin_angle_;
-  float ya=XR*sin_angle_+YL* cos_angle_;
-  float xb=XL*cos_angle_-YU*sin_angle_;
-  float yb=XL*sin_angle_+YU* cos_angle_;
+  float x0= XL*cos_angle_-YL*sin_angle_;
+  float y0= XL*sin_angle_+YL*cos_angle_;
+  float xa= XR*cos_angle_-YL*sin_angle_;
+  float ya= XR*sin_angle_+YL*cos_angle_;
+  float xb= XL*cos_angle_-YU*sin_angle_;
+  float yb= XL*sin_angle_+YU*cos_angle_;
   // projection sur la sphere
   float tet0,phi0,teta,phia,tetb,phib;
-  PixProjToAngle(x0, y0,tet0, phi0);
-  PixProjToAngle(xa, ya,teta, phia);
-  PixProjToAngle(xb, yb,tetb, phib);
+  PixProjToAngle(x0, y0, tet0, phi0);
+  PixProjToAngle(xa, ya, teta, phia);
+  PixProjToAngle(xb, yb, tetb, phib);
   // angle solide
-  float sol=fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
+  float sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
   return r_8(sol);
 }
 
-
-//++
-void  LocalMap::SetOrigin(float theta0, float phi0, float angle)
-
-//    
-//--
-{
-  theta0_=theta0;
-  phi0_=phi0;
-  x0_=mSzX_/2;
-  y0_=mSzY_/2;
-  angle=angle*Pi/180.;
-  cos_angle_=cos(angle);
-  sin_angle_=sin(angle);
-  originFlag_ =1; 
-}
-
+//++
+template<class T>
+void  LocalMap<T>::SetOrigin(float theta0, float phi0, float angle)
+//
+// Fixe la repere de reference ( angles en degres) 
+//--
+{
+  theta0_= (double)theta0;
+  phi0_  = (double)phi0;
+  angle_ = (double)angle;
+  x0_= nSzX_/2;
+  y0_= nSzY_/2;
+  cos_angle_= cosdf(angle);
+  sin_angle_= sindf(angle);
+  originFlag_= true; 
+  cout << " LocalMap:: set origin 1 done" << endl;
+}
  
 //++
-void  LocalMap::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
-
-//    
-//--
-{
-  theta0_=theta0;
-  phi0_=phi0;
-  x0_=x0;
-  y0_=y0;
-  angle=angle*Pi/180.;
-  cos_angle_=cos(angle);
-  sin_angle_=sin(angle);
-  originFlag_ =1; 
+template<class T>
+void  LocalMap<T>::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
+//
+// Fixe le repere de reference (angles en degres)  
+//--
+{
+  theta0_= (double)theta0;
+  phi0_  = (double)phi0;
+  angle_ = (double)angle;
+  x0_= x0;
+  y0_= y0;
+  cos_angle_= cosdf(angle);
+  sin_angle_= sindf(angle);
+  originFlag_= true;
+  cout << " LocalMap:: set origin 2 done" << endl;
 } 
 
 //++
-void  LocalMap::SetSize(float angleX, float angleY)
-
-//    
-//--
-{
- tgAngleX_=r_8(tan(angleX*Pi/360.));
- tgAngleY_=r_8(tan(angleY*Pi/360.));
- SzFlag_=1;
-}
-//++
-void  LocalMap::Project(SphericalMap& sphere) const
-
-//    Projection to spherical map   
-//--
-{
-  for (int m=0; m<mNPix_; m++) {
-    float theta,phi;
-    PixThetaPhi(m,theta,phi);
-    sphere(theta,phi)=mPix_[m];
-    //    cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
-  }
-}
-
-void LocalMap::InitNul()
-//
-//    initialise à zéro certaines variables internes
-{ 
-originFlag_=0;
-SzFlag_=0;
-cos_angle_=1.;
-sin_angle_=0.;
-mPix_=NULL;
-}
-
-void LocalMap::Clear()
-//
-//    
-{ 
-  if (mPix_)   delete[] mPix_;
-}
-
-void    LocalMap::Getij(int k, int& i, int& j) const
-{
-  i=(k+1)%mSzX_-1;
-  if (i==-1) i=mSzX_-1;
-  j=(k-i+2)/mSzX_;
-}
-
-void  LocalMap::ReferenceToUser(float &theta, float &phi) const
-
-//       
-{
-  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
-    //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
-  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-
-    throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
-  } 
+template<class T>
+void LocalMap<T>::SetSize(float angleX, float angleY)
+//
+// Fixe l'extension de la carte (angles en degres)
+//--
+{
+  angleX_= (double)angleX;
+  angleY_= (double)angleY;
+
+  //tgAngleX_= T(tan(angleX*Pi/360.));
+  //tgAngleY_= T(tan(angleY*Pi/360.));
+
+
+  // tangente de la moitie de l'ouverture angulaire totale
+  tgAngleX_= tand(0.5*angleX_);
+  tgAngleY_= tand(0.5*angleY_);
+
+  extensFlag_= true;
+  cout << " LocalMap:: set extension done" << endl;
+}
+
+//++
+template<class T>
+void LocalMap<T>::Project(SphericalMap<T>& sphere) const
+//
+// Projection to spherical map   
+//--
+{
+  for(int m = 0; m < nPix_; m++) 
+    {
+      float theta,phi;
+      PixThetaPhi(m,theta,phi);
+      sphere(theta,phi)= pixels_(m);
+      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
+    }
+}
+
+//++
+template<class T>
+void LocalMap<T>::Getij(int k, int& i, int& j) const
+//
+//--
+{
+  i= (k+1)%nSzX_-1;
+  if(i == -1) i= nSzX_-1;
+  j= (k-i+2)/nSzX_;
+}
+
+//++
+template<class T>
+void  LocalMap<T>::ReferenceToUser(float &theta, float &phi) const
+//
+// --     
+{
+  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi) 
+    {
+      //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
+      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
+    } 
+
   //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
-  theta=(r_4)theta0_+theta-(r_4)Pi*0.5;
-  if (theta<0.) {
-    theta=-theta;
-    phi+=(r_4)Pi;
-  } else 
-    if (theta>(r_4)Pi) {
-      theta=2.*(r_4)Pi-theta;
-      phi+=(r_4)Pi;
+  theta= (r_4)theta0_*Pi/180.+theta-(r_4)Pi*0.5;
+  if(theta < 0.) 
+    {
+      theta= -theta;
+      phi += (r_4)Pi;
     } 
-  
-  phi=(r_4)phi0_+phi;
-  while (phi>=2.*(r_4)Pi) phi-=2.*(r_4)Pi;
-  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
-    cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
-    cout << " theta= " << theta << " phi= " << phi << endl;
-    exit(0);
-  }
+  else 
+    {
+      if(theta > (r_4)Pi) 
+        {
+          theta= 2.*(r_4)Pi-theta;
+          phi += (r_4)Pi;
+        } 
+    }
+
+  phi= (r_4)phi0_*Pi/180.+phi;
+  while(phi >= 2.*(r_4)Pi) phi-=2.*(r_4)Pi;
+
+  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi) 
+    {
+      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
+      cout << " theta= " << theta << " phi= " << phi << endl;
+      exit(0);
+    }
   //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
 }
 
-void  LocalMap::UserToReference(float &theta, float &phi) const
-
-//       
-{
-  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
-    cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
-  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
-
-    throw "LocalMap::UserToReference (theta,phi) out of range"; 
-  } 
-  float phi1=phi-(r_4)phi0_;
-  if (phi1<0.) phi1+=2.*(r_4)Pi;
-
-  float theta1=theta-(r_4)theta0_+(r_4)Pi*0.5;
-  if (theta1<0.) {
-    theta=-theta1;
-    phi1+=(r_4)Pi;
-  } else 
-    if (theta1>(r_4)Pi) {
-      theta=2.*(r_4)Pi-theta1;
-      phi1+=(r_4)Pi;
-    }
-  while (phi1>=2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
-  phi=phi1;
-  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
-    cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
-    cout << " theta= " << theta << " phi= " << phi << endl;
-    exit(0);
-  }
-}
-
-void LocalMap::PixProjToAngle(float x, float y,float &theta, float &phi) const {
-
-// (x,y) representent les coordonnees en unites de pixels d'un point DANS 
-// LE PLAN DE REFERENCE. On recupere (theta,phi) par rapport au repere "absolu"
-// theta=pi/2,phi=0.
-
-theta=(r_4)Pi*0.5-atan(2*y* tgAngleY_/(float) mSzY_);
-phi=atan2(2*x* tgAngleX_,(float) mSzX_);
-if( phi< 0. ) phi+= DeuxPi;
-}
-
-void  LocalMap::AngleProjToPix(float theta, float phi, float& x, float& y) const {
-
+//++
+template<class T>
+void  LocalMap<T>::UserToReference(float &theta, float &phi) const
+//
+// --     
+{
+  if(theta > (r_4)Pi || theta < 0.  || phi<0. || phi >= 2*(r_4)Pi ) 
+    {
+      cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
+      //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+      throw "LocalMap::UserToReference (theta,phi) out of range"; 
+    }
+
+  float phi1=phi-(r_4)phi0_*Pi/180.;
+  if(phi1 < 0.) phi1+=2.*(r_4)Pi;
+
+  float theta1= theta-(r_4)theta0_*Pi/180.+(r_4)Pi*0.5;
+  if(theta1 < 0.) 
+    {
+      theta= -theta1;
+      phi1+= (r_4)Pi;
+    } 
+  else
+    { 
+      if(theta1 > (r_4)Pi) 
+        {
+          theta= 2.*(r_4)Pi-theta1;
+          phi1+= (r_4)Pi;
+        }
+    }
+
+  while(phi1 >= 2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
+  phi= phi1;
+  if(theta > (r_4)Pi || theta < 0.  || phi < 0. || phi >= 2*(r_4)Pi ) 
+    {
+      cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
+      cout << " theta= " << theta << " phi= " << phi << endl;
+      exit(0);
+    }
+}
+
+//++
+template<class T>
+void LocalMap<T>::PixProjToAngle(float x, float y, float& theta, float& phi) const
+//
+// (x,y) representent les coordonnees en unites de pixels d'un point DANS  LE PLAN DE REFERENCE. 
+// On recupere (theta,phi) par rapport au repere "absolu" theta=pi/2 et phi=0.
+//--
+{
+  theta= (r_4)Pi*0.5-atan(2*y*tgAngleY_/(float)nSzY_);
+  phi= atan2(2*x*tgAngleX_,(float)nSzX_);
+  if(phi < 0.) phi += DeuxPi;
+}
+
+//++
+template<class T>
+void LocalMap<T>::AngleProjToPix(float theta, float phi, float& x, float& y) const
+//
 // (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere 
 // (i,j)  DANS LE PLAN DE REFERENCE. 
-  if (phi>=(r_4)Pi) phi-= DeuxPi;
+//--
+{
+  if(phi >= (r_4)Pi) phi-= DeuxPi;
   //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
-  y=0.5*mSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
-  x=0.5*mSzX_*tan(phi)/tgAngleX_;
-}
-
-
-
-
-
+  y= 0.5*nSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
+  x= 0.5*nSzX_*tan(phi)/tgAngleX_;
+}
+
+template<class T>
+void LocalMap<T>::print(ostream& os) const
+{
+  os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;
+  if(LocalMap_isDone())
+    {
+      os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;
+      os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;
+      os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;
+      os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;
+      os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;
+    }
+
+  os << " contenu de pixels : ";
+  for(int i=0; i < nPix_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os <<  pixels_(i) <<", ";
+    }
+  os << endl;
+}
+
+//*******************************************************************
+// class FIO_LocalMap<T>
+//  Les objets delegues pour la gestion de persistance 
+//*******************************************************************
+
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap()
+{
+  dobj= new LocalMap<T>;
+  ownobj= true;
+}
+
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(string const& filename)
+{
+  dobj= new LocalMap<T>;
+  ownobj= true;
+  Read(filename);
+}
+
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
+{
+  dobj= new LocalMap<T>(obj);
+  ownobj= true;
+}
+
+template <class T>
+FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)
+{
+  dobj= obj;
+  ownobj= false;
+}
+
+template <class T>
+FIO_LocalMap<T>::~FIO_LocalMap()
+{
+  if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_LocalMap<T>::DataObj()
+{
+  return(dobj);
+}
+
+template <class T>
+void FIO_LocalMap<T>::ReadSelf(PInPersist& is)
+{
+  cout << " FIO_LocalMap:: ReadSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      dobj= new LocalMap<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 nSzX;
+  is.GetI4(nSzX);
+  dobj->setSize_x(nSzX);
+
+  int_4 nSzY;
+  is.GetI4(nSzY);
+  dobj->setSize_y(nSzY);
+
+  int_4 nPix;
+  is.GetI4(nPix);
+  dobj->setNbPixels(nPix);
+
+  string ss("local mapping is done");
+  string sso;
+  is.GetStr(sso);
+  if(sso == ss)
+    {
+      cout<<" ReadSelf:: local mapping"<<endl;
+      int_4 x0, y0;
+      float theta, phi, angle;
+      is.GetI4(x0);
+      is.GetI4(y0);
+      is.GetR4(theta);
+      is.GetR4(phi);
+      is.GetR4(angle);
+      dobj->SetOrigin(theta, phi, x0, y0, angle);
+
+      float angleX, angleY;
+      is.GetR4(angleX);
+      is.GetR4(angleY);
+      dobj->SetSize(angleX, angleY);
+    }
+
+  T* pixels= new T[nPix];
+  PIOSReadArray(is, pixels, nPix);
+  dobj->setDataBlock(pixels, nPix);
+  delete [] pixels;
+}
+
+template <class T>
+void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
+{
+  cout << " FIO_LocalMap:: WriteSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      cout << " WriteSelf:: dobj= null " << endl;
+      return;
+    }
+
+  char strg[256];
+  int_4 nSzX= dobj->Size_x();
+  int_4 nSzY= dobj->Size_y();
+  int_4 nPix= dobj->NbPixels();
+ 
+  if(dobj->ptrInfo()) 
+    {
+      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d HasInfo",nSzX,nSzY);
+      os.PutLine(strg);
+      os << dobj->Info();
+    }
+  else 
+    { 
+      sprintf(strg,"LocalMap: NPixX=%6d  NPixY=%9d ",nSzX,nSzY);
+      os.PutLine(strg);  
+    }
+
+  os.PutI4(nSzX);
+  os.PutI4(nSzY);
+  os.PutI4(nPix);
+
+  if(dobj->LocalMap_isDone())
+    {
+      string ss("local mapping is done");
+      os.PutStr(ss);
+      int_4 x0, y0;
+      float theta, phi, angle;
+      dobj->Origin(theta, phi, x0, y0, angle);
+      os.PutI4(x0);
+      os.PutI4(y0);
+      os.PutR4(theta);
+      os.PutR4(phi);
+      os.PutR4(angle);
+
+      float angleX, angleY;
+      dobj->Aperture(angleX, angleY);
+      os.PutR4(angleX);
+      os.PutR4(angleY);
+    }
+  else
+    {
+      string ss("no local mapping");
+      os.PutStr(ss);
+    }
+
+  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
+}
+

trunk/SophyaLib/Samba/localmap.h

@@ -4 +4 @@
 #include "pixelmap.h"
 #include "sphericalmap.h"
+#include "ndatablock.h"
+
+#include "anydataobj.h"
+#include "ppersist.h"
 
 // A local map of a region of the sky, in cartesian coordinates.
@@ -13 +17 @@
 //    indice de colonne et j indice de ligne. La carte est supposee resider 
 //    dans un plan tangent, dont le point de tangence est repere (x0,y0) dans 
-//    la carte et (theta0, phi0) sur la sphere celeste. L'extension de la 
+//    la carte et (theta0, phi0) sur la sphere celeste. L extension de la 
 //    carte est definie par les valeurs de deux angles couverts respectivement
 //    par la totalite des pixels en x de la carte et la totalite des pixels 
 //    en y. (SetSize()).
 //    On considere un "plan de reference" : plan tangent a la sphere celeste 
-//    aux angles theta=Pi/2 et phi=0. Dans ce plan L'origine des coordonnees 
-//    est le point de tangence. L'axe Ox est la tangente parallele a 
-//    l'equateur, dirige vers les phi croissants, l'axe Oy est parallele 
+//    aux angles theta=Pi/2 et phi=0. Dans ce plan L origine des coordonnees 
+//    est le point de tangence. L axe Ox est la tangente parallele a 
+//    lequateur, dirige vers les phi croissants, l axe Oy est parallele 
 //    au meridien, dirige vers le pole nord.
 //    De maniere interne a la classe une carte est definie dans ce plan de 
-//    reference et transportee  jusqu'au point (theta0, phi0) de sorte que les //    axes restent paralleles aux meridiens et paralleles. L'utilisateur peut 
+//    reference et transportee  jusqu au point (theta0, phi0) de sorte que les //    axes restent paralleles aux meridiens et paralleles. L utilisateur peut 
 //    definir sa carte selon un repere en rotation par rapport au repere de 
-//    reference (par l'angle entre le parallele et l'axe Ox souhaite -- 
+//    reference (par l angle entre le parallele et l axe Ox souhaite -- 
 //    methode SetOrigin(...))
 
 
+// ***************** Class LocalMap *****************************
 
+template<class T>
+class LocalMap : public PixelMap<T>, public AnyDataObj 
+{
 
+public:
 
+LocalMap();
+LocalMap(int_4 nx, int_4 ny);
+LocalMap(const LocalMap<T>& lm);
+virtual ~LocalMap();
 
+// ---------- Overloading of () to access pixel number k ----
 
-class LocalMap : public PixelMap {
-public:
-   LocalMap();
-   LocalMap(int_4 nx, int_4 ny);
-   virtual ~LocalMap();
-   // Overloading of () to access pixel number k.
-inline  r_8&        operator()(int_4 k)
-                          {return(PixVal(k));}
-inline  r_8 const&  operator()(int_4 k) const
-                          {return(PixVal(k));}
-inline  r_8&        operator()(int ix, int iy)
-                { return PixVal(iy*mSzX_+ix) ; };
-inline  r_8 const&        operator()(int ix, int iy) const
-                { return PixVal(iy*mSzX_+ix) ; };
+inline T& operator()(int_4 k) {return(PixVal(k));}
+inline T const& operator()(int_4 k) const {return(PixVal(k));}
+inline T& operator()(int ix, int iy) {return PixVal(iy*nSzX_+ix);};
+inline T const& operator()(int ix, int iy) const {return PixVal(iy*nSzX_+ix);};
    
-// ---------- Persistence handling
+// ---------- Definition of PixelMap abstract methods -------
+
+/* return/set the number of pixels */
+virtual int_4 NbPixels() const;
+inline void setNbPixels(int_4 n) {nPix_= n;}
+  
+/* return the value of pixel number k */
+virtual T& PixVal(int_4 k);
+virtual T const& PixVal(int_4 k) const;
+
+/* return the index of pixel at (theta,phi) */
+virtual int_4 PixIndexSph(float theta, float phi) const;
    
-   enum {classId = 0xF003 };
-int_4               ClassId() const        { return classId; }
+/* return the spherical coordinates of center of pixel number k */
+virtual void PixThetaPhi(int_4 k, float& theta, float& phi) const;
 
-virtual void        WriteSelf(POutPersist&) const;
-virtual void        ReadSelf(PInPersist&);
+/* return the Pixel Solid angle  (steradians) */
+virtual r_8 PixSolAngle(int_4 k) const; 
 
-// ---------- Definition of PixelMap abstract methods
+// ---------- Specific methods ------------------------------
 
-   // Number of pixels
-   virtual int_4       NbPixels() const;
-   
-   // Value of pixel number k
-   virtual r_8&        PixVal(int_4 k);
-   virtual r_8 const&  PixVal(int_4 k) const;
+void ReSize(int_4 nx, int_4 ny);
 
-   // Index of pixel at (theta,phi)
-   virtual int_4       PixIndexSph(float theta, float phi) const;
-   
-   // Spherical coordinates of center of pixel number k
-   virtual void        PixThetaPhi(int_4 k, float& theta, float& phi) const;
-   // Pixel Solid angle  (steradians)
-   virtual r_8         PixSolAngle(int_4 k) const; 
+inline virtual char* TypeOfMap() const {return "LOCAL";};
+ 
+/* Origin (with angle between x axis and phi axis, in degrees)  x0,y0  the default: middle of map*/
+virtual void SetOrigin(float theta=90., float phi=0., float angle=0.);
+virtual void SetOrigin(float theta,float phi,int_4 x0,int_4 y0,float angle=0.);
 
-// ---------- Specific methods
+/* Pixel size (degres) */ 
+virtual void SetSize(float angleX, float angleY);
 
-   // Origin (with angle between x axis and phi axis, in degrees)
-    virtual void  SetOrigin(float theta0, float phi0, float angle=0.); // x0,y0 default: middle of map
-  virtual void        SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle=0.);
-   // Pixel size (degres) 
-   virtual void        SetSize(float angleX, float angleY);
+/* Check to see if the local mapping is done */
+inline bool LocalMap_isDone() const {return(originFlag_ && extensFlag_);};
 
-   // Projection to/from spherical map
-   //virtual void        Extract(SphericalMap const& sphere);
-   virtual void        Project(SphericalMap& sphere) const;
+/* Projection to/from spherical map */
+virtual void Project(SphericalMap<T>& sphere) const;
   
-   // There should be a more complex algorithm somewhere to combine *several*
-   // local maps to a full sphere.
-   //   -> static method, or separate class
+/* There should be a more complex algorithm somewhere to combine *several* local maps to a full sphere.
+      -> static method, or separate class */
+  
+/* provides a integer characterizing the pixelization refinement  (here : number of pixels) */
+inline virtual int_4 SizeIndex() const {return(nPix_);}
+inline int_4 Size_x() const {return nSzX_;}
+inline void setSize_x(int_4 n) {nSzX_= n;}
+inline int_4 Size_y() const {return nSzY_;}
+inline void setSize_y(int_4 n) {nSzY_= n;}
 
-   void Pixelize(int_4,int_4); // Allocate pixel array
+inline void Origin(float& theta, float& phi,int& x0, int& y0, float& angle) const {theta= (float)theta0_; phi= (float)phi0_; x0= x0_; y0= y0_;angle= (float)angle_;}
 
-// ------------- méthodes internes ----------------------
+inline void Aperture(float& anglex, float& angley) const {anglex= (float)angleX_; angley= (float)angleY_;}
+
+/* retourne le pointeur vers/remplit  le vecteur des contenus des pixels */ 
+inline const NDataBlock<T>* getDataBlock() const {return (&pixels_);}
+inline void setDataBlock(T* data, int_4 n) {pixels_.FillFrom(n,data);}
+
+/* impression */ 
+void print(ostream& os) const;
+
+// ---------- Méthodes internes -----------------------------
           
 private :
 
-   void      InitNul();
-   void      Clear();
-   void      Getij(int k, int& i, int& j) const;
-   void      ReferenceToUser(float &theta, float &phi) const; 
-   void      UserToReference(float &theta, float &phi) const;
-   void      PixProjToAngle(float x, float y,float &theta, float &phi) const; 
-   void      AngleProjToPix(float theta, float phi, float& x, float& y) const; 
-// ------------- variables internes -----------------------
+void InitNul();
+void Getij(int k, int& i, int& j) const;
+void ReferenceToUser(float &theta, float &phi) const; 
+void UserToReference(float &theta, float &phi) const;
+void PixProjToAngle(float x, float y,float &theta, float &phi) const; 
+void AngleProjToPix(float theta, float phi, float& x, float& y) const; 
 
-   int_4   mSzX_, mSzY_;
-   int_4   mNPix_;
-   int_4   x0_;
-   int_4   y0_;
-   int_4    originFlag_;
-   int_4    SzFlag_;
-   r_4   cos_angle_;
-   r_4   sin_angle_;
-   r_8     theta0_;
-   r_8     phi0_;
-   r_8     omeg_;
-   r_8     tgAngleX_;
-   r_8     tgAngleY_;
-   r_8*    mPix_;
+// ---------- Variables internes ----------------------------
+
+int_4 nSzX_;
+int_4 nSzY_;
+int_4 nPix_;
+bool originFlag_;
+bool extensFlag_;
+int_4 x0_;
+int_4 y0_;
+r_8 theta0_;
+r_8 phi0_;
+r_8 angle_;
+r_4 cos_angle_;
+r_4 sin_angle_;
+r_8 angleX_;
+r_8 angleY_;
+r_8 tgAngleX_;
+r_8 tgAngleY_;
+NDataBlock<T> pixels_;
+};
+
+// ------------- Classe pour la gestion de persistance --
+template <class T>
+class FIO_LocalMap : public PPersist  
+{
+
+public:
+
+FIO_LocalMap();
+FIO_LocalMap(string const & filename); 
+FIO_LocalMap(const LocalMap<T>& obj);
+FIO_LocalMap(LocalMap<T>* obj);
+virtual ~FIO_LocalMap();
+virtual AnyDataObj* DataObj();
+inline operator LocalMap<T>() { return(*dobj); }
+inline LocalMap<T> getObj() { return(*dobj); }
+
+protected :
+
+virtual void ReadSelf(PInPersist&);           
+virtual void WriteSelf(POutPersist&) const;  
+LocalMap<T>* dobj;
+bool ownobj;
 };
 

trunk/SophyaLib/Samba/mlobe.cc

@@ -133 +133 @@
 
 /* --Methode-- */
-double MainLobe::Convol(SphericalMap& sph)
+double MainLobe::Convol(SphericalMap<double>& sph)
 {
 double ret=0.;

trunk/SophyaLib/Samba/mlobe.h

@@ -29 +29 @@
                void     SetDirection(float teta, float phi);
                double   Value(int kpx, float& teta, float& phi);
-               double   Convol(SphericalMap& sph);
+               double   Convol(SphericalMap<double>& sph);
         friend ostream& operator << (ostream& s, MainLobe const& lob);
         inline void     Print() { cout << (*this); } 

trunk/SophyaLib/Samba/pixelmap.h

@@ -15 +15 @@
 //      LocalMap
 
+template<class T>
+class PixelMap  
+{
 
-class PixelMap : public PPersist {
 public:
-   PixelMap():mInfo_(NULL) {}
-   virtual ~PixelMap() {}
+  
+PixelMap():mInfo_(NULL) {};
+virtual ~PixelMap() {};
+  
+// Number of pixels
+virtual int_4 NbPixels() const=0;
    
-   // Number of pixels
-   virtual int_4       NbPixels() const=0;
+// Value of pixel number k
+virtual T& PixVal(int_4 k)=0;
+virtual T const& PixVal(int_4 k) const=0;
    
-   // Value of pixel number k
-   virtual r_8&        PixVal(int_4 k)=0;
-   virtual r_8 const&  PixVal(int_4 k) const=0;
-   
-   // Index of pixel at (theta,phi)
-   virtual int_4       PixIndexSph(float theta, float phi) const=0;
+// Index of pixel at (theta,phi)
+virtual int_4 PixIndexSph(float theta, float phi) const=0;
 
-   // Value of pixel number at (theta,phi)
-   virtual r_8&        PixValSph(float theta, float phi)
+// Value of pixel number at (theta,phi)
+virtual T& PixValSph(float theta, float phi)
                            {return PixVal(PixIndexSph(theta,phi));}
-   virtual r_8 const&  PixValSph(float theta, float phi) const
-                           {return PixVal(PixIndexSph(theta,phi));}                           
-   
-   // Spherical coordinates of center of pixel number k
-   virtual void        PixThetaPhi(int_4 k, float& theta, float& phi) const=0;
+virtual T const& PixValSph(float theta, float phi) const
+                           {return PixVal(PixIndexSph(theta,phi));}
 
-   // Pixel  (steradians)
-   virtual r_8         PixSolAngle(int_4 k) const =0; 
+// Spherical coordinates of center of pixel number k
+virtual void PixThetaPhi(int_4 k, float& theta, float& phi) const=0;
 
-   // Overloading of () to access pixel number k.
-   inline  r_8&        operator()(int_4 k)
-                          {return(PixVal(k));}
-   inline  r_8 const&  operator()(int_4 k) const
-                          {return(PixVal(k));}
+// provides a integer characterizing the pixelization refinement 
+// (depending of the type of the map)              
+virtual int_4 SizeIndex() const=0;
+virtual char* TypeOfMap() const =0;
 
-   // Note : no overloading of (float,float) to access pixel (theta,phi).
-   // overloading of (float,float) in SphericalMap
-   // overloading of (int,int)     in CartesianMap
+// Pixel  (steradians)
+virtual r_8 PixSolAngle(int_4 k) const =0; 
+
+// Overloading of () to access pixel number k.
+inline T& operator()(int_4 k) {return(PixVal(k));}
+inline T const& operator()(int_4 k) const {return(PixVal(k));}
+
+// Note : no overloading of (float,float) to access pixel (theta,phi).
+// overloading of (float,float) in SphericalMap
+// overloading of (int,int)     in CartesianMap
 
 //++
-DVList&  Info()
+DVList& Info()
 //
-//      Renvoie une rM-ifM-irence sur l'objet DVList AssociM-i
+// Renvoie une reference sur l''objet DVList Associe
 //--
-{
-if (mInfo_ == NULL)  mInfo_ = new DVList;
-return(*mInfo_);
-}
+  {
+    if (mInfo_ == NULL)  mInfo_ = new DVList;
+    return(*mInfo_);
+  }
 
-   protected :
+const DVList*  ptrInfo() const
+  {
+    return mInfo_;
+  }
 
-  DVList* mInfo_;        // Infos (variables) attachees 
+protected :
+
+    DVList* mInfo_;        // Infos (variables) attachees 
 };
-
 #endif

trunk/SophyaLib/Samba/spheregorski.cc

@@ -2 +2 @@
 #include "spheregorski.h"
 #include "strutil.h"
-extern "C" {
+#include <complex>
+#include "piocmplx.h"
+
+extern "C" 
+{
 #include <stdio.h>
 #include <stdlib.h>
@@ -8 +12 @@
 }
       
-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 +185 @@
 //++
 
-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 +203 @@
 //--
 {
-  //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 +293 @@
 //--
 {
-  
   // 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 +363 @@
 //--
 {
-
-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 +480 @@
 //--
 {
-  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 +491 @@
 //--
 {
-  return ang2pix_nest(nSide_, double(theta), double(phi));
+  return ang2pix_nest(nSide_,double(theta),double(phi));
 }
 
@@ -497 +497 @@
 /* --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 +505 @@
   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 +522 @@
 }
 
-
-//++
-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 +533 @@
   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 +547 @@
   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 +558 @@
   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 +581 @@
 //
 { 
-  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 +592 @@
   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 +615 @@
   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 +663 @@
   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 +674 @@
   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 +683 @@
   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 +871 @@
   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 +950 @@
   }
     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 +967 @@
     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 +1114 @@
 }
 
-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
+

trunk/SophyaLib/Samba/spheregorski.h

@@ -3 +3 @@
 
 #include "sphericalmap.h"
-#include "cvector.h"
+#include "tvector.h"
+#include "ndatablock.h"
+
+#include "anydataobj.h"
+#include "ppersist.h"
 
 // les attributs de classe pix2x_ et pix2y_ sont relatifs a la traduction des 
-// indices  RING en indices NESTED (ou l'inverse). Ce sont des tableaux 
+// indices  RING en indices NESTED (ou l inverse). Ce sont des tableaux 
 // d'entiers initialises et remplis par le constructeur. Dans la version 
 // FORTRAN de healpix, il s'agissait de tableaux mis en COMMON. Ils etaient 
-// initialises au premier appel d'une conversion d'indice. Je ne peux pas 
+// initialises au premier appel d'une conversion d indice. Je ne peux pas 
 // garder cette procedure car, par exemple, la fonction PixValNest() const 
 // n'autorisera pas la constitution de ces tableaux (a cause du const). Ainsi, 
-// des la creation d'un objet SphereGorski ces tableaux sont calcules, ce qui 
+// des la creation d un objet SphereGorski ces tableaux sont calcules, ce qui 
 // est, certes, inutile si on ne se sert pas des indices NESTED. Ca ne me 
-// rejouit pas, mais c'est le seul moyen que j'ai trouve pour temir compte de 
+// rejouit pas, mais c est le seul moyen que j ai trouve pour temir compte de 
 // toutes les demandes et des options prises.
 // 
 //                                                G. Le Meur 
 
-class SphereGorski : public SphericalMap {
+// ***************** CLASSE SphereGorski *****************************
+
+template<class T>
+class SphereGorski : public SphericalMap<T>,  public AnyDataObj 
+{
+
 public :
-
 
 SphereGorski();
 SphereGorski(int_4 m);
-SphereGorski(char* flnm);
+SphereGorski(const SphereGorski<T>& s);
 virtual ~SphereGorski();
-
-// ------------ Persistence handling
-     enum {classId = 0xF002 };
-
-int_4                  ClassId() const        { return classId; }
-
-virtual void           WriteSelf(POutPersist&) const;
-virtual void           ReadSelf(PInPersist&);
-
-// -----------   FITS IO (1D FITS ARRAY)
-void ReadFits(char flnm[] );
-void WriteFits(char flnm[]);
 
 // ------------------ Definition of PixelMap abstract methods
 
-/* Nombre de pixels du decoupage                                         */
-virtual  int_4           NbPixels() const;
+/* Nombre de pixels du decoupage */
+virtual int_4 NbPixels() const;
+inline void setNbPixels(int_4 n) {nPix_= n;}
 
-/* Valeur du contenu du pixel d'indice "RING" k                                 */
-virtual  r_8&       PixVal(int_4 k);
-virtual  r_8  const&     PixVal(int_4 k) const;
+/* Valeur du contenu du pixel d'indice "RING" k  */
+virtual T& PixVal(int_4 k);
+virtual T const& PixVal(int_4 k) const;
 
+/* Nombre de tranches en theta */ 
 int_4 NbThetaSlices() const;
-void  GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const;
+void GetThetaSlice(int_4 index,r_4& theta,TVector<float>& phi,TVector<T>& value) const;
 
 /* Indice "RING" du pixel vers lequel pointe une direction definie par 
 ses  coordonnees spheriques */                                    
-virtual  int_4           PixIndexSph(float theta, float phi) const;
+virtual int_4 PixIndexSph(float theta, float phi) const;
 
-/* Coordonnees spheriques du milieu du pixel d'indice "RING" k       */
-virtual  void          PixThetaPhi(int_4 k, float& teta, float& phi) const;
+/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
+virtual void PixThetaPhi(int_4 k, float& teta, float& phi) const;
 
-// Pixel Solid angle  (steradians)
-virtual r_8         PixSolAngle(int_4 dummy) const; 
-
+/* Pixel Solid angle  (steradians) */
+virtual r_8 PixSolAngle(int_4 dummy) const;
+inline void setPixSolAngle(r_8 x) {omeg_= x;}
 
 // --------------- Specific methods 
 
-    // NEST indexing
+virtual void Resize(int_4 m);
+
+inline virtual char* TypeOfMap() const {return "RING";};
+
 /* Valeur du contenu du pixel d'indice "NEST" k                                 */
-virtual  r_8&       PixValNest(int_4 k);
-virtual  r_8 const&       PixValNest(int_4 k) const;
+virtual T& PixValNest(int_4 k);
+virtual T const& PixValNest(int_4 k) const;
 
 /* Indice "NEST" du pixel vers lequel pointe une direction definie par 
 ses  coordonnees spheriques */                                    
-virtual  int_4      PixIndexSphNest(float theta, float phi) const;
+virtual int_4 PixIndexSphNest(float theta, float phi) const;
 
 /* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
-virtual  void       PixThetaPhiNest(int_4 k, float& teta, float& phi) const;
+virtual void PixThetaPhiNest(int_4 k, float& teta, float& phi) const;
 
-   // algorithme de pixelisation
+/* algorithme de pixelisation */
 void Pixelize(int_4); 
 
-
 /* convertit index nested en ring  */
-int_4 NestToRing(int_4 );
-
+int_4 NestToRing(int_4 ) const;
 
 /* convertit index ring en nested" */
-int_4 RingToNest(int_4 );
-
+int_4 RingToNest(int_4 ) const;
 
 /*    analyse en harmoniques spheriques des valeurs des pixels de la 
    sphere : appel du module anafast (Gorski-Hivon) */
-void anharm(int, float, float*);
+//void anharm(int, float, float*);
  
 /*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) */
-void synharm(int , int, float, float* ); 
+//void synharm(int, int, float, float*); 
 
+/* retourne/fixe la valeur du parametre Gorski */
+inline virtual int_4 SizeIndex() const {return(nSide_);}
+inline void setSizeIndex(int_4 n) {nSide_= n;}
+
+/* retourne les pointeurs /remplit les tableaux */
+inline const NDataBlock<T>* getDataBlock() const { return (&pixels_); }
+inline void setDataBlock(T* data, int_4 m) { pixels_.FillFrom(m,data); }
+
+/* retourne/fixe les parametres des modules anafast et synfast */
+void getParafast(int_4& nlmax,int_4& nmmax,int_4& iseed,float& fwhm,float& quadr,float& cut) const;
+void setParafast(int_4 nlmax,int_4 nmmax,int_4 iseed,float fwhm,float quadr,float cut,char* filename);
+
+/* retourne/fixe le nom du fichier qui contient le spectre de puissance */
+void powfile(char filename[]) const;
+
+/* impression */ 
+void print(ostream& os) const;
+         
+private :
 
 // ------------- méthodes internes ----------------------
-          
-private :
+void InitNul();
 
-void          InitNul();
-void          Clear();
-int           nest2ring(int nside, int ipnest) const;
-void          mk_pix2xy(int *pix2x,int *pix2y);
-int           ring2nest(int nside, int ipring) const;
-void          mk_xy2pix(int *x2pix, int *y2pix);
-int           ang2pix_ring(int nside, double theta, double phi) const;
-int           ang2pix_nest(int nside, double theta, double phi) const;
-void          pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const;
-void          pix2ang_nest(int nside, int ipix, double& theta, double& phi) const;
+int  nest2ring(int nside, int ipnest) const;
+int  ring2nest(int nside, int ipring) const;
+
+int  ang2pix_ring(int nside, double theta, double phi) const;
+int  ang2pix_nest(int nside, double theta, double phi) const;
+void pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const;
+void pix2ang_nest(int nside, int ipix, double& theta, double& phi) const;
+
 // ------------- variables internes -----------------------
-                   private :
-class Pix2XY {
-  public :
-Pix2XY() { 
-    pix2x_=new int[1024];
-    pix2y_=new int[1024];
-    mk_pix2xy(pix2x_,pix2y_); 
-  }
-~Pix2XY() {
-  if (pix2x_)  delete[] pix2x_;
-  if (pix2y_)  delete[] pix2y_;
-  }
-void          pix2ang_nest(int nside, int ipix, double& theta, double& phi) const;
+int_4 nSide_;
+int_4 nPix_;
+r_8 omeg_;
 
-private : 
-int *pix2x_;
-int *pix2y_;
-void          mk_pix2xy(int *pix2x,int *pix2y);
-
-};
+NDataBlock<T> pixels_;
 
 int nlmax_;
 int nmmax_;
 int iseed_;
-int  nSide_;
-int_4  nPix_;
-int *pix2x_;
-int *pix2y_;
-int *x2pix_;
-int *y2pix_;
 float fwhm_;
 float quadrupole_;
 float sym_cut_deg_;
-r_8  omeg_;
-r_8*    mPix_;
 char powFile_[128];
-Pix2XY *pix2xy_;
+};
+
+//
+// ------------- Classe pour la gestion de persistance --
+//
+template <class T>
+class FIO_SphereGorski : public PPersist  
+{
+public:
+
+FIO_SphereGorski();
+FIO_SphereGorski(string const & filename); 
+FIO_SphereGorski(const SphereGorski<T>& obj);
+FIO_SphereGorski(SphereGorski<T>* obj);
+virtual ~FIO_SphereGorski();
+virtual AnyDataObj* DataObj();
+inline operator SphereGorski<T>() { return(*dobj); }
+inline SphereGorski<T> getObj() { return(*dobj); }
+
+protected :
+
+virtual void ReadSelf(PInPersist&);           
+virtual void WriteSelf(POutPersist&) const;  
+SphereGorski<T>* dobj;
+bool ownobj;
+};
+
+//
+// ------------- Classe PIXELS_XY -----------------------
+//
+class PIXELS_XY
+{
+
+public :
+
+static PIXELS_XY& instance();
+
+NDataBlock<int> pix2x_;
+NDataBlock<int> pix2y_;
+NDataBlock<int> x2pix_;
+NDataBlock<int> y2pix_;
+
+private :
+
+PIXELS_XY();
+void mk_pix2xy();
+void mk_xy2pix();
 };
 #endif

trunk/SophyaLib/Samba/spherethetaphi.cc

@@ -1 @@
-//
-// 
 #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 +61 @@
 //++
 
-
-SphereThetaPhi::SphereThetaPhi()
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi()
 
 //--
@@ -49 +71 @@
 /* --Methode-- */
 //++
-SphereThetaPhi::SphereThetaPhi(char* flnm)
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(char* flnm)
 
 //    Constructeur : charge une image à partir d'un fichier 
@@ -55 +78 @@
 {
   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 +96 @@
 {
   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 +124 @@
 //--
 //++
-SphereThetaPhi::~SphereThetaPhi()
+template <class T>
+SphereThetaPhi<T>::~SphereThetaPhi()
 
 //--
@@ -89 +135 @@
 // 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 +231 @@
   
   // 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 +253 @@
   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 +280 @@
 //--
 {
-  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 +295 @@
   bool fgzn= false;
   
-  if( (k< 0) || (k >= mNPix) ) {
+  if( (k< 0) || (k >= NPix_) ) {
     tetMin= -99999.; 
     phiMin= -99999.;  
@@ -322 +304 @@
   
   // 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 +338 @@
 /* --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 +345 @@
 {
   int_4 nbslices;
-  nbslices= 2*mNTheta;
+  nbslices= 2*NTheta_;
   return(nbslices);
 }
@@ -368 +351 @@
 /* --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 +374 @@
 /* --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 +389 @@
 
   // 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 +406 @@
 /* --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 +420 @@
   
   // 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 +429 @@
   }
   
-  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 +437 @@
 /* --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 +470 @@
 /* --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 +512 @@
 {
   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 +518 @@
   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 +537 @@
   
   //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 +586 @@
 {
   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);
+}

trunk/SophyaLib/Samba/spherethetaphi.h

@@ -3 +3 @@
 
 #include "sphericalmap.h"
-#include "pdataarray.h"
-#include "cvector.h"
+#include "ndatablock.h"
+#include "tvector.h"
 
-// ***************** CLASSE SphereThetaPhi *****************************
+#include "anydataobj.h"
+#include "ppersist.h"
 
+// ***************** Class SphereThetaPhi *****************************
 
-class SphereThetaPhi : public SphericalMap {
+template <class T>
+class SphereThetaPhi : public SphericalMap<T>, public AnyDataObj
+{
+
 public :
 
-
-SphereThetaPhi(int_4 m, int_4 pet);
 SphereThetaPhi();
+SphereThetaPhi(int_4 m);
 SphereThetaPhi(char* flnm);
+SphereThetaPhi(const SphereThetaPhi<T>& s);
 virtual ~SphereThetaPhi();
 
-// ------------ Persistence handling
+// ------------ Definition of PixelMap abstract methods -
 
-     enum {classId = 0xF001 };
-int_4                  ClassId() const        { return classId; }
+/* retourne/fixe le nombre de pixels */ 
+virtual int_4 NbPixels() const;
+inline void setNbPixels(int_4 nbpix) { NPix_= nbpix; }
 
-virtual void           WriteSelf(POutPersist&) const;
-virtual void           ReadSelf(PInPersist&);
+/* retourne la valeur du pixel d'indice k */ 
+virtual T&       PixVal(int_4 k);
+virtual T const& PixVal(int_4 k) const;
 
-// ------------------ Definition of PixelMap abstract methods
+/* retourne l'indice du pixel a (theta,phi) */ 
+virtual int_4 PixIndexSph(float theta, float phi) const;
 
-// number of pixels
-virtual  int_4           NbPixels() const;
+/* retourne les coordonnees Spheriques du centre du pixel d'indice k */ 
+virtual void PixThetaPhi(int_4 k, float& theta, float& phi) const;
 
-// Value of pixel number k
-virtual  r_8&       PixVal(int_4 k);
-virtual  r_8 const&       PixVal(int_4 k) const;
+/* retourne/fixe l'angle Solide de Pixel   (steradians) */ 
+virtual r_8 PixSolAngle(int_4 dummy) const;
+inline void setPixSolAngle(r_8 omega) { Omega_= omega; }
+ 
+/* retourne/fixe la valeur du parametre de decoupage m */ 
+inline virtual int_4 SizeIndex() const { return( NTheta_); }
+inline void setSizeIndex(int_4 nbindex) { NTheta_= nbindex; }
 
-// Index of pixel at (theta,phi)
-virtual  int_4           PixIndexSph(float theta, float phi) const;
+// ------------- Specific methods  ----------------------
 
-// Spherical coordinates of center of pixel number k
-virtual  void          PixThetaPhi(int_4 k, float& theta, float& phi) const;
+virtual void Resize(int_4 m);
 
-// Pixel Solid angle  (steradians)
-virtual r_8         PixSolAngle(int_4 dummy) const; 
-// --------------- Specific methods 
+inline virtual char* TypeOfMap() const {return "TETAFI";};
 
+/* Valeurs de theta des paralleles et phi des meridiens limitant  le pixel d'indice k */
+virtual void Limits(int_4 k,float& tet1,float& tet2,float& phi1,float& phi2);
 
-/* Valeurs de theta des paralleles et phi des meridiens limitant          */
-/* le pixel d'indice k                                                   */
-virtual  void          Limits(int_4 k, float& tet1, float& tet2, 
-                                float& phi1, float& phi2 );
+/* Nombre de tranches en theta */ 
+int_4 NbThetaSlices() const;
 
-/* Nombre de tranches en theta                                           */
-int_4           NbThetaSlices() const;
-
-/* Nombre de pixels en phi de la tranche d'indice kt                     */
-int_4           NPhi(int_4 kt) const;
+/* Nombre de pixels en phi de la tranche d'indice kt */ 
+int_4 NPhi(int_4 kt) const;
 
 /* Renvoie dans t1,t2 les valeurs respectives de theta min et theta max  */
-/* de la tranche d'indice kt                                             */
-void          Theta(int_4 kt, float& t1, float& t2);
+/* de la tranche d'indice kt  */ 
+void Theta(int_4 kt, float& t1, float& t2);
 
 /* Renvoie dans p1,p2 les valeurs phimin et phimax du pixel d'indice jp  */
-/* dans la tranche d'indice kt                                           */
-void          Phi(int_4 kt, int_4 jp, float& p1, float& p2);
+/* dans la tranche d'indice kt  */ 
+void Phi(int_4 kt, int_4 jp, float& p1, float& p2);
 
 /* Renvoie l'indice k du pixel d'indice jp dans la tranche d'indice kt   */
-int_4           Index(int_4 kt, int_4 jp) const;
+int_4 Index(int_4 kt, int_4 jp) const;
 
+/* Indice kt de la tranche et indice jp du pixel d'indice k  */ 
+void ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp);
 
-/* Indice kt de la tranche et indice jp du pixel d'indice k            */
-void          ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp);
+void Pixelize(int_4); 
 
-void Pixelize(int_4,int_4); 
-void GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const; 
+void GetThetaSlice(int_4 index,r_4& theta,TVector<float>& phi,TVector<T>& value) const; 
 
+/*retourne le tableau contenant le nombre de pixels en phi de chacune des bandes en theta */
+inline const int_4* getmNPhi() const { return(NPhi_); }
+void setmNPhi(int_4* array, int_4 m);
+
+/* retourne/remplit le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes */
+inline const int_4* getmTNphi() const { return(TNphi_); }
+void setmTNphi(int_4* array, int_4 m);
+
+/* retourne/remplit le tableau contenant les valeurs limites de theta */
+inline const r_4* getmTheta() const { return(Theta_); }
+void setmTheta(r_4* array, int_4 m);
+
+/* retourne le pointeur vers/remplit  le vecteur des contenus des pixels */ 
+inline const NDataBlock<T>* getDataBlock() const { return (&pixels_); }
+void setDataBlock(T* data, int_4 m);
+
+/* impression */ 
+void print(ostream& os) const;
 
 private :
-// ------------- méthodes internes ----------------------
-          
 
-void          InitNul();
-void          Clear();
+// ------------- méthodes internes ----------------------          
+void InitNul();
+void Clear();
 
-// ------------- variables internes -----------------------
-int_4 mNTheta,mNPet;
-int_4 mNPix;
-r_4* mTheta;
-r_8  mOmeg;
-int_4* mNphi;
-int_4* mTNphi;
-//r_8* mPix;
-PDataArray<r_8>* _pixel;
-
-
+// ------------- variables internes ---------------------
+int_4 NTheta_;
+int_4 NPix_;
+r_8 Omega_;
+int_4* NPhi_;
+int_4* TNphi_;
+r_4* Theta_;
+NDataBlock<T> pixels_;
 };
 
+// ------------- Classe pour la gestion de persistance --
+template <class T>
+class FIO_SphereThetaPhi : public PPersist  
+{
+public:
 
+FIO_SphereThetaPhi();
+FIO_SphereThetaPhi(string const & filename); 
+FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj);
+FIO_SphereThetaPhi(SphereThetaPhi<T>* obj);
+virtual ~FIO_SphereThetaPhi();
+virtual AnyDataObj* DataObj();
+inline operator SphereThetaPhi<T>() { return(*dobj); }
+inline SphereThetaPhi<T> getObj() { return(*dobj); }
+
+protected :
+
+virtual void ReadSelf(PInPersist&);           
+virtual void WriteSelf(POutPersist&) const;  
+SphereThetaPhi<T>* dobj;
+bool ownobj;
+};
 
 #endif

trunk/SophyaLib/Samba/sphericalmap.h

@@ -6 +6 @@
 #include <math.h>
 #include "pixelmap.h"
-#include "cvector.h" 
+#include "tvector.h" 
+
 // Map of pixels on a whole sphere.
 // Class hierarchy :
@@ -16 +17 @@
 //      LocalMap
 
-class SphericalMap : public PixelMap {
+template<class T>
+class SphericalMap : public PixelMap<T> 
+{
+  
+public :
 
-public :
 SphericalMap() {};
 virtual  ~SphericalMap() {};
 
-   // Overloading of () to access pixel number k.
-   inline  r_8&        operator()(int_4 k)
-                          {return(PixVal(k));}
-   inline  r_8 const&  operator()(int_4 k) const
-                          {return(PixVal(k));}
-inline  r_8&        operator()(float theta, float phi)
-                { return(PixValSph(theta, phi)) ; };
-inline  r_8 const&        operator()(float theta, float phi) const
-                { return(PixValSph(theta, phi)) ; };
+// Overloading of () to access pixel number k.
+inline T& operator()(int_4 k) {return(PixVal(k));}
+inline T  const& operator()(int_4 k) const {return(PixVal(k));}
+inline T& operator()(float theta,float phi) {return(PixValSph(theta,phi));};
+inline T  const& operator()(float theta,float phi) const {return(PixValSph(theta,phi));};
 
-
+// index characterizing the size pixelization : m for SphereThetaPhi
+// nside for Gorski sphere...
+virtual void Resize(int_4 m)=0;
 virtual int_4 NbThetaSlices() const=0;
-virtual void GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const=0; 
+virtual void GetThetaSlice(int_4 index, r_4& theta, TVector<float>& phi, TVector<T>& value) const=0; 
 };
-
 #endif