Index: trunk/SophyaLib/Samba/alm.h
===================================================================
--- trunk/SophyaLib/Samba/alm.h	(revision 757)
+++ trunk/SophyaLib/Samba/alm.h	(revision 758)
@@ -1,59 +1,9 @@
 #ifndef ALM_SEEN
 #define ALM_SEEN
+
 #include "nbrandom.h"
 #include "nbmath.h"
+#include "triangmtx.h"
 #include "tvector.h"
-#include "pexceptions.h"
-/*! Class for inferior triangular matrix (base class for the class Alm) */
-template <class T>
-class TriangularMatrix 
-  {
-    
-  public :
-
-TriangularMatrix()   {};
-  /* instanciate a triangular matrix from the number of rows */
-TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) {elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); };
-TriangularMatrix(const TriangularMatrix<T>& a,  bool share=false)  : elem_(a.elem_, share),  long_diag_(a.long_diag_) {;}
-/*! resize the matrix with a new number of rows */
-inline void ReSizeRow(int rowSize) 
-  {
-    long_diag_=(uint_4)rowSize;
-    elem_.ReSize(long_diag_*(long_diag_+1)/2);
-  }
- inline void SetTemp(bool temp=false) const {elem_.SetTemp(temp);}
-
-inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
-{
-  elem_=a.elem_;
-  long_diag_ = a.long_diag_;
-  return *this;
-}
-inline T& operator()(int l, int m) 
-  {
-      return  elem_(adr_ij(l,m));
-  }
-inline T const& operator()(int l, int m) const 
-  {
-      return *(elem_.Begin()+ adr_ij(l,m));
-  }
- inline  int_4  rowNumber() const {return (int_4)long_diag_;}
-  private: 
- /*! compute the address of an element in the single array representing the matrix */
-inline uint_4 adr_ij(int i,int j) const 
-{
-  int adr= i*(i+1)/2+j;
-  if ( adr >= elem_.Size() || adr <0 )
-{
- cout << " attention depassement dans triangularMatrix " << endl;
- cout << " l= " << i << " m= " << j << " tableau reserve longueur " << elem_.Size() << endl; 
-}
-  return(i*(i+1)/2+j);
-}
-
- uint_4 long_diag_;
- NDataBlock<T> elem_;
-
-  };
 
 
Index: trunk/SophyaLib/Samba/localmap.cc
===================================================================
--- trunk/SophyaLib/Samba/localmap.cc	(revision 757)
+++ 	(revision )
@@ -1,755 +1,0 @@
-#include "localmap.h"
-#include "nbmath.h"
-#include <complex>
-#include "piocmplx.h"
-
-#include <string.h>
-#include <iostream.h>
-
-
-//*****************************************************************************
-//++ 
-// Class	LocalMap
-//
-// include      localmap.h nbmath.h
-//
-//    A local map of a region of the sky, in cartesian coordinates.
-//    It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
-//    (x0, y0 might be outside of this local map)
-//    default value of (x0, y0) is middle of the map, center of 
-//    pixel(nx/2, ny/2)
-//
-//    A local map is a 2 dimensional array, with i as column index and j
-//    as row index. The map is supposed to lie on a plan tangent to the 
-//    celestial sphere in a point whose coordinates are (x0,y0) on the local
-//    map and (theta0, phi0) on the sphere. The range of the map is defined
-//    by two values of angles covered respectively by all the pixels in 
-//    x direction and all the pixels in y direction (SetSize()).
-//
-//    A "reference plane" is considered : this plane is tangent to the 
-//    celestial sphere in a point with angles theta=Pi/2 and phi=0. This 
-//    point is the origine of coordinates is of the reference plane. The 
-//    x-axis is the tangent parallel to the equatorial line and oriented 
-//    toward the increasing phi's ; the y-axis is parallel to the meridian 
-//    line and oriented toward the north pole.
-//    
-//    Internally, a map is first defined within this reference plane and
-//    tranported until the point (theta0, phi0) in such a way that both 
-//    axes are kept parallel to meridian and parallel lines of the sphere.
-//    The user can define its own map with axes rotated with respect to
-//    reference axes (this rotation is characterized by angle between 
-//    the local parallel line and the wanted x-axis-- see method 
-//    SetOrigin(...))
-//    
-//
-
-//
-//-- 
-//++
-//
-// Links	Parents
-//
-//    PixelMap
-//    
-//-- 
-//++
-//
-// Links	Childs
-//
-//     
-//-- 
-//++
-// Titre	Constructors
-//--
-//++
-template<class T>
-LocalMap<T>::LocalMap()
-//
-//     
-//--
-{
-  InitNul();
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny)
-//
-//    
-//--
-{
-  InitNul();
-  nPix_= nx*ny;
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)
-  : pixels_(lm.pixels_, share)
-
-//
-//    copy constructor
-//--
-{
-  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_;
-}
-
-//++
-// Titre	Destructor
-//--
-//++
-template<class T>
-LocalMap<T>::~LocalMap()
-//
-//--
-{
-  InitNul();
-}
-
-
-
-//++
-// Titre	Public Methods
-//--
-
-//++
-template<class T>
-void LocalMap<T>::ReSize(int_4 nx, int_4 ny) 
-//
-//    Resize storage area for pixels 
-//--
-{
-  InitNul();
-  nSzX_ = nx;
-  nSzY_ = ny;
-  nPix_= nx*ny;
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-int_4 LocalMap<T>::NbPixels() const
-//
-//    Return number of pixels
-//--
-{
-  return(nPix_);
-}
-
-//++
-template<class T>
-T& LocalMap<T>::PixVal(int_4 k)
-//
-//    Return value of pixel with index k
-//--
-{
-  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
-//
-//   const version of previous method
-//--
-{
-  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>
-bool LocalMap<T>::ContainsSph(double theta, double phi) const
-{
-// $CHECK$  Reza 11/11/99 - A modifier
-return(true);
-}
-
-//++
-template<class T>
-int_4 LocalMap<T>::PixIndexSph(double theta,double phi) const
-//
-//    Return index of the pixel with spherical coordinates (theta,phi)
-//--
-{
-  int_4 i,j;
-  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)
-  double theta_aux= theta;
-  double phi_aux  = phi;
-  UserToReference(theta_aux, phi_aux);
-
-  // coordonnees dans le plan local en unites de pixels
-  double x,y;
-  AngleProjToPix(theta_aux,phi_aux, x, y);
-
-  double xmin= -x0_-0.5;
-  double xmax= xmin+nSzX_;
-  if((x > xmax) || (x < xmin)) return(-1);
-  double xcurrent= xmin;
-  for(i = 0; i < nSzX_; i++ ) 
-    {
-      xcurrent += 1.;
-      if( x < xcurrent ) break;
-    }
-  double ymin= -y0_-0.5;
-  double ymax= ymin+nSzY_;
-  if((y >  ymax) || (y < ymin)) return(-1);
-  double 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,double& theta,double& phi) const
-//
-//    Return (theta, phi) coordinates of pixel with index k
-//--
-{
-  if(!(originFlag_) || !(extensFlag_)) 
-    {
-      cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
-      exit(0);
-    }
-
-  int_4 i,j;
-  Getij(k,i,j);
-
-  double X= double(i-x0_);
-  double Y= double(j-y0_);
-  // situation de ce pixel dans le plan de reference
-  double x= X*cos_angle_-Y*sin_angle_;
-  double y= X*sin_angle_+Y* cos_angle_;
-  // projection sur la sphere
-  PixProjToAngle(x, y, theta, phi);
-  // retour au plan utilisateur
-  ReferenceToUser(theta, phi);
-}
-
-template <class T>
-T LocalMap<T>::SetPixels(T v)
-{
-pixels_.Reset(v);
-return(v);
-}
- 
-//++
-template<class T>
-double 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_4 i,j;
-  Getij(k,i,j);
-  double X= double(i-x0_);
-  double Y= double(j-y0_);
-  double XR= X+double(i)*0.5;
-  double XL= X-double(i)*0.5;
-  double YU= Y+double(j)*0.5;
-  double YL= Y-double(j)*0.5;
-
-  // situation  dans le plan de reference
-  double x0= XL*cos_angle_-YL*sin_angle_;
-  double y0= XL*sin_angle_+YL*cos_angle_;
-  double xa= XR*cos_angle_-YL*sin_angle_;
-  double ya= XR*sin_angle_+YL*cos_angle_;
-  double xb= XL*cos_angle_-YU*sin_angle_;
-  double yb= XL*sin_angle_+YU*cos_angle_;
-
-  // projection sur la sphere
-  double thet0,phi0,theta,phia,thetb,phib;
-  PixProjToAngle(x0, y0, thet0, phi0);
-  PixProjToAngle(xa, ya, theta, phia);
-  PixProjToAngle(xb, yb, thetb, phib);
-
-  // angle solide
-  double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
-  return sol;
-}
-
-//++
-template<class T>
-void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)
-//
-//    set the referential of the map (angles in degrees)
-//    (default x0=siz_x/2,  y0=siz_y/2) 
-//--
-{
-  theta0_= theta0;
-  phi0_  = phi0;
-  angle_ = angle;
-  x0_= nSzX_/2;
-  y0_= nSzY_/2;
-  cos_angle_= cos(angle*Pi/180.);
-  sin_angle_= sin(angle*Pi/180.);
-  originFlag_= true; 
-  cout << " LocalMap:: set origin 1 done" << endl;
-}
- 
-//++
-template<class T>
-void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle)
-//
-//    set the referential of the map (angles in degrees)
-//--
-{
-  theta0_= theta0;
-  phi0_  = phi0;
-  angle_ = angle;
-  x0_= x0;
-  y0_= y0;
-  cos_angle_= cos(angle*Pi/180.);
-  sin_angle_= sin(angle*Pi/180.);
-  originFlag_= true;
-  cout << " LocalMap:: set origin 2 done" << endl;
-} 
-
-//++
-template<class T>
-void LocalMap<T>::SetSize(double angleX,double angleY)
-//
-//    angle range of tthe map (angles in degrees)
-//--
-{
-  angleX_= angleX;
-  angleY_= angleY;
-
-  // tangente de la moitie de l'ouverture angulaire totale
-  tgAngleX_= tan(0.5*angleX_*Pi/180.);
-  tgAngleY_= tan(0.5*angleY_*Pi/180.);
-
-  extensFlag_= true;
-  cout << " LocalMap:: set extension done" << endl;
-}
-
-//++
-template<class T>
-void LocalMap<T>::Project(SphericalMap<T>& sphere) const
-//
-//    Projection to a spherical map   
-//--
-{
-  for(int m = 0; m < nPix_; m++) 
-    {
-      double theta,phi;
-      PixThetaPhi(m,theta,phi);
-      sphere(theta,phi)= pixels_(m);
-      // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
-    }
-}
-// Titre	Private Methods
-//++
-template<class T>
-void LocalMap<T>::InitNul()
-//
-//    set some attributes to zero
-//--
-{ 
-  originFlag_= false;
-  extensFlag_= false;
-  cos_angle_= 1.0;
-  sin_angle_= 0.0;
-//  pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
-}
-
-//++
-template<class T>
-void LocalMap<T>::Getij(int_4 k,int_4& i,int_4& j) const
-//
-//    Return 2 indices corresponding to the pixel number k
-//--
-{
-  i= (k+1)%nSzX_-1;
-  if(i == -1) i= nSzX_-1;
-  j= (k-i+2)/nSzX_;
-}
-
-//++
-template<class T>
-void  LocalMap<T>::ReferenceToUser(double& theta,double& phi) const
-//
-//    Transform a pair of coordinates (theta, phi) given in 
-//    reference coordinates into map coordinates
-//--     
-{
-  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
-    {
-      throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
-    } 
-
-  theta= theta0_*Pi/180.+theta-Pi*0.5;
-  if(theta < 0.) 
-    {
-      theta= -theta;
-      phi += Pi;
-    } 
-  else 
-    {
-      if(theta > Pi) 
-	{
-	  theta= 2.*Pi-theta;
-	  phi += Pi;
-	} 
-    }
-
-  phi= phi0_*Pi/180.+phi;
-  while(phi >= 2.*Pi) phi-= 2.*Pi;
-
-  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 
-    {
-      cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
-      cout << " theta= " << theta << " phi= " << phi << endl;
-      exit(0);
-    }
-}
-
-//++
-template<class T>
-void  LocalMap<T>::UserToReference(double& theta,double& phi) const
-//
-//    Transform a pair of coordinates (theta, phi) given in 
-//    map coordinates into reference coordinates
-//--     
-{
-  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*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"; 
-    }
-
-  double phi1= phi-phi0_*Pi/180.;
-  if(phi1 < 0.) phi1+= 2.*Pi;
-
-  double theta1= theta-theta0_*Pi/180.+Pi*0.5;
-  if(theta1 < 0.) 
-    {
-      theta= -theta1;
-      phi1+= Pi;
-    } 
-  else
-    { 
-      if(theta1 > Pi) 
-	{
-	  theta= 2.*Pi-theta1;
-	  phi1+= Pi;
-	}
-    }
-
-  while(phi1 >= 2.*Pi) phi1-= 2.*Pi;
-  phi= phi1;
-  if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*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(double x,double y,double& theta,double& phi) const
-//
-//  
-//    Given coordinates in pixel units in the REFERENCE PLANE, return
-//    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
-//--
-{
-  theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);
-  phi= atan2(2.*x*tgAngleX_,(double)nSzX_);
-  if(phi < 0.) phi += DeuxPi;
-}
-
-//++
-template<class T>
-void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const
-//
-//    Given coordinates  (theta, phi) in "absolute" referential 
-//    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
-//--
-{
-  if(phi >= Pi) phi-= DeuxPi;
-  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
-  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;
-}
-//++
-// Titre	class FIO_LocalMap
-//    Delegated objects for persitance management
-//--
-
-//*******************************************************************
-// 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>;
-  dobj->DataBlock().SetTemp(true);
-  ownobj= true;
-  Read(filename);
-}
-
-//++
-template <class T>
-FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)
-//
-//--
-{
-  dobj= new LocalMap<T>(obj, true);
-  dobj->DataBlock().SetTemp(true);
-  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)
-//
-//--
-{
-
-  if(dobj == NULL) 
-    {
-      dobj= new LocalMap<T>;
-      dobj->DataBlock().SetTemp(true);
-      ownobj= true;
-    }
-
-  // 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);
-  dobj->ReSize(nSzX, nSzY); 
-  string ss("local mapping is done");
-  string sso;
-  is.GetStr(sso);
-  if(sso == ss)
-    {
-      cout<<" ReadSelf:: local mapping"<<endl;
-      int_4 x0, y0;
-      double theta, phi, angle;
-      is.GetI4(x0);
-      is.GetI4(y0);
-      is.GetR8(theta);
-      is.GetR8(phi);
-      is.GetR8(angle);
-      dobj->SetOrigin(theta, phi, x0, y0, angle);
-
-      double angleX, angleY;
-      is.GetR8(angleX);
-      is.GetR8(angleY);
-      dobj->SetSize(angleX, angleY);
-    }
-
-// On lit le DataBlock;
-  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
-  fio_nd.Read(is);
-}
-
-//++
-template <class T>
-void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const
-//
-//--
-{
-  if(dobj == NULL) 
-    {
-      cout << " FIO_LocalMap::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;
-      double theta, phi, angle;
-      dobj->Origin(theta, phi, x0, y0, angle);
-      os.PutI4(x0);
-      os.PutI4(y0);
-      os.PutR8(theta);
-      os.PutR8(phi);
-      os.PutR8(angle);
-
-      double angleX, angleY;
-      dobj->Aperture(angleX, angleY);
-      os.PutR8(angleX);
-      os.PutR8(angleY);
-    }
-  else
-    {
-      string ss("no local mapping");
-      os.PutStr(ss);
-    }
-
-// On ecrit le dataBlock 
-  FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
-  fio_nd.Write(os);
-}
-
-#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
Index: trunk/SophyaLib/Samba/localmap.h
===================================================================
--- trunk/SophyaLib/Samba/localmap.h	(revision 757)
+++ 	(revision )
@@ -1,233 +1,0 @@
-// This may look like C code, but it is really -*- C++ -*-
-#ifndef LOCALMAP_SEEN
-#define LOCALMAP_SEEN
-
-#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.
-/*! A local map has an origin in (theta0, phi0), mapped to pixel(x0, y0)
-   (x0, y0 might be outside of this local map)
- default value of (x0, y0) is middle of the map, center of pixel(nx/2, ny/2)
-    A local map is a 2 dimensional array, with i as column index and j
-    as row index. The map is supposed to lie on a plan tangent to the 
-    celestial sphere in a point whose coordinates are (x0,y0) on the local
-   map and (theta0, phi0) on the sphere. The range of the map is defined
-    by two values of angles covered respectively by all the pixels in 
-    x direction and all the pixels in y direction (SetSize()).
-
-    A "reference plane" is considered : this plane is tangent to the 
-    celestial sphere in a point with angles theta=Pi/2 and phi=0. This 
-    point is the origine of coordinates is of the reference plane. The 
-    x-axis is the tangent parallel to the equatorial line and oriented 
-    toward the increasing phi's ; the y-axis is parallel to the meridian 
-    line and oriented toward the north pole.
-    
-    Internally, a map is first defined within this reference plane and
-    tranported until the point (theta0, phi0) in such a way that both 
-    axes are kept parallel to meridian and parallel lines of the sphere.
-    The user can define its own map with axes rotated with respect to
-    reference axes (this rotation is characterized by angle between 
-    the local parallel line and the wanted x-axis-- see method 
-    SetOrigin(...))
-*/
-//
-//    la carte est consideree comme un tableau a deux indices i et j, i etant 
-//    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 
-//    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 
-//    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 
-//    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 -- 
-//    methode SetOrigin(...))
-
-
-// ***************** Class LocalMap *****************************
-
-template<class T>
-class LocalMap : public PixelMap<T>
-{
-
-public:
-
-LocalMap();
-LocalMap(int_4 nx, int_4 ny);
-LocalMap(const LocalMap<T>& lm, bool share=false);
-virtual ~LocalMap();
-
-/*! Setting blockdata to temporary (see ndatablock documentation) */
-inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};
-
-
-// ---------- 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()(int_4 ix, int_4 iy) {return PixVal(iy*nSzX_+ix);};
-inline T const& operator()(int_4 ix, int_4 iy) const {return PixVal(iy*nSzX_+ix);};
-   
-// ---------- Definition of PixelMap abstract methods -------
-
-/* return/set the number of pixels */
-/*!    Return number of pixels */
-virtual int_4 NbPixels() const;   // D.Y. int change en int_4 rationalisation Mac
-  
-/* return the value of pixel number k */
-/*!    Return value of pixel with index k */
-virtual T& PixVal(int_4 k);
-/*!   const version of previous method */
-virtual T const& PixVal(int_4 k) const;
-
-/* Return true if teta,phi in map  */
-virtual bool ContainsSph(double theta, double phi) const;
-/* return the index of pixel at (theta,phi) */
-/*!    Return index of the pixel with spherical coordinates (theta,phi) */
-virtual int_4 PixIndexSph(double theta,double phi) const;
-   
-/* return the spherical coordinates of center of pixel number k */
-/*!    Return (theta, phi) coordinates of pixel with index k */
-virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;
-
-/*! Set all pixels to value v */
-virtual T SetPixels(T v);
-
-/* return the Pixel Solid angle  (steradians) */
-/*!    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. 
-*/  
-virtual double PixSolAngle(int_4 k) const; 
-
-// ---------- Specific methods ------------------------------
-
-/*!    Resize storage area for pixels */
-void ReSize(int_4 nx, int_4 ny);
-
-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*/
-/*!    set the referential of the map (angles in degrees)
-
-    (default x0=siz_x/2,  y0=siz_y/2) 
-*/
-virtual void SetOrigin(double theta=90.,double phi=0.,double angle=0.);
-/*!    set the referential of the map (angles in degrees) */
-virtual void SetOrigin(double theta,double phi,int_4 x0,int_4 y0,double angle=0.);
-
-/* Pixel size (degres) */ 
-/*!    angle range of tthe map (angles in degrees) */
-virtual void SetSize(double angleX,double angleY);
-
-/* Check to see if the local mapping is done */
-inline bool LocalMap_isDone() const {return(originFlag_ && extensFlag_);};
-
-/*! Projection to a 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 */
-  
-/* 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 int_4 XSize() const {return nSzX_;}
-inline int_4 Size_y() const {return nSzY_;}
-inline int_4 YSize() const {return nSzY_;}
-
-inline void Origin(double& theta,double& phi,int_4& x0,int_4& y0,double& angle) const {theta= theta0_; phi= phi0_; x0= x0_; y0= y0_;angle= angle_;}
-
-inline void Aperture(double& anglex,double& angley) const {anglex= angleX_; angley= angleY_;}
-
-
-/*  Acces to the DataBlock  */
-inline       NDataBlock<T>& DataBlock()       {return pixels_;}
-inline const NDataBlock<T>& DataBlock() const {return pixels_;}
-
-/* impression */ 
-void print(ostream& os) const;
-
-// ---------- Méthodes internes -----------------------------
-          
-private :
-
-void InitNul();
-/*!    Return 2 indices corresponding to the pixel number k */
-void Getij(int_4 k,int_4& i,int_4& j) const;
-/*!    Transform a pair of coordinates (theta, phi) given in 
-    reference coordinates into map coordinates
-*/
-void ReferenceToUser(double& theta,double& phi) const; 
-/*!    Transform a pair of coordinates (theta, phi) given in 
-   map coordinates into reference coordinates
-*/
-void UserToReference(double& theta,double& phi) const;
-/*!   Given coordinates in pixel units in the REFERENCE PLANE, return
-    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
-*/
-void PixProjToAngle(double x,double y,double& theta,double& phi) const; 
-/*!    Given coordinates  (theta, phi) in "absolute" referential 
-    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
-*/
-void AngleProjToPix(double theta,double phi,double& x,double& y) const; 
-
-// ---------- Variables internes ----------------------------
-
-int_4 nSzX_;
-int_4 nSzY_;
-int_4 nPix_;
-bool originFlag_;
-bool extensFlag_;
-int_4 x0_;
-int_4 y0_;
-double theta0_;
-double phi0_;
-double angle_;
-double cos_angle_;
-double sin_angle_;
-double angleX_;
-double angleY_;
-double tgAngleX_;
-double tgAngleY_;
-NDataBlock<T> pixels_;
-};
-
-// ------------- Classe pour la gestion de persistance --
-//!    Delegated objects for persitance management
-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;
-};
-
-#endif
Index: trunk/SophyaLib/Samba/longlat.cc
===================================================================
--- trunk/SophyaLib/Samba/longlat.cc	(revision 757)
+++ 	(revision )
@@ -1,82 +1,0 @@
-#include "longlat.h"
-//++ 
-// Class	LongLat
-//
-// include      longlat.h
-//
-//    translation from  coordinates (longitude, latitude) 
-//    to theta,phi
-//
-//    longitude=phi
-//
-//    latitude=Pi/2-theta
-//
-//--
-//++
-// Titre	Constructors
-//--
-//++
-LongLat::LongLat()
-//
-//--
-{
-  _lon=0.;
-  _lat=0.;
-}
-//++
-LongLat::LongLat(double longitude, double latitude)
-//
-//--
-{
-  _lon=mod(longitude,pi2);
-  if( latitude == pi_over_2 ) _lat=latitude;
-  else _lat=mod(latitude+pi_over_2,M_PI)-pi_over_2; // dans [-Pi/2,Pi/2]
-}
-//++
-LongLat::LongLat(double x, double y, double z)
-//
-//--
-{
-  double norm=sqrt(x*x+y*y+z*z);
-  double theta;
-  if( norm != 0. ) 
-    {
-      theta=acos(z/norm); // dans [0,Pi]
-      if( mod(theta,M_PI) == 0. ) _lon=0.; // on est sur +-Oz, le vecteur z est en phi=0
-      //      else _lon=acos(x/sin(theta)/norm)+M_PI*(y<0);
-      else _lon=scangle(y/sin(theta)/norm,x/sin(theta)/norm);
-    }
-  else // vecteur nul
-    {
-      theta=0.;
-      _lon=0.;
-    }
-  _lat=pi_over_2-theta;
-}
-//++
-// Titre	Public Methods
-//--
-//++
-void LongLat::Set(double longitude, double latitude)
-//
-//--
-{
-  _lon=mod(longitude,pi2);
-  if( latitude == pi_over_2 ) _lat=latitude;
-  else _lat=mod(latitude+pi_over_2,M_PI)-pi_over_2; // dans [-Pi/2,Pi/2]
-}
-//++
-void LongLat::Print(ostream& os) const
-//
-//--
-{
-  os << "LongLat : longitude = " << _lon << " phi   = " << this->Phi() << endl;
-  os << "LongLat : latitude    = " << _lat << " theta = " << this->Theta() << endl;
-}
-//++
-//
-//  double Longitude() const 
-//  double Phi() const 
-//  double Latitude() 
-//  double Theta() const 
-//--
Index: trunk/SophyaLib/Samba/longlat.h
===================================================================
--- trunk/SophyaLib/Samba/longlat.h	(revision 757)
+++ 	(revision )
@@ -1,47 +1,0 @@
-#ifndef LONGLAT_H_SEEN
-#define LONGLAT_H_SEEN
-
-#include <math.h>
-#include <iostream.h>
-#include <stdio.h>
-#include <string.h>
-#include "utilgeom.h"
-
-/*
-  passage des coordonnees longitude, latitude en theta,phi
-  longitude=phi
-  latitude=Pi/2-theta
-*/
-
-class LongLat
-{
-  
- public:
-  
-  LongLat();
-  LongLat(double longitude, double latitude);
-  LongLat(double x, double y, double z);
-  LongLat(const LongLat&);
-  
-  void Set(double longitude, double latitude);
-  double Longitude() const {return _lon;}
-  double Phi() const {return _lon;}
-  double Latitude() const {return _lat;}
-  double Theta() const {return pi_over_2-_lat;}
-
-  void Print(ostream& os) const;
-  
- private:
-  
-  double _lon; // lambda
-  double _lat; // beta
-  
-};
-
-inline ostream& operator<<(ostream& s, const LongLat& ll) 
-{  
-  ll.Print(s);  
-  return s;  
-}
-
-#endif
Index: trunk/SophyaLib/Samba/pixelmap.h
===================================================================
--- trunk/SophyaLib/Samba/pixelmap.h	(revision 757)
+++ 	(revision )
@@ -1,138 +1,0 @@
-// This may look like C code, but it is really -*- C++ -*-
-#ifndef PIXELMAP_SEEN
-#define PIXELMAP_SEEN
-
-#include "ppersist.h" 
-#include "dvlist.h"
-#include "spherepos.h"
-#include <iostream.h>
-
-//! General map of pixels on part of sphere or whole sphere
-/*! Class hierarchy :
-\verbatim
-   PixelMap
-     SphericalMap
-         SphereThetaPhi
-         SphereGorski
-         SphereIco
-      LocalMap
-\endverbatim
-*/
-template<class T>
-class PixelMap  : public AnyDataObj 
-{
-
-public:
-  
-              PixelMap(SphereCoordSys* cs = NULL) : mInfo_(NULL) 
-                    { if (cs) cs_ = cs; else cs_ = new SphereCoordSys; }
-virtual       ~PixelMap() 
-                    { if (mInfo_) delete mInfo_; if (cs_) delete cs_; }
-
-// Set/Change/Get the coordinate system
-virtual void  SetCoordSys(SphereCoordSys* cs)
-                    { if (cs) { delete cs_;  cs_ = cs; } }
-inline SphereCoordSys* GetCoordSys() const { return(cs_); }
-
-/*! 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;
-
-// Map s coverage
-virtual bool ContainsSph(double theta, double phi) const=0;
-virtual bool Contains(const SpherePosition& spos) const;
-
-/*! Index of pixel at (theta,phi) */
-virtual int_4 PixIndexSph(double theta, double phi) const=0;
-// Index of pixel at a sky-position 
-virtual int_4 PixIndex(const SpherePosition& spos);
-
-/*! Value of pixel number at (theta,phi) */
-virtual T& PixValSph(double theta, double phi)
-                           {return PixVal(PixIndexSph(theta,phi));}
-virtual T const& PixValSph(double theta, double phi) const
-                           {return PixVal(PixIndexSph(theta,phi));}
-
-/*! Spherical coordinates of center of pixel number k */
-virtual void PixThetaPhi(int_4 k, double& theta, double& phi) const=0;
-
-/*! 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;
-
-/*! Pixel  (steradians) */
-virtual double PixSolAngle(int_4 k) const =0; 
-
-/*! Setting blockdata to temporary (see ndatablock documentation) */
-virtual void SetTemp(bool temp=false) 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));}
-
-// Overloading of () to access pixel at a sky position .
-inline T& operator()(const SpherePosition& spos) 
-                     { return(PixVal(PixIndex(spos))); }
-inline T const& operator()(const SpherePosition& spos) const 
-                     { return(PixVal(PixIndex(spos))); }
-
-
-// Note : no overloading of (double,double) to access pixel (theta,phi).
-// overloading of (double,double) in SphericalMap
-// overloading of (int,int)     in CartesianMap
-
-/*! Setting pixel values to a constant */
-virtual T SetPixels(T v);
-inline T operator = (T v) { return(SetPixels(v)); }
-
-
-                     
-//++
-DVList& Info()
-//
-// Renvoie une reference sur l''objet DVList Associe
-//--
-  {
-    if (mInfo_ == NULL)  mInfo_ = new DVList;
-    return(*mInfo_);
-  }
-
-const DVList*  ptrInfo() const
-  {
-    return mInfo_;
-  }
-
-protected :
-    SphereCoordSys *cs_;   // Coordinate system used in the map
-    DVList* mInfo_;        // Infos (variables) attachees 
-};
-
-
-template <class T>
-int_4 PixelMap<T>::PixIndex(const SpherePosition& spos)
-{
-UnitVector v = spos.Transform(*cs_);
-return(PixIndexSph(v.Theta(), v.Phi()));
-}
-
-template <class T>
-bool PixelMap<T>::Contains(const SpherePosition& spos) const
-{
-UnitVector v = spos.Transform(*cs_);
-return(ContainsSph(v.Theta(), v.Phi()));
-}
-template <class T>
-T PixelMap<T>::SetPixels(T v)
-{
-int k;
-for(k=0; k<NbPixels(); k++) PixVal(k) = v;
-return(v);
-}
-
-#endif
Index: trunk/SophyaLib/Samba/sambainit.cc
===================================================================
--- trunk/SophyaLib/Samba/sambainit.cc	(revision 757)
+++ trunk/SophyaLib/Samba/sambainit.cc	(revision 758)
@@ -1,35 +1,14 @@
-// Classe d'initialisation du module Outils++ 
-#include <unistd.h>
+// Classe d'initialisation du module Samba 
 
 #include "sambainit.h"
-#include "ppersist.h"
-#include "objfio.h"
-#include "spherethetaphi.h"
-#include "spheregorski.h"
-#include "spherepos.h"
 
 int SambaInitiator::FgInit = 0;
 
 SambaInitiator::SambaInitiator()
-  : PeidaOutilsInitiator()
+  : SkyMapInitiator()
 {
   SambaInitiator::FgInit++;
   if (SambaInitiator::FgInit > 1)  return;
-
-//   Enregistrement des classes PPersist du modules Outils++
-  
-  PPRegister(FIO_SphereCoordSys);
-
-  PPRegister(FIO_SphereThetaPhi<float>);
-  PPRegister(FIO_SphereThetaPhi<double>);
-  PPRegister(FIO_SphereThetaPhi< complex<float> >);
-  PPRegister(FIO_SphereThetaPhi< complex<double> >);
-
-  PPRegister(FIO_SphereGorski<float>);
-  PPRegister(FIO_SphereGorski<double>);
-  PPRegister(FIO_SphereGorski< complex<float> >);
-  PPRegister(FIO_SphereGorski< complex<double> >);
-
-
+  ntinit = new NToolsInitiator;
 }
 
@@ -37,4 +16,5 @@
 {
   SambaInitiator::FgInit--;
+  if (FgInit < 1)  delete ntinit;
 }
 
Index: trunk/SophyaLib/Samba/sambainit.h
===================================================================
--- trunk/SophyaLib/Samba/sambainit.h	(revision 757)
+++ trunk/SophyaLib/Samba/sambainit.h	(revision 758)
@@ -7,9 +7,11 @@
 
 // Classe d''initialisation pour les modules de PEIDA++
-#include "outilsinit.h"
+#include "skymapinit.h"
+#include "ntoolsinit.h"
 
-class SambaInitiator : public PeidaOutilsInitiator {
+class SambaInitiator : public SkyMapInitiator {
 private:
   static int FgInit;
+  NToolsInitiator * ntinit; 
 public:
   			SambaInitiator();
@@ -17,6 +19,6 @@
 };
 
-#undef PeidaInit
-#define PeidaInit() SambaInitiator _sambainitiator_
+#undef SophyaInit
+#define SophyaInit() SambaInitiator sambainitiator
 
 #endif
Index: trunk/SophyaLib/Samba/spheregorski.cc
===================================================================
--- trunk/SophyaLib/Samba/spheregorski.cc	(revision 757)
+++ 	(revision )
@@ -1,1315 +1,0 @@
-#include "machdefs.h"
-#include <math.h>
-#include <complex>
-
-#include "pexceptions.h"
-#include "spheregorski.h"
-#include "strutil.h"
-
-extern "C" 
-{
-#include <stdio.h>
-#include <stdlib.h>
-#include <unistd.h>
-}
-      
-
-//*******************************************************************
-// Class PIXELS_XY 
-// Construction des tableaux necessaires a la traduction des indices RING en 
-// indices NESTED (ou l'inverse)
-//*******************************************************************
-
-PIXELS_XY::PIXELS_XY()
-{
-  pix2x_.ReSize(1024);
-  pix2x_.Reset();
-  pix2y_.ReSize(1024);
-  pix2y_.Reset();
-  x2pix_.ReSize(128);
-  x2pix_.Reset();
-  y2pix_.ReSize(128);
-  y2pix_.Reset();
-  mk_pix2xy(); 
-  mk_xy2pix(); 
-}
-
-PIXELS_XY& PIXELS_XY::instance()
-{
-  static PIXELS_XY single;
-  return (single);
-}
-
-void PIXELS_XY::mk_pix2xy() 
-{
-  /*
-    ==================================================
-    subroutine mk_pix2xy
-    ==================================================
-    c     constructs the array giving x and y in the face from pixel number
-    c     for the nested (quad-cube like) ordering of pixels
-    c
-    c     the bits corresponding to x and y are interleaved in the pixel number
-    c     one breaks up the pixel number by even and odd bits
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int kpix, jpix, IX, IY, IP, ID;
-  
-  for(kpix = 0; kpix < 1024; kpix++) 
-    {
-      jpix = kpix;
-      IX = 0;
-      IY = 0;
-      IP = 1 ;//  ! bit position (in x and y)
-      while( jpix!=0 )
-	{ // ! go through all the bits
-	  ID=jpix%2;//  ! bit value (in kpix), goes in ix
-	  jpix = jpix/2;
-	  IX = ID*IP+IX;
-      
-	  ID=jpix%2;//  ! bit value (in kpix), goes in iy
-	  jpix = jpix/2;
-	  IY = ID*IP+IY;
-     
-	  IP = 2*IP;//        ! next bit (in x and y)
-	}
-      pix2x_(kpix) = IX;//     ! in 0,31
-      pix2y_(kpix) = IY;//     ! in 0,31
-    }
-}
-
-void PIXELS_XY::mk_xy2pix() 
-{
-  /*
-    =================================================
-    subroutine mk_xy2pix
-    =================================================
-    c     sets the array giving the number of the pixel lying in (x,y)
-    c     x and y are in {1,128}
-    c     the pixel number is in {0,128**2-1}
-    c
-    c     if  i-1 = sum_p=0  b_p * 2^p
-    c     then ix = sum_p=0  b_p * 4^p
-    c          iy = 2*ix
-    c     ix + iy in {0, 128**2 -1}
-    =================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int K,IP,I,J,ID;
-  for(I = 1; I <= 128; I++) 
-    {
-      J  = I-1;//            !pixel numbers
-      K  = 0;//
-      IP = 1;//
-      truc : if( J==0 ) 
-	{
-	  x2pix_(I-1) = K;
-	  y2pix_(I-1) = 2*K;
-	}
-      else 
-	{
-	  ID = (int)fmod(J,2);
-	  J  = J/2;
-	  K  = IP*ID+K;
-	  IP = IP*4;
-	  goto truc;
-	}
-    }     
-}
-
-//*******************************************************************
-//++ 
-// Class	SphereGorski
-//
-// include      spheregorski.h strutil.h
-//
-//    Pixelisation Gorski  
-//
-//
-//|    -----------------------------------------------------------------------
-//|     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
-//|    -----------------------------------------------------------------------
-//
-//    the sphere is split in 12 diamond-faces containing nside**2 pixels each
-//
-//    the numbering of the pixels (in the nested scheme) is similar to
-//    quad-cube
-//    In each face the first pixel is in the lowest corner of the diamond
-//
-//    the faces are                    (x,y) coordinate on each face
-//|        .   .   .   .   <--- North Pole
-//|       / \ / \ / \ / \                          ^        ^     
-//|      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      /      
-//|       \ / \ / \ / \ /                        y   \    /  x   
-//|      4 . 5 . 6 . 7 . 4 <--- equator               \  /                     
-//|       / \ / \ / \ / \                              \/        
-//|      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner  
-//|       \ / \ / \ / \ /                                  
-//|        .   .   .   .   <--- South Pole
-//|
-//    phi:0               2Pi                                 
-//
-//    in the ring scheme pixels are numbered along the parallels
-//    the first parallel is the one closest to the north pole and so on
-//    on each parallel, pixels are numbered starting from the one closest
-//    to phi = 0
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-//++
-//
-// Links	Parents
-//
-//    SphericalMap 
-//-- 
-
-/* --Methode-- */
-//++
-// Titre	Constructors
-//--
-//++
-
-template<class T>
-SphereGorski<T>::SphereGorski()
-
-//--
-{
-  InitNul();
-  pixels_.Reset();
-}
-
-//++
-template<class T>
-SphereGorski<T>::SphereGorski(int_4 m)
-
-//    m is the "nside" of the Gorski algorithm
-//
-//    The total number of pixels will be Npix =  12*nside**2
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-{
-
-  if(m <= 0 || m > 8192) 
-    {
-      cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
-      throw RangeCheckError("SphereGorski<T>::SphereGorski() - Out of bound nside (< 8192)!");
-    }
-  // verifier que m est une puissance de deux 
-  int x= m;
-  while(x%2 == 0) x/=2;
-  if(x != 1) 
-    {  
-      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
-      throw ParmError("SphereGorski<T>::SphereGorski() - nside != 2^n !");
-    }
-  InitNul();
-  Pixelize(m);
-  SetThetaSlices();
-}
-//++
-template<class T>
-SphereGorski<T>::SphereGorski(const SphereGorski<T>& s, bool share)
-  : pixels_(s.pixels_, share)
-//    copy constructor
-//--
-{
-  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
-
-  nSide_= s.nSide_;
-  nPix_ = s.nPix_;
-  omeg_ = s.omeg_;
-}
-
-//++
-// Titre	Destructor
-//--
-//++
-template<class T>
-SphereGorski<T>::~SphereGorski()
-
-//--
-{
-}
-
-//++
-// Titre	Public Methods
-//--
-
-//++
-template<class T>
-void SphereGorski<T>::Resize(int_4 m)
-
-//    m is the "nside" of the Gorski algorithm
-//
-//    The total number of pixels will be Npix =  12*nside**2
-//
-//    nside MUST be a power of 2 (<= 8192)
-//--
-{
-  if (m<=0 || m> 8192) {
-    cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
-    exit(1);
-  }
-  // verifier que m est une puissance de deux 
-  int x= m;
-  while (x%2==0) x/=2;
-  if(x != 1) 
-    {  
-      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
-      exit(1);
-    }
-  InitNul();
-  Pixelize(m);
-  SetThetaSlices();
-}
-
-template<class T>
-void  SphereGorski<T>::Pixelize( int_4 m) 
-
-//    prépare la pixelisation Gorski (m a la même signification 
-//    que pour le constructeur)
-//
-//    
-{
-  // On memorise les arguments d'appel
-  nSide_= m;  
-
-  // Nombre total de pixels sur la sphere entiere
-  nPix_= 12*nSide_*nSide_;
-
-  // pour le moment les tableaux qui suivent seront ranges dans l'ordre 
-  // de l'indexation GORSKY "RING"
-  // on pourra ulterieurement changer de strategie et tirer profit
-  // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra 
-  // de pourquoi c'est faire
-
-  // Creation et initialisation du vecteur des contenus des pixels 
-  pixels_.ReSize(nPix_);
-  pixels_.Reset();
-
-  // solid angle per pixel   
-  omeg_= 4.0*Pi/nPix_;
-}
-
-template<class T>
-void SphereGorski<T>::InitNul()
-//
-//    initialise à zéro les variables de classe 
-{
-  nSide_= 0;
-  nPix_ = 0;
-  omeg_ = 0.;
-//  pixels_.Reset();  -  Il ne faut pas mettre les pixels a zero si share !
-}
-
-/* --Methode-- */
-//++
-template<class T>
-int_4 SphereGorski<T>::NbPixels() const
-
-//    Retourne le nombre de pixels du découpage
-//--
-{  
-  return(nPix_);
-}
-
-//++
-template<class T>
-uint_4 SphereGorski<T>::NbThetaSlices() const 
-
-//    Return number of slices in theta direction on the  sphere
-//--
-{
-  uint_4 nbslices = uint_4(4*nSide_-1);
-  if (nSide_<=0) 
-    {
-      nbslices = 0;
-      throw PException(" sphere not pixelized, NbSlice=0 ");
-    }
-  return nbslices;
-}
-
-//++
-template<class T>
-void SphereGorski<T>::GetThetaSlice(int_4 index, r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 
-
-//    For a theta-slice with index 'index', return : 
-//
-//    the corresponding "theta" 
-//
-//    a vector containing the phi's of the pixels of the slice
-//
-//    a vector containing the corresponding values of pixels 
-//
-//--
-{
-
-  if (index<0 || index >= NbThetaSlices()) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::GetThetaSlice : Pixel index out of range" <<endl;
-      THROW(rangeCheckErr);
-    }
-
-
-  int_4 iring= sliceBeginIndex_(index);
-  int_4 lring  =  sliceLenght_(index);
-
-
-  phi.ReSize(lring);
-  value.ReSize(lring);
-  double TH= 0.;
-  double FI= 0.;
-  for(int_4 kk = 0; kk < lring;kk++) 
-    {
-      PixThetaPhi(kk+iring,TH,FI);
-      phi(kk)= FI;
-      value(kk)= PixVal(kk+iring);
-    }
-  theta= TH;
-}
-//++
-//++
-
-template<class T>
-void SphereGorski<T>::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const 
-
-//    For a theta-slice with index 'sliceIndex', return : 
-//
-//    the corresponding "theta" 
-//    the corresponding "phi" for first pixel of the slice 
-//
-//    a vector containing the indices of the pixels of the slice
-//    (equally distributed in phi)
-//
-//    a vector containing the corresponding values of pixels 
-//
-//--
-{
-
-  if (sliceIndex<0 || sliceIndex >= NbThetaSlices()) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::GetThetaSlice : Pixel index out of range" <<endl;
-      THROW(rangeCheckErr);
-    }
-
-  int_4 iring= sliceBeginIndex_(sliceIndex);
-  int_4 lring  =  sliceLenght_(sliceIndex);
-  pixelIndices.ReSize(lring);
-  value.ReSize(lring);
-  for(int_4 kk = 0; kk < lring;kk++) 
-    {
-      pixelIndices(kk)= kk+iring;
-      value(kk)= PixVal(kk+iring);
-    }
-  PixThetaPhi(iring, theta, phi0);
-}
-//++
-template<class T>
-void SphereGorski<T>::SetThetaSlices()  
-
-//--
-{
-  sliceBeginIndex_.ReSize(4*nSide_-1);
-  sliceLenght_.ReSize(4*nSide_-1);
-  for (int sliceIndex=0; sliceIndex<  nSide_-1; sliceIndex++)
-    {
-      sliceBeginIndex_(sliceIndex)  = 2*sliceIndex*(sliceIndex+1);
-      sliceLenght_(sliceIndex) = 4*(sliceIndex+1);
-    }
-  for (int sliceIndex= nSide_-1; sliceIndex<  3*nSide_; sliceIndex++)
-    {
-      sliceBeginIndex_(sliceIndex)  = 2*nSide_*(2*sliceIndex-nSide_+1);
-      sliceLenght_(sliceIndex) = 4*nSide_;
-    }
-  for (int sliceIndex= 3*nSide_; sliceIndex< 4*nSide_-1; sliceIndex++)
-    {
-      int_4 nc= 4*nSide_-1-sliceIndex;
-      sliceBeginIndex_(sliceIndex)  = nPix_-2*nc*(nc+1); 
-      sliceLenght_(sliceIndex) = 4*nc;
-    }
-}
-
-/* --Methode-- */
-//++
-template<class T>
-T& SphereGorski<T>::PixVal(int_4 k)
-
-//    Return value of  pixel with "RING" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::PIxVal : exception  a mettre en place, k= " << k <<endl;
-      THROW(rangeCheckErr);
-    }
-  return pixels_(k);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-T const& SphereGorski<T>::PixVal(int_4 k) const
-
-//    Return value of  pixel with "RING" index k 
-//--
-{
-  if((k < 0) || (k >= nPix_)) 
-    {
-      //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
-      cout << " SphereGorski::PIxVal : exception  a mettre en place, k= " << k <<endl;
-      THROW(rangeCheckErr);
-    }
-  return *(pixels_.Data()+k);
-}
-
-//++
-template<class T>
-T& SphereGorski<T>::PixValNest(int_4 k)
-
-//    Return value of  pixel with "NESTED" index 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 pixels_(nest2ring(nSide_,k));
-}
-//++
-
-template<class T>
-T const& SphereGorski<T>::PixValNest(int_4 k) const
-
-//    Return value of  pixel with "NESTED" index 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 *(pixels_.Data()+pix);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-bool SphereGorski<T>::ContainsSph(double /*theta*/, double /*phi*/) const
-//--
-{
-return(true);
-}
-
-/* --Methode-- */
-//++
-template<class T>
-int_4 SphereGorski<T>::PixIndexSph(double theta,double phi) const
-
-//    Return "RING" index of the pixel corresponding to 
-//    direction (theta, phi).
-//--
-{
-  return ang2pix_ring(nSide_,theta,phi);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::PixIndexSphNest(double theta,double  phi) const
-
-//    Return "NESTED" index of the pixel corresponding to 
-//    direction (theta, phi).
-//--
-{
-  return ang2pix_nest(nSide_,theta,phi);
-}
-
-
-/* --Methode-- */
-//++
-template<class T>
-void SphereGorski<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
-
-//   Return (theta,phi) coordinates of middle of  pixel with "RING" index k
-//--
-{
-  pix2ang_ring(nSide_,k,theta,phi);
-}
-
-template <class T>
-T SphereGorski<T>::SetPixels(T v)
-{
-pixels_.Reset(v);
-return(v);
-}
-
-//++
-template<class T>
-double SphereGorski<T>::PixSolAngle(int_4 /*dummy*/) const
-//    Pixel Solid angle  (steradians)
-//    All the pixels have the same solid angle. The dummy argument is
-//    for compatibility with eventual pixelizations which would not 
-//    fulfil this requirement.
-//--
-{
-  return omeg_;
-}
-
-//++
-template<class T>
-void SphereGorski<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
-
-//   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
-//--
-{   
-  pix2ang_nest(nSide_,k,theta,phi);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::NestToRing(int_4 k) const
-
-//    translation from NESTED index  into RING index 
-//
-//--
-{
-  return  nest2ring(nSide_,k);
-}
-
-//++
-template<class T>
-int_4 SphereGorski<T>::RingToNest(int_4 k) const
-//
-//    translation from  RING index  into NESTED index 
-//
-//--
-{
-  return  ring2nest(nSide_,k);
-}
-
-
-template<class T>
-int_4 SphereGorski<T>::nest2ring(int_4 nside, int_4 ipnest) const 
-{
-  /*
-    ====================================================
-    subroutine nest2ring(nside, ipnest, ipring)
-    ====================================================
-    c     conversion from NESTED to RING pixel number
-    ====================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  int npix, npface, face_num, ncap, n_before;
-  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
-  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
-  int ns_max=8192;
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
-  
-  if(  nside<1 ||  nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 *  nside* nside;
-  if( ipnest<0 || ipnest>npix-1 ) {
-    cout << "ipnest out of range" << endl;
-    exit(0);
-  }
-
-  ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
-  nl4   = 4* nside;
-  
-  //c     finds the face, and the number in the face
-  npface =  nside* nside;
-  //cccccc      ip = ipnest - 1         ! in {0,npix-1}
-  
-  face_num = ipnest/npface;//  ! face number in {0,11}
-  ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
-  //c     finds the x,y on the face (starting from the lowest corner)
-  //c     from the pixel number
-  ip_low=ipf%1024;                //   ! content of the last 10 bits
-  ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
-  ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
-  ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
-  
-  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
-  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
-  
-  //c     transforms this in (horizontal, vertical) coordinates
-  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
-  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
-  
-  //c     computes the z coordinate on the sphere
-  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
-  jr =  jrll[face_num]*nside - jrt - 1;
-  nr = nside;//                  ! equatorial region (the most frequent)
-  n_before = ncap + nl4 * (jr - nside);
-  kshift=(jr - nside)%2;
-  if( jr<nside ) {//then     ! north pole region
-    nr = jr;
-    n_before = 2 * nr * (nr - 1);
-    kshift = 0;
-  }
-  else if( jr>3*nside ) {//then ! south pole region
-    nr = nl4 - jr;
-    n_before = npix - 2 * (nr + 1) * nr;
-    kshift = 0;
-  }
-  
-  //c     computes the phi coordinate on the sphere, in [0,2Pi]
-  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
-  
-  if( jp>nl4 ) jp = jp - nl4;
-  if( jp<1 )   jp = jp + nl4;
-  
-  int aux=n_before + jp - 1;
-  return (n_before + jp - 1);// ! in {0, npix-1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ring2nest(int_4 nside, int_4 ipring) const 
-{
-  /*
-    ==================================================
-    subroutine ring2nest(nside, ipring, ipnest)
-    ==================================================
-    c     conversion from RING to NESTED pixel number
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  double fihip, hip;
-  int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
-  int     kshift, face_num, nr;
-  int irn, ire, irm, irs, irt, ifm , ifp;
-  int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
-  int ns_max(8192);
-  
-  // coordinate of the lowest corner of each face
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 * nside*nside;
-  if( ipring<0 || ipring>npix-1 ) {
-    cout << "ipring out of range" << endl;
-    exit(0);
-  }
-  
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  npix = 12*nside*nside;//      ! total number of points
-  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
-  ipring1 = ipring + 1;
-  
-  //c     finds the ring number, the position of the ring and the face number
-  if( ipring1<=ncap ) {//then
-    
-    hip   = ipring1/2.;
-    fihip = floor ( hip );
-    irn   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
-    iphi  = ipring1 - 2*irn*(irn - 1);
-    
-    kshift = 0;
-    nr = irn   ;//               ! 1/4 of the number of points on the current ring
-    face_num = (iphi-1) / irn;// ! in {0,3}
-  }
-  else if( ipring1<=nl2*(5*nside+1) ) {//then
-    
-    ip    = ipring1 - ncap - 1;
-    irn   = (int)floor( ip / nl4 ) + nside;//               ! counted from North pole
-    iphi  = (int)fmod(ip,nl4) + 1;
-    
-    kshift  = (int)fmod(irn+nside,2);//  ! 1 if irn+nside is odd, 0 otherwise
-    nr = nside;
-    ire =  irn - nside + 1;// ! in {1, 2*nside +1}
-    irm =  nl2 + 2 - ire;
-    ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
-    ifp = (iphi - irm/2 + nside -1) / nside;
-    if( ifp==ifm ) {//then          ! faces 4 to 7
-      face_num = (int)fmod(ifp,4) + 4;
-    }
-    else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
-      face_num = ifp;
-    }
-    else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
-      face_num = ifp + 7;
-    }
-  }
-  else {
-    
-    ip    = npix - ipring1 + 1;
-    hip   = ip/2.;
-    fihip = floor ( hip );
-    irs   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//  ! counted from South pole
-    iphi  = 4*irs + 1 - (ip - 2*irs*(irs-1));
-    
-    kshift = 0;
-    nr = irs;
-    irn   = nl4 - irs;
-    face_num = (iphi-1) / irs + 8;// ! in {8,11}
-  }
-  
-  //c     finds the (x,y) on the face
-  irt =   irn  - jrll[face_num]*nside + 1;//       ! in {-nside+1,0}
-  ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
-
-
-  if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
-  
-  ix =  (ipt - irt ) / 2;
-  iy = -(ipt + irt ) / 2;
-  
-  ix_low = (int)fmod(ix,128);
-  ix_hi  = ix/128;
-  iy_low = (int)fmod(iy,128);
-  iy_hi  = iy/128;
-  ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
-
-  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ang2pix_ring(int_4 nside, double theta, double phi) const 
-{
-  /*
-    ==================================================
-    c     gives the pixel number ipix (RING) 
-    c     corresponding to angles theta and phi
-    c==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  int nl2, nl4, ncap, npix, jp, jm, ipix1;
-  double  z, za, tt, tp, tmp;
-  int ir, ip, kshift;
-
-  double piover2(Pi/2.);
-  double twopi(2.*Pi);
-  double z0(2./3.);
-  int ns_max(8192);
-
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-
-  if( theta<0. || theta>Pi) {
-    cout << "theta out of range" << endl;
-    exit(0);
-  }
-
-  z = cos(theta);
-  za = fabs(z);
-  if( phi >= twopi)  phi = phi - twopi;
-  if (phi < 0.)     phi = phi + twopi;
-  tt = phi / piover2;//  ! in [0,4)
-
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  ncap  = nl2*(nside-1);// ! number of pixels in the north polar cap
-  npix  = 12*nside*nside;
-
-  if( za <= z0 ) {
-    
-    jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of  ascending edge line 
-    jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
-
-    ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
-    kshift = 0;
-    if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
-    
-    ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
-    if( ip>nl4 ) ip = ip - nl4;
-    
-    ipix1 = ncap + nl4*(ir-1) + ip ;
-  }
-  else {
-
-    tp = tt - floor(tt);//      !MOD(tt,1.d0)
-    tmp = sqrt( 3.*(1. - za) );
-
-    jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
-    jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
-
-    ir = jp + jm + 1;//        ! ring number counted from the closest pole
-    ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
-    if( ip>4*ir ) ip = ip - 4*ir;
-      
-    ipix1 = 2*ir*(ir-1) + ip;
-      if( z<=0. ) {
-	ipix1 = npix - 2*ir*(ir+1) + ip;
-      }
-  }
-    return (ipix1 - 1);// ! in {0, npix-1}
-}
-
-template<class T>
-int_4 SphereGorski<T>::ang2pix_nest(int_4 nside, double theta, double phi) const 
-{
-  /*
-    ==================================================
-    subroutine ang2pix_nest(nside, theta, phi, ipix)
-    ==================================================
-    c     gives the pixel number ipix (NESTED) 
-    c     corresponding to angles theta and phi
-    c
-    c     the computation is made to the highest resolution available (nside=8192)
-    c     and then degraded to that required (by integer division)
-    c     this doesn't cost more, and it makes sure 
-    c     that the treatement of round-off will be consistent 
-    c     for every resolution
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  double    z, za, z0, tt, tp, tmp;
-  int face_num,jp,jm;
-  int ifp, ifm;
-  int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
-  double piover2(Pi/2.), twopi(2.*Pi);
-  int ns_max(8192);
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  if( theta<0 || theta>Pi ) {
-    cout << "theta out of range" << endl;
-    exit(0);
-  }
-  z  = cos(theta);
-  za = fabs(z);
-  z0 = 2./3.;
-  if( phi>=twopi ) phi = phi - twopi;
-  if( phi<0. )    phi = phi + twopi;
-  tt = phi / piover2;// ! in [0,4[
-  if( za<=z0 ) { // then ! equatorial region
-    
-    //(the index of edge lines increase when the longitude=phi goes up)
-    jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
-    jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
-    
-    //c        finds the face
-    ifp = jp / ns_max;//  ! in {0,4}
-    ifm = jm / ns_max;
-    if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then  ! faces 4 to 7
-    else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
-    else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
-    
-    ix = (int)fmod(jm, ns_max);
-    iy = ns_max - (int)fmod(jp, ns_max) - 1;
-  }
-  else { //! polar region, za > 2/3
-    
-    ntt = (int)floor(tt);
-    if( ntt>=4 ) ntt = 3;
-    tp = tt - ntt;
-    tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
-    
-    //(the index of edge lines increase when distance from the closest pole goes up)
-    jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
-    jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
-    jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
-    jm = (int)min(ns_max-1, jm);
-    
-    // finds the face and pixel's (x,y)
-    if( z>=0 ) {
-      face_num = ntt;//  ! in {0,3}
-      ix = ns_max - jm - 1;
-      iy = ns_max - jp - 1;
-    }
-    else {
-      face_num = ntt + 8;// ! in {8,11}
-      ix =  jp;
-      iy =  jm;
-    }
-  }
-  
-  ix_low = (int)fmod(ix,128);
-  ix_hi  =     ix/128;
-  iy_low = (int)fmod(iy,128);
-  iy_hi  =     iy/128;
-  ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
-  //  ipf = ipf / pow(ns_max/nside,2.);//  ! in {0, nside**2 - 1}
-  //  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
-  // $CHECK$  Reza 25/10/99 , pow remplace par *
-  ipf = ipf / ((ns_max/nside)*(ns_max/nside));
-  return (ipf + face_num*nside*nside);
-}
-
-template<class T>
-void SphereGorski<T>::pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const {
-  /*
-    ===================================================
-    c     gives theta and phi corresponding to pixel ipix (RING) 
-    c     for a parameter nside
-    ===================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-  
-  int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
-  double  fact1, fact2, fodd, hip, fihip;
-  
-  int ns_max(8192);
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12*nside*nside;      // ! total number of points
-  if( ipix<0 || ipix>npix-1 ) {
-    cout << "ipix out of range" << endl;
-    exit(0);
-  }
-  
-  ipix1 = ipix + 1; // in {1, npix}
-  nl2 = 2*nside;
-  nl4 = 4*nside;
-  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
-  fact1 = 1.5*nside;
-  fact2 = 3.0*nside*nside;
-  
-  if( ipix1 <= ncap ) {  //! North Polar cap -------------
-    
-    hip   = ipix1/2.;
-    fihip = floor(hip);
-    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
-    iphi  = ipix1 - 2*iring*(iring - 1);
-    
-    theta = acos( 1. - iring*iring / fact2 );
-    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
-    //    cout << theta << " " << phi << endl;
-  }
-  else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
-    
-    ip    = ipix1 - ncap - 1;
-    iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
-    iphi  = ip%nl4 + 1;
-    
-    fodd  = 0.5 * (1 + (iring+nside)%2 );//  ! 1 if iring+nside is odd, 1/2 otherwise
-    theta = acos( (nl2 - iring) / fact1 );
-    phi   = ((double)iphi - fodd) * Pi /(2.*nside);
-  }
-  else {//! South Polar cap -----------------------------------
-    
-    ip    = npix - ipix1 + 1;
-    hip   = ip/2.;
-    fihip = floor(hip);
-    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//     ! counted from South pole
-    iphi  = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
-    
-    theta = acos( -1. + iring*iring / fact2 );
-    phi   = ((double)iphi - 0.5) * Pi/(2.*iring);
-    //    cout << theta << " " << phi << endl;
-  }
-}
-
-template<class T>
-void SphereGorski<T>::pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const {
-  /*
-    ==================================================
-    subroutine pix2ang_nest(nside, ipix, theta, phi)
-    ==================================================
-    c     gives theta and phi corresponding to pixel ipix (NESTED) 
-    c     for a parameter nside
-    ==================================================
-  */
-  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
-  //  (16/12/98)
-
-  const PIXELS_XY& PXY= PIXELS_XY::instance();
-    
-  int npix, npface, face_num;
-  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
-  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
-  double z, fn, fact1, fact2;
-  double piover2(Pi/2.);
-  int ns_max(8192);
-          
-  // ! coordinate of the lowest corner of each face
-  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
-  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
-  
-  if( nside<1 || nside>ns_max ) {
-    cout << "nside out of range" << endl;
-    exit(0);
-  }
-  npix = 12 * nside*nside;
-  if( ipix<0 || ipix>npix-1 ) {
-    cout << "ipix out of range" << endl;
-    exit(0);
-  }
-  
-  fn = 1.*nside;
-  fact1 = 1./(3.*fn*fn);
-  fact2 = 2./(3.*fn);
-  nl4   = 4*nside;
-
-  //c     finds the face, and the number in the face
-  npface = nside*nside;
-  
-  face_num = ipix/npface;//  ! face number in {0,11}
-  ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
-  
-  //c     finds the x,y on the face (starting from the lowest corner)
-  //c     from the pixel number
-  ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
-  ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
-  ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
-  ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
-  
-  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
-  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
-  
-  //c     transforms this in (horizontal, vertical) coordinates
-  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
-  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
-  
-  //c     computes the z coordinate on the sphere
-  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
-  jr =  jrll[face_num]*nside - jrt - 1;
-  nr = nside;//                  ! equatorial region (the most frequent)
-  z  = (2*nside-jr)*fact2;
-  kshift = (int)fmod(jr - nside, 2);
-  if( jr<nside ) { //then     ! north pole region
-    nr = jr;
-    z = 1. - nr*nr*fact1;
-    kshift = 0;
-  }
-  else {
-    if( jr>3*nside ) {// then ! south pole region
-      nr = nl4 - jr;
-      z = - 1. + nr*nr*fact1;
-      kshift = 0;
-    }
-  }
-  theta = acos(z);
-  
-  //c     computes the phi coordinate on the sphere, in [0,2Pi]
-  //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
-  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
-  if( jp>nl4 ) jp = jp - nl4;
-  if( jp<1 )   jp = jp + nl4;
-  phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
-}
-
-
-
-template <class T>
-void 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 << " content of pixels : ";
-  for(int i=0; i < nPix_; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os <<  pixels_(i) <<", ";
-    }
-  os << endl;
-
-  os << endl;
-  //const PIXELS_XY& PXY= PIXELS_XY::instance();
-
-  //os << endl;  os << " contenu des tableaux conversions "<<endl;
-  //for(int i=0; i < 5; i++)
-  //  {
-  //   os<<PXY.pix2x_(i)<<", "<<PXY.pix2y_(i)<<", "<<PXY.x2pix_(i)<<", "<<PXY.y2pix_(i)<<endl;
-  // }
-  os << endl;
-
-}
-
-//*******************************************************************
-// 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>;
-  dobj->pixels_.SetTemp(true);
-  ownobj= true;
-  Read(filename);
-}
-
-template <class T>
-FIO_SphereGorski<T>::FIO_SphereGorski(const SphereGorski<T>& obj)
-{
-  dobj= new SphereGorski<T>(obj, true);
-  dobj->pixels_.SetTemp(true);
-  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)
-{
-
-  if(dobj == NULL) 
-    {
-      dobj= new SphereGorski<T>;
-      dobj->pixels_.SetTemp(true);
-      ownobj= true;      
-    }
-
-// On lit les 3 premiers uint_8
-uint_8 itab[3];
-is.Get(itab, 3);
-// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Read(is);
-  // 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);
-  int_4 nPix;
-  is.GetI4(nPix);
-  double Omega;
-  is.GetR8(Omega);
-  dobj->setParameters(nSide,nPix, Omega);
-
-// On lit les DataBlocks;
-  FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-  fio_nd.Read(is);
-  FIO_NDataBlock<int> fio_ndi1(&dobj->sliceBeginIndex_);
-  fio_ndi1.Read(is);
-  FIO_NDataBlock<int> fio_ndi2(&dobj->sliceLenght_);
-  fio_ndi2.Read(is);
-}
-
-template <class T>
-void FIO_SphereGorski<T>::WriteSelf(POutPersist& os) const
-{
-  if(dobj == NULL) 
-    {
-      cout << " WriteSelf:: dobj= null " << endl;
-      return;
-    }
-
-//  On ecrit 3 uint_8 
-//  0 : Numero de version,  1 : Size index,  2  reserve a l
-uint_8 itab[3];
-itab[0] = 1;
-itab[1] = dobj->SizeIndex();
-itab[2] = 0;
-os.Put(itab, 3);
-// Let's write the SphereCoordSys object
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Write(os);
-
-  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());
-// On ecrit les dataBlocks 
-  FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-  fio_nd.Write(os);
-  FIO_NDataBlock<int> fio_ndi1(&dobj->sliceBeginIndex_);
-  fio_ndi1.Write(os);
-  FIO_NDataBlock<int> fio_ndi2(&dobj->sliceLenght_);
-  fio_ndi2.Write(os);
-}
-
-#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
-
Index: trunk/SophyaLib/Samba/spheregorski.h
===================================================================
--- trunk/SophyaLib/Samba/spheregorski.h	(revision 757)
+++ 	(revision )
@@ -1,275 +1,0 @@
-#ifndef SPHEREGORSKI_SEEN
-#define SPHEREGORSKI_SEEN
-
-#include "sphericalmap.h"
-#include "tvector.h"
-#include "ndatablock.h"
-
-#include "anydataobj.h"
-#include "ppersist.h"
-
-
-// ***************** CLASSE SphereGorski *****************************
-
-//!  class SphereGorski
-/*!
-   Pixelisation Gorski  
-
-
-    -----------------------------------------------------------------------
-     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
-    -----------------------------------------------------------------------
-
-    the sphere is split in 12 diamond-faces containing nside**2 pixels each
-
-    the numbering of the pixels (in the nested scheme) is similar to
-    quad-cube
-    In each face the first pixel is in the lowest corner of the diamond
-
-    the faces are                    (x,y) coordinate on each face
-\verbatim
-        .   .   .   .   <--- North Pole
-       / \ / \ / \ / \                          ^        ^     
-      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      / 
-       \ / \ / \ / \ /                        y   \    /  x  
-      4 . 5 . 6 . 7 . 4 <--- equator               \  /      
-       / \ / \ / \ / \                              \/      
-      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner 
-       \ / \ / \ / \ /      
-        .   .   .   .   <--- South Pole
-\endverbatim
-    phi:0               2Pi                                 
-
-   in the ring scheme pixels are numbered along the parallels
-   the first parallel is the one closest to the north pole and so on  
-   on each parallel, pixels are numbered starting from the one closest
-    to phi = 0
-
-    nside MUST be a power of 2 (<= 8192)
-
-*/
-
-template<class T>
-class FIO_SphereGorski;
-
-template<class T>
-class SphereGorski : public SphericalMap<T>
-{
-
-  friend class FIO_SphereGorski<T>;
-
-public :
-
-SphereGorski();
-/*!    
-  m is the "nside" of the Gorski algorithm
-
-  The total number of pixels will be Npix =  12*nside**2
-
-  nside MUST be a power of 2 (<= 8192)
-*/
-SphereGorski(int_4 m);
-SphereGorski(const SphereGorski<T>& s, bool share=false);
-//!	Destructor
-virtual ~SphereGorski();
-
-/*! Setting blockdata to temporary (see ndatablock documentation) */
-inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};
-// ------------------ Definition of PixelMap abstract methods
-
-/* Nombre de pixels du decoupage */
-/*!    Return number of  pixels of the  splitting */
-virtual int_4 NbPixels() const;
-
-/* Valeur du contenu du pixel d'indice "RING" k  */
-/*!    Return value of  pixel with "RING" index k */
-virtual T& PixVal(int_4 k);
-virtual T const& PixVal(int_4 k) const;
-
-/* Nombre de tranches en theta */ 
-/*!    Return number of slices in theta direction on the  sphere */
-uint_4 NbThetaSlices() const;
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-    a vector containing the phi's of the pixels of the slice
-
-    a vector containing the corresponding values of pixels 
-*/
-void GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const;
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-   the corresponding "phi" for first pixel of the slice 
-
-    a vector containing indices of the pixels of the slice
-
-   (equally distributed in phi)
-
-    a vector containing the corresponding values of pixels 
-*/
- void GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const ;
-
-/* Return true if teta,phi in map  */
-virtual bool ContainsSph(double theta, double phi) const;
-/* Indice "RING" du pixel vers lequel pointe une direction definie par 
-ses  coordonnees spheriques */                                    
-/*!  Return "RING" index of the pixel corresponding to direction (theta, phi).
- */
-virtual int_4 PixIndexSph(double theta,double phi) const;
-
-/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
-virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;
-
-/*! Set all pixels to value v */
-virtual T SetPixels(T v);
-
-/* Pixel Solid angle  (steradians) */
-/*!    Pixel Solid angle  (steradians)
-
-    All the pixels have the same solid angle. The dummy argument is
-    for compatibility with eventual pixelizations which would not 
-   fulfil this requirement.
-*/
-virtual double PixSolAngle(int_4 dummy=0) const;
-
-/*  Acces to the DataBlock  */
-inline       NDataBlock<T>& DataBlock()       {return pixels_;}
-inline const NDataBlock<T>& DataBlock() const {return pixels_;}
-
-// --------------- Specific methods 
-
-/*!    
-  m is the "nside" of the Gorski algorithm
-
-  The total number of pixels will be Npix =  12*nside**2
-
-  nside MUST be a power of 2 (<= 8192)
-*/
-virtual void Resize(int_4 m);
-
-// pour l'instant le tableau est ordonne selon RING, uniquement
-inline virtual char* TypeOfMap() const {return "RING";};
-
-
-/* Valeur du contenu du pixel d'indice "NEST" k */
-/*!    Return value of  pixel with "NESTED" index k */
-virtual T& PixValNest(int_4 k);
-/*!    Return value of  pixel with "NESTED" index k */
-virtual T const& PixValNest(int_4 k) const;
-
-/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
-ses  coordonnees spheriques */                                    
-/*! Return "NESTED" index of the pixel corresponding to direction (theta, phi).
- */
-virtual int_4 PixIndexSphNest(double theta,double phi) const;
-
-/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
-/*!   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
- */
-virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;
-
-/* algorithme de pixelisation */
-void Pixelize(int_4); 
-
-/* convertit index nested en ring  */
-/*!    translation from NESTED index  into RING index */
-int_4 NestToRing(int_4) const;
-
-/* convertit index ring en nested" */
-/*!    translation from  RING index  into NESTED index */
-int_4 RingToNest(int_4) const;
-
-
-/* retourne la valeur du parametre Gorski */
-inline virtual int_4 SizeIndex() const {return(nSide_);}
-
-/* retourne les pointeurs /remplit les tableaux */
-
-
-
-
-/* impression */ 
-void print(ostream& os) const;
-         
-private :
-
-// ------------- méthodes internes ----------------------
-void InitNul();
- void SetThetaSlices(); 
-
-int_4  nest2ring(int_4 nside,int_4 ipnest) const;
-int_4  ring2nest(int_4 nside,int_4 ipring) const;
-
-int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
-int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
-void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
-void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
-inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 
-  {
-    nSide_= nside;
-    nPix_= nbpixels;
-    omeg_= solangle;
-  }
-
-// ------------- variables internes -----------------------
-
-int_4 nSide_;
-int_4 nPix_;
-double omeg_;
-
-NDataBlock<T> pixels_;
-NDataBlock<int> sliceBeginIndex_;
- NDataBlock<int> sliceLenght_;
-
-};
-
-//
-// ------------- 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_4> pix2x_;
-NDataBlock<int_4> pix2y_;
-NDataBlock<int_4> x2pix_;
-NDataBlock<int_4> y2pix_;
-
-private :
-
-PIXELS_XY();
-void mk_pix2xy();
-void mk_xy2pix();
-};
-#endif
Index: trunk/SophyaLib/Samba/spherepos.cc
===================================================================
--- trunk/SophyaLib/Samba/spherepos.cc	(revision 757)
+++ 	(revision )
@@ -1,131 +1,0 @@
-// 04/01/00 : implantation de la persistance par classe deleguee - Guy Le Meur
-
-#include "spherepos.h"
-
-static char *head_spherepos_cc_ = "$Header: /Users/garnier/temp/CVSSophya/SophyaLib/Samba/spherepos.cc,v 1.4 2000-01-04 14:27:15 ansari Exp $";
-
-
-//................. SphereCoordSys class .................
-
-SphereCoordSys::SphereCoordSys(){
-  id_ = SphereCoordSys_NEUTRAL;
-  description_ = "NEUTRAL SphereCoordSystem";
-}
-SphereCoordSys::SphereCoordSys(const SphereCoordSys& a )
-{
-  id_ = a.id_;
-  description_ = a.description_;
-}
-
-SphereCoordSys::SphereCoordSys(int id, const string& description){
-  id_ = id;
-  description_ = description;
-}
-
-SphereCoordSys::~SphereCoordSys()
-{
-}
-SphereCoordSys& SphereCoordSys::operator = (const SphereCoordSys& a)
-{
-  id_ = a.id_;
-  description_ = a.description_;
-  return *this;
-}
-
-
-////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------------------
-//   Les objets delegues pour la gestion de persistance 
-// -------------------------------------------------------------------------
-
-FIO_SphereCoordSys::FIO_SphereCoordSys()
-  {
-    dobj=new SphereCoordSys;
-    ownobj=true;
-  }
-FIO_SphereCoordSys::FIO_SphereCoordSys(string const & filename)
-  {
-    dobj=new SphereCoordSys;
-    ownobj=true; 
-    Read(filename);
-  }
-FIO_SphereCoordSys::FIO_SphereCoordSys(const SphereCoordSys & obj)
-  {
-    dobj = new SphereCoordSys(obj);
-    ownobj=true; 
-  }
-FIO_SphereCoordSys::FIO_SphereCoordSys(SphereCoordSys * obj)
-  {
-    dobj = obj;
-    ownobj=false; 
-  }
-FIO_SphereCoordSys::~FIO_SphereCoordSys()
-  {
-    if (ownobj && dobj) delete dobj;
-  }
-AnyDataObj* FIO_SphereCoordSys::DataObj()
-  {
-    return(dobj);
-  }
-
-
-void  FIO_SphereCoordSys::ReadSelf(PInPersist& is)
-  {
-uint_8 itab[3];
- string description;
-is.Get(itab, 3);
-is.GetStr(description);
-if (dobj == NULL) dobj = new SphereCoordSys(itab[1], description);
- else *dobj= SphereCoordSys(itab[1], description);
-  }           
-void       FIO_SphereCoordSys::WriteSelf(POutPersist& os) const
-{
-if (dobj == NULL)   return;  // Attention - $CHECK$ Guy 04/01/00
-//  On ecrit 3 uint_8 
-//  0 : Numero de version,  1 : Id,  2  reserve a l
-uint_8 itab[3];
-itab[0] = 1;
-itab[1] = dobj->Id();
-itab[2] = 0;
-os.Put(itab, 3);
-os.PutStr(dobj->description());
-} 
-
-
-
-//................. SpherePosition class .................
-
-SpherePosition::SpherePosition() 
-        : UnitVector()
-{
-  cs_ = new SphereCoordSys;
-}
-
-SpherePosition::SpherePosition(double theta, double phi, SphereCoordSys* cs)
-        : UnitVector(theta, phi)
-{
-if (cs == NULL) cs_ = new SphereCoordSys;
-else cs_ = cs;
-}
-
-SpherePosition::SpherePosition(double x, double y, double z, SphereCoordSys* cs)
-        : UnitVector(x, y, z)
-{
-if (cs == NULL) cs_ = new SphereCoordSys;
-else cs_ = cs;
-}
-
-SpherePosition::~SpherePosition()
-{
-if (cs_) delete cs_;
-}
-
-UnitVector SpherePosition::Transform(const SphereCoordSys& /*cs*/) const
-{
-return(*this);
-}
-
-double SpherePosition::Separation(const SpherePosition& gamma) const 
-{
-return(this->SepAngle(gamma));
-}
Index: trunk/SophyaLib/Samba/spherepos.h
===================================================================
--- trunk/SophyaLib/Samba/spherepos.h	(revision 757)
+++ 	(revision )
@@ -1,82 +1,0 @@
-// This may look like C code, but it is really -*- C++ -*-
-#ifndef SPHEREPOS_H_SEEN
-#define SPHEREPOS_H_SEEN
-
-#include "machdefs.h"
-#include "ppersist.h"
-#include <string>
-#include "unitvector.h"
-#include "anydataobj.h"
-
-static char *head_spherepos_h_ = "$Header: /Users/garnier/temp/CVSSophya/SophyaLib/Samba/spherepos.h,v 1.4 2000-01-04 14:27:11 ansari Exp $";
-
-namespace SOPHYA {
-
-enum SphereCoordSysIds { SphereCoordSys_NEUTRAL =0, 
-                         SphereCoordSys_ROTATION =1,
-                         SphereCoordSys_OTHER = 0xFFFF
-};
-
- 
-class SphereCoordSys : public  AnyDataObj {
-public:
-  
-  SphereCoordSys();
-  SphereCoordSys(const SphereCoordSys& a);
-                 SphereCoordSys(int id, const string& description);
-  virtual        ~SphereCoordSys();
-
-  virtual int    Id() const {return id_; }
-  virtual string description() const {return description_; }
-  SphereCoordSys& operator = (const SphereCoordSys& a);
-
-
-private:
-  int id_;
-  string description_;
-};
-
-// Classe pour la gestion de persistance
-
-class FIO_SphereCoordSys : public  PPersist  {
-public:
-  FIO_SphereCoordSys();
-  FIO_SphereCoordSys(string const & filename);
-  FIO_SphereCoordSys(const SphereCoordSys & obj);
-  FIO_SphereCoordSys(SphereCoordSys * obj);
-  virtual ~FIO_SphereCoordSys();
-  virtual AnyDataObj* DataObj();
-  inline operator SphereCoordSys() { return(*dobj); }
-protected :
-  virtual void       ReadSelf(PInPersist&);
-  virtual void       WriteSelf(POutPersist&) const;  
-  SphereCoordSys * dobj;
-  bool ownobj;
-};
-
-
-class SpherePosition : public UnitVector {
-public:
-  
-                 SpherePosition();
-// cs should be created by new SphereCoordSys...
-// It will be deleted by the destructor
-                 SpherePosition(double theta, double phi, SphereCoordSys* cs = NULL);
-                 SpherePosition(double x, double y, double z, SphereCoordSys* cs = NULL);
-                 SpherePosition(const UnitVector& v, SphereCoordSys* cs = NULL);
-  virtual        ~SpherePosition();
-
-  inline double  phi() const {return Phi(); }
-  inline double  theta() const { return(Theta()); }
-  inline SphereCoordSys *   GetCoordSys() const { return(cs_); } 
-  virtual UnitVector Transform(const SphereCoordSys& cs) const ;
-
-  virtual double     Separation(const SpherePosition& gamma) const;
-
-protected:
-  SphereCoordSys* cs_;
-};
-
-} // Fin du namespace
-
-#endif
Index: trunk/SophyaLib/Samba/spherethetaphi.cc
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 757)
+++ 	(revision )
@@ -1,815 +1,0 @@
-#include "spherethetaphi.h"
-#include "nbmath.h"
-#include <complex>
-#include "piocmplx.h"
-#include <iostream.h>
-
-
-//***************************************************************
-//++ 
-// Class	SphereThetaPhi
-// 
-//
-// include      spherethetaphi.h nbmath.h 
-//
-//    sphere splitted with respect to theta, phi : each hemisphere is 
-//    splitted into (m-1) parallels (equator does not enter into account).
-//    This operation defines m slices, each of which is splitted into 
-//    equidistant meridians. This splitting is realized in such a way that 
-//    all pixels have the same area and are as square as possible.
-
-//    One begins with the hemisphere with positive z, starting from the pole
-//    toward the equator. The first pixel is the polar cap ; it is circular
-//    and centered on theta=0. 
-//--
-//++
-//
-// Links	Parents
-//
-//    SphericalMap 
-//-- 
-
-/* --Methode-- */
-//++
-// Titre	Constructors
-//--
-//++
-
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi()
-
-//--
-{
-  InitNul();
-  pixels_.Reset();
-}
-
-
-/* --Methode-- */
-
-//++
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi(int_4 m)
-
-//    m is the number of slices in theta on an hemisphere (the polar cap
-//    forms the first slice).
-//    pet is a dummy parameter at the moment. 
-//--
-{
-  InitNul();
-  Pixelize(m);
-}
-
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s, bool share)
-  : pixels_(s.pixels_ , share)
-{
-  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
-  NTheta_= s.NTheta_;
-  NPix_  = s.NPix_;
-  NPhi_.ReSize(NTheta_);
-  Theta_.ReSize(NTheta_+1);
-  TNphi_.ReSize(NTheta_+1);
-
-  NPhi_  = s.NPhi_;
-  TNphi_ = s.TNphi_;
-  Theta_ = s.Theta_;
-  Omega_ = s.Omega_;
-}
-
-//++
-// Titre	Destructor
-//--
-//++
-template <class T>
-SphereThetaPhi<T>::~SphereThetaPhi()
-
-//--
-{}
-
-//++
-// Titre	Public Méthods
-//--
-template <class T>
-void SphereThetaPhi<T>::InitNul()
-//
-{
-  NTheta_= 0;
-  NPix_  = 0;
-//  pixels_.Reset();  Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
-}
-
-
-//++
-template <class T>
-void SphereThetaPhi<T>::Resize(int_4 m)
-//   re-pixelize the sphere
-//--
-{
-  InitNul();
-  Pixelize(m);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int_4 SphereThetaPhi<T>::NbPixels() const
-
-//    Return total number of pixels 
-//--
-{
-  return(NPix_);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-T& SphereThetaPhi<T>::PixVal(int_4 k)
-
-//    Return value of pixel with index k
-//--
-{
-  if((k < 0) || (k >= NPix_)) 
-    {
-      throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
-    }
-  return pixels_(k);
-}
-
-//++
-template <class T>
-T const& SphereThetaPhi<T>::PixVal(int_4 k) const
-
-//    Return value of pixel with index k
-//--
-{
-  if((k < 0) || (k >= NPix_)) 
-    {
-       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
-    } 
-  return *(pixels_.Data()+k);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-bool SphereThetaPhi<T>::ContainsSph(double /*theta*/, double /*phi*/) const
-//--
-{
-  return(true);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int_4 SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
-
-//    Return index of the pixel corresponding to 
-//    direction (theta, phi).
-//--
-{
-  double dphi;
-  int i,k;
-  bool fgzn = false;
-
-  if((theta > Pi) || (theta < 0.)) return(-1);
-  if((phi < 0.) || (phi > DeuxPi)) return(-1);
-  if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
-  
-  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
-  // Theta_[kt] et Theta_[kt+1]
-  for( i=1; i< NTheta_; i++ ) 
-    if( theta < Theta_(i) ) break;
-  
-  dphi= DeuxPi/(double)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-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
-
-//   Return (theta,phi) coordinates of middle of  pixel with  index k
-//--
-{
-  int i; 
-  bool fgzn = false;
-  
-  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< NTheta_; i++ ) 
-    if( k < TNphi_(i+1) ) break;
-
-  // angle theta
-  theta= 0.5*(Theta_(i)+Theta_(i+1));
-  if (fgzn) theta= Pi-theta;
-  
-  // angle phi
-  k -= TNphi_(i);
-  phi= DeuxPi/(double)NPhi_(i)*(double)(k+.5);
-  if (fgzn) phi= DeuxPi-phi;
-}
-
-template <class T>
-T SphereThetaPhi<T>::SetPixels(T v)
-{
-pixels_.Reset(v);
-return(v);
-}
-
-//++
-template <class T>
-double SphereThetaPhi<T>::PixSolAngle(int_4 /*dummy*/) const
-
-//    Pixel Solid angle  (steradians)
-//    All the pixels have the same solid angle. The dummy argument is
-//    for compatibility with eventual pixelizations which would not 
-//    fulfil this requirement.
-//--
-{
-  return Omega_;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Limits(int_4 k,double& tetMin,double& tetMax,double& phiMin,double& phiMax) 
-
-//   Return values of theta,phi which limit the pixel with  index k
-//--
-  {
-  int j;
-  double dphi;
-  bool fgzn= false;
-  
-  if((k < 0) || (k >= NPix_)) {
-    tetMin= -99999.; 
-    phiMin= -99999.;  
-    tetMax= -99999.; 
-    phiMax= -99999.;  
-    return;
-  }
-  
-  // si on se trouve dans l'hémisphère Sud
-  if(k >= NPix_/2) { 
-    fgzn= true; 
-    k= NPix_-1-k; 
-  }
-  
-  // on recupere l'indice i de la tranche qui contient le pixel k
-  int i;
-  for( i=0; i< NTheta_; i++ ) 
-    if(k < TNphi_(i+1)) break;
-  
-  // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= Theta_(i);
-  tetMax= Theta_(i+1);
-  // valeurs limites de theta dans l'hemisphere Sud
-  if (fgzn) {
-    tetMin= Pi - Theta_(i+1);
-    tetMax= Pi - Theta_(i);
-  }
-  
-  // pixel correspondant dans l'hemisphere Nord
-  if (fgzn) k= TNphi_(i+1)-k+TNphi_(i)-1;
-  
-  // indice j de discretisation ( phi= j*dphi )
-  j= k-TNphi_(i);
-  dphi= DeuxPi/(double)NPhi_(i);
-  
-  // valeurs limites de phi 
-  phiMin= dphi*(double)(j);
-  phiMax= dphi*(double)(j+1);
-  return;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-uint_4 SphereThetaPhi<T>::NbThetaSlices() const
-
-//    Return number of theta-slices on the sphere
-//--
-{
-  if (NTheta_<=0) 
-    {
-      throw PException(" sphere not pixelized, NbSlice=0 ");
-    }
-  return( 2*NTheta_);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int_4 SphereThetaPhi<T>::NPhi(int_4 kt) const
-
-//    Return number of pixels in phi-direction of the kt-th slice
-//-- 
-{
-  int nbpix;  
-  // verification
-  if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) {
-    kt= 2*NTheta_-1-kt; 
-  }
-  
-  // nombre de pixels
-  nbpix= NPhi_(kt);
-  return(nbpix);
-}
-
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Theta(int_4 kt,double& tetMin,double& tetMax) 
-
-//    Return  theta values which limit the slice kt
-//--
-{
-  bool fgzn= false;
-  // verification
-  if( (kt< 0) || (kt>= 2*NTheta_) ) {
-    tetMin= -99999.; 
-    tetMax= -99999.; 
-    return;
-  }
-
-  // si on se trouve dans l'hemisphere Sud
-  if( kt >= NTheta_ ) {
-    fgzn= true;
-    kt= 2*NTheta_-1-kt; 
-  }
-  
-  // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= Theta_(kt);
-  tetMax= Theta_(kt+1);
-  // valeurs limites de theta dans l'hemisphere Sud
-  if (fgzn) {
-    tetMin= Pi - Theta_(kt+1);
-    tetMax= Pi - Theta_(kt);
-  }  
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Phi(int_4 kt,int_4 jp,double& phiMin,double& phiMax) 
-
-//   Return values of phi which limit the jp-th pixel of the kt-th slice
-//--
-{
-  // verification
-  if((kt < 0) || (kt >= 2*NTheta_)) {
-    phiMin= -99999.; 
-    phiMax= -99999.; 
-    return;
-  }
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
-  
-  // verifie si la tranche kt contient au moins jp pixels
-  if( (jp< 0) || (jp >= NPhi_(kt)) ) {
-    phiMin= -88888.; 
-    phiMax= -88888.; 
-    return;
-  }
-  
-  double dphi= DeuxPi/(double)NPhi_(kt);
-  phiMin= dphi*(double)(jp);
-  phiMax= dphi*(double)(jp+1);
-  return;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int_4 SphereThetaPhi<T>::Index(int_4 kt,int_4 jp) const
-
-//    Return pixel index  with sequence index jp in the slice kt 
-//--
-{
-  int k;
-  bool fgzn= false;
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) {
-    fgzn= true; 
-    kt= 2*NTheta_-1-kt;
-  }
-  
-  // si la tranche kt contient au moins i pixels
-  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);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 
-
-//    Return indices kt (theta) and jp (phi) of  pixel with index k
-//--
-{
-  bool fgzn= false;  
-  // si on se trouve dans l'hemisphere Sud
-  if(k >= NPix_/2) 
-    { 
-      fgzn= true; 
-      k= NPix_-1-k; 
-    }
-  
-  // on recupere l'indice kt de la tranche qui contient le pixel k
-  int i;
-  for(i = 0; i < NTheta_; i++) 
-    if(k < TNphi_(i+1)) break;
-  
-  // indice  kt de tranche
-  if (fgzn) kt= 2*NTheta_-1-i;
-  else kt= i;
-  
-  // indice jp de pixel
-  if (fgzn) jp= TNphi_(i+1)-k-1;
-  else jp= k-TNphi_(i);  
-}
-//++
-template <class T>
-void  SphereThetaPhi<T>::Pixelize(int_4 m) 
-
-//    achieve the splitting into pixels (m has the same signification 
-//    as for the constructor)
-//
-//    Each theta-slice of the north hemisphere will be spitted starting f
-//    from  phi=0 ...
-//
-//    South hemisphere is scanned in the same direction according to phi
-//    and from equator to the pole (the pixel following the last one of
-//    the slice closest to the equator with z>0, is the pixel with lowest 
-//    phi of the slice closest of the equator with z<0).
-//--
-{
-  int ntotpix,j;
-  
-  // Decodage et controle des arguments d'appel :
-  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
-  if (m < 2) m = 2;
-  if (m > 16384) m = 16384;
-  
-  // On memorise les arguments d'appel
-  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...)
-  //  Theta_= new double[m+1];
-  Theta_.ReSize(m+1);
-  
-  // Le nombre de pixels en phi de chacune des bandes en theta
-  //  NPhi_ = new int_4[m];
-  NPhi_.ReSize(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)
-  //  TNphi_= new int_4[m+1];
-    TNphi_.ReSize(m+1);
-  
-  // Calcul du nombre total de pixels dans chaque bande pour optimiser
-  // le rapport largeur/hauteur des pixels 
-  
-  //calotte polaire
-  TNphi_(0)= 0;
-  NPhi_(0) = 1;
-  
-  //bandes jusqu'a l'equateur
-  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  = TNphi_(m-1)+NPhi_(m-1);
-  TNphi_(m)= ntotpix;
-  // et sur la sphere entiere
-  NPix_= 2*ntotpix;
-  
-  // Creation et initialisation du vecteur des contenus des pixels 
-  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 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 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_8& theta, TVector<r_8>& phi, TVector<T>& value) const 
-
-//    For a theta-slice with index 'index', return : 
-//    the corresponding "theta" 
-//    a vector containing the phi's of the pixels of the slice
-//    a vector containing the corresponding values of pixels 
-//--
-
-{
-
-  if(index < 0 || index >= NbThetaSlices()) 
-    {
-       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range"); 
-    }
-
-  int iring= Index(index,0);
-  int lring  = NPhi(index);
-
-  phi.ReSize(lring);
-  value.ReSize(lring);
-  double Te= 0.;
-  double 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>::GetThetaSlice(int_4 index,r_8& theta, r_8& phi0,TVector<int_4>& pixelIndices, TVector<T>& value) const 
-
-//    For a theta-slice with index 'index', return : 
-//    the corresponding "theta" 
-//    the corresponding "phi" for first pixel of the slice 
-//    a vector containing the indices of the pixels of the slice
-//    (equally distributed in phi)
-//    a vector containing the corresponding values of pixels 
-//--
-
-{
-
-  if(index < 0 || index >= NbThetaSlices()) 
-    {
-       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range"); 
- 
-    }
-
-  int iring= Index(index,0);
-  int lring  = NPhi(index);
-
-  pixelIndices.ReSize(lring);
-  value.ReSize(lring);
-  double Te= 0.;
-  double Fi= 0.;
-  for(int kk = 0; kk < lring; kk++) 
-    {
-      pixelIndices(kk)=kk+iring ;
-      value(kk)= PixVal(kk+iring);
-    }
-      PixThetaPhi(iring,theta,phi0);
-}
-
-
-
-
-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_ : ";
-  int i;
-  for(i=0; i < NTheta_; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << NPhi_(i) <<", ";
-    }
-  os << endl;
-
-  os << " contenu de Theta_ : ";
-  for(i=0; i < NTheta_+1; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << Theta_(i) <<", ";
-    }
-  os << endl;
-
-  os << " contenu de TNphi_ : ";
-  for(i=0; i < NTheta_+1; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << TNphi_(i) <<", ";
-    }
-  os << endl;
-
-  os << " contenu de pixels : ";
-  for(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>;
-  dobj->pixels_.SetTemp(true);
-  ownobj= true;
-  Read(filename);
-}
-
-template <class T>
-FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
-{
-  dobj= new SphereThetaPhi<T>(obj, true);
-  dobj->pixels_.SetTemp(true);
-  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)
-{
-
-  if(dobj == NULL) 
-    {
-      dobj= new SphereThetaPhi<T>;
-      dobj->pixels_.SetTemp(true);
-      ownobj= true;      
-    }
-
-  // On lit les 3 premiers uint_8
-  uint_8 itab[3];
-  is.Get(itab, 3);
-// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Read(is);
-
-  // 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);  
-  int_4 mNPix;
-  is.GetI4(mNPix);
-  double mOmeg;
-  is.GetR8(mOmeg);
-  dobj->setParameters(mNPix, mOmeg, mNTheta); 
-  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
-  fio_nd_nphi.Read(is);
-  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
-  fio_nd_Tnphi.Read(is);
-  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
-  fio_nd_Theta.Read(is);
-
-FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-fio_nd.Read(is);
-}
-
-template <class T>
-void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
-{
-
-  if(dobj == NULL) 
-    {
-      cout << " WriteSelf:: dobj= null " << endl;
-      return;
-    }
-
-//  On ecrit 3 uint_8 
-//  0 : Numero de version,  1 : Size index,  2  reserve a l
-uint_8 itab[3];
-itab[0] = 1;
-itab[1] = dobj->SizeIndex();
-itab[2] = 0;
-os.Put(itab, 3);
-// Let's write the SphereCoordSys object
-  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
-  fio_scs.Write(os);
-  
-  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());
-  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
-  fio_nd_nphi.Write(os);
-  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
-  fio_nd_Tnphi.Write(os);
-  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
-  fio_nd_Theta.Write(os);
-FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
-fio_nd.Write(os);
-}
-
-#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
Index: trunk/SophyaLib/Samba/spherethetaphi.h
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.h	(revision 757)
+++ 	(revision )
@@ -1,207 +1,0 @@
-#ifndef SPHERETHETAPHI_SEEN
-#define SPHERETHETAPHI_SEEN
-
-#include "sphericalmap.h"
-#include "ndatablock.h"
-#include "tvector.h"
-
-#include "anydataobj.h"
-#include "ppersist.h"
-
-template <class T>
-class FIO_SphereThetaPhi;  
-
-
-// ***************** Class SphereThetaPhi *****************************
-/*!    sphere splitted with respect to theta, phi : each hemisphere is 
-    splitted into (m-1) parallels (equator does not enter into account).
-   This operation defines m slices, each of which is splitted into 
-    equidistant meridians. This splitting is realized in such a way that 
-    all pixels have the same area and are as square as possible.
-
-    One begins with the hemisphere with positive z, starting from the pole
-    toward the equator. The first pixel is the polar cap ; it is circular
-    and centered on theta=0. 
-*/
-template <class T>
-class SphereThetaPhi : public SphericalMap<T>
-{
-
-  friend class FIO_SphereThetaPhi<T>;
-
-public :
-
-SphereThetaPhi();
-/*!    m is the number of slices in theta on an hemisphere (the polar cap
-   forms the first slice).
-    pet is a dummy parameter at the moment. 
-*/
-SphereThetaPhi(int_4 m);
-SphereThetaPhi(const SphereThetaPhi<T>& s, bool share=false);
-virtual ~SphereThetaPhi();
-
-/*! Setting blockdata to temporary (see ndatablock documentation) */
-inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};
-
-// ------------ Definition of PixelMap abstract methods -
-
-/* retourne le nombre de pixels */ 
-/*!    Return total number of pixels  */
-virtual int_4 NbPixels() const;
-
-/* retourne la valeur du pixel d'indice k */ 
-/*!    Return value of pixel with index k */
-virtual T&       PixVal(int_4 k);
-virtual T const& PixVal(int_4 k) const;
-
-/* Return true if teta,phi in map  */
-virtual bool ContainsSph(double theta, double phi) const;
-/* retourne l'indice du pixel a (theta,phi) */ 
-/*    Return index of the pixel corresponding to direction (theta, phi). */
-virtual int_4 PixIndexSph(double theta, double phi) const;
-
-/* retourne les coordonnees Spheriques du centre du pixel d'indice k */ 
-/*!   Return (theta,phi) coordinates of middle of  pixel with  index k */
-virtual void PixThetaPhi(int_4 k, double& theta, double& phi) const;
-
-/*! Setting pixel values to a constant */
-virtual T SetPixels(T v);
-
-/* retourne/fixe l'angle Solide de Pixel   (steradians) */ 
-/*!   Pixel Solid angle  (steradians)
-
-    All the pixels have the same solid angle. The dummy argument is
-   for compatibility with eventual pixelizations which would not 
-    fulfil this requirement.
-*/
-virtual double PixSolAngle(int_4 dummy=0) const;
- 
-/* retourne/fixe la valeur du parametre de decoupage m */ 
-inline virtual int_4 SizeIndex() const { return( NTheta_); }
-
-/*  Acces to the DataBlock  */
-inline       NDataBlock<T>& DataBlock()       {return pixels_;}
-inline const NDataBlock<T>& DataBlock() const {return pixels_;}
-
-// ------------- Specific methods  ----------------------
-
-/*!   re-pixelize the sphere */
-virtual void Resize(int_4 m);
-
-inline virtual char* TypeOfMap() const {return "TETAFI";};
-
-/* Valeurs de theta des paralleles et phi des meridiens limitant le pixel d'indice k */
-/*   Return values of theta,phi which limit the pixel with  index k */
-virtual void Limits(int_4 k,double& th1,double& th2,double& phi1,double& phi2);
-
-/* Nombre de tranches en theta */ 
-/*!    Return number of theta-slices on the sphere */
-uint_4 NbThetaSlices() 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  */ 
-/*!    Return  theta values which limit the slice kt */
-void Theta(int_4 kt, double& t1, double& t2);
-
-/* Renvoie dans p1,p2 les valeurs phimin et phimax du pixel d'indice jp  */
-/* dans la tranche d'indice kt  */ 
-/*!   Return values of phi which limit the jp-th pixel of the kt-th slice */
-void Phi(int_4 kt, int_4 jp, double& p1, double& p2);
-
-/* Renvoie l'indice k du pixel d'indice jp dans la tranche d'indice kt   */
-/*!   Return pixel index  with sequence index jp in the slice kt */
-int_4 Index(int_4 kt, int_4 jp) const;
-
-/* Indice kt de la tranche et indice jp du pixel d'indice k  */ 
-/*!    Return indices kt (theta) and jp (phi) of  pixel with index k */
-void ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp);
-
-/*!    achieve the splitting into pixels (m has the same signification 
-    as for the constructor)
-
-    Each theta-slice of the north hemisphere will be spitted starting f
-    from  phi=0 ...
-
-    South hemisphere is scanned in the same direction according to phi
-    and from equator to the pole (the pixel following the last one of
-    the slice closest to the equator with z>0, is the pixel with lowest 
-    phi of the slice closest of the equator with z<0).
-*/
-void Pixelize(int_4); 
-
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-    a vector containing the phi's of the pixels of the slice
-
-    a vector containing the corresponding values of pixels 
-*/
-void GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const; 
-
-/*!   For a theta-slice with index 'index', return : 
-
-   the corresponding "theta" 
-
-   the corresponding "phi" for first pixel of the slice 
-
-    a vector containing indices of the pixels of the slice
-
-   (equally distributed in phi)
-
-    a vector containing the corresponding values of pixels 
-*/
- void GetThetaSlice(int_4 index, r_8& theta, r_8& phi0,TVector<int_4>& pixelIndices, TVector<T>& value) const ;
-
-
-/* impression */ 
-void print(ostream& os) const;
-
-private :
-
-// ------------- méthodes internes ----------------------          
-void InitNul();
-inline void setParameters(int nbpix, double omega, int nbThetaIndex) 
-  {
-    NPix_= nbpix;
-    Omega_= omega;
-    NTheta_= nbThetaIndex;
-  }
-
-// ------------- variables internes ---------------------
- int_4 NTheta_;   // nombre de tranches en theta, pour une demi-sphere
- int_4 NPix_;     // nombre total de pixels
- double Omega_; // angle solide constant pour chaque pixel
- NDataBlock<int_4>  NPhi_;  // tableau donnant, pour chaque bande en theta, 
-                         //le nombre de pixels selon phi
- NDataBlock<int_4> TNphi_;
- NDataBlock<r_8> 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); }
-
-protected :
-
-virtual void ReadSelf(PInPersist&);           
-virtual void WriteSelf(POutPersist&) const;  
-SphereThetaPhi<T>* dobj;
-bool ownobj;
-};
-
-#endif
Index: trunk/SophyaLib/Samba/sphericalmap.h
===================================================================
--- trunk/SophyaLib/Samba/sphericalmap.h	(revision 757)
+++ 	(revision )
@@ -1,47 +1,0 @@
-#ifndef SPHERICALMAP_SEEN
-#define SPHERICALMAP_SEEN
-
-// valeurs de Pi, 2*Pi, etc                                              
-#include "nbmath.h"
-#include <math.h>
-#include "pixelmap.h"
-#include "tvector.h" 
-
-// Map of pixels on a whole sphere.
-// Class hierarchy :
-//   PixelMap
-//      SphericalMap
-//         SphereThetaPhi
-//         SphereGorski
-//         SphereIco
-//      LocalMap
-
-template<class T>
-class SphericalMap : public PixelMap<T> 
-{
-  
-public :
-
-SphericalMap() {};
-virtual	 ~SphericalMap() {};
-
-// Overloading of () to access pixel number k.
-inline T& operator()(int k) {return(PixVal(k));}
-inline T  const& operator()(int k) const {return(PixVal(k));}
-inline T& operator()(double theta,double phi) {return(PixValSph(theta,phi));};
-inline T  const& operator()(double theta,double phi) const {return(PixValSph(theta,phi));};
-
-// index characterizing the size pixelization : m for SphereThetaPhi
-// nside for Gorski sphere...
-virtual void Resize(int m)=0;
-virtual uint_4 NbThetaSlices() const=0;
-virtual void GetThetaSlice(int_4 index,r_8& theta, TVector<r_8>& phi, TVector<T>& value) const=0; 
-virtual void GetThetaSlice(int_4 sliceIndex, r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const=0 ;
-};
-#endif
-
-
-
-
-
-
Index: trunk/SophyaLib/Samba/unitvector.cc
===================================================================
--- trunk/SophyaLib/Samba/unitvector.cc	(revision 757)
+++ 	(revision )
@@ -1,102 +1,0 @@
-#include <math.h>
-#include "unitvector.h"
-//*****************************************************************************
-//++ 
-// Class	UnitVector
-//
-// include      unitvector.h vector3d.h math.h
-//--
-//++
-//
-// Links	Parents
-//
-//    Vector3d
-//    
-//-- 
-//++
-// Titre	Constructors
-//--
-//++
-UnitVector::UnitVector() : Vector3d(1.,0.,0.) 
-//
-//--
-{
-}
-//++
-UnitVector::UnitVector(double x, double y, double z) : Vector3d(x,y,z) 
-//
-//--
-{
-  this->Normalize();
-}
-//++
-UnitVector::UnitVector(double theta, double phi) : Vector3d(theta,phi) 
-//
-//--
-{
-  this->Normalize();
-}
-//++
-UnitVector::UnitVector(const Vector3d& v) : Vector3d(v) 
-//
-//--
-{
-  this->Normalize();
-}
-//++
-// Titre	Public Method
-//--
-//++
-void UnitVector::Print(ostream& os) const 
-//
-//--
-{
-  os << "UnitVector : (X,Y,Z)= (" << _x << ", " << _y << ", " << _z 
-     << ") Theta,Phi= " << _theta << ", " << _phi << "\n"
-     << "norme =" << this->Norm() << endl;
-}
-//++
-// Titre	Operators
-//--
-//++
-Vector3d& UnitVector::operator = (const Vector3d& v)
-//
-//--
-
-{
-  Setxyz(v.X(),v.Y(),v.Z());
-  this->Normalize();
-  return *this;
-}
-//++
-Vector3d& UnitVector::operator += (const Vector3d& v)
-//
-//--
-{
-  Setxyz(_x+v.X(),_y+v.Y(),_z+v.Z());
-  this->Normalize();
-  return *this;
-}
-//++
-Vector3d& UnitVector::operator -= (const Vector3d& v)
-//
-//--
-{
-  Setxyz(_x-v.X(),_y-v.Y(),_z-v.Z());
-  this->Normalize();
-  return *this;
-}
-//++
-Vector3d UnitVector::operator + (const Vector3d& v) const
-//
-//--
-{
-  return UnitVector(_x+v.X(),_y+v.Y(),_z+v.Z());
-}
-//++
-Vector3d UnitVector::operator - (const Vector3d& v) const
-//
-//--
-{
-  return UnitVector(_x-v.X(),_y-v.Y(),_z-v.Z());
-}
Index: trunk/SophyaLib/Samba/unitvector.h
===================================================================
--- trunk/SophyaLib/Samba/unitvector.h	(revision 757)
+++ 	(revision )
@@ -1,35 +1,0 @@
-#ifndef UNITVECTOR_H_SEEN
-#define UNITVECTOR_H_SEEN
-
-#include <iostream.h>
-#include <stdio.h>
-#include <string.h>
-#include "vector3d.h"
-
-class UnitVector : public Vector3d 
-{
-  
- public:
-  
-  UnitVector();
-  UnitVector(double x, double y, double z);
-  UnitVector(double theta, double phi=0);
-  UnitVector(const Vector3d&);
-  
-  /*virtual*/ Vector3d& operator=(const Vector3d&);  //$CHECK$ EA 101299 operator = non virtual
-  virtual Vector3d& operator+=(const Vector3d&);
-  virtual Vector3d& operator-=(const Vector3d&);
-  virtual Vector3d operator+(const Vector3d&) const;
-  virtual Vector3d operator-(const Vector3d&) const;
-
-  virtual Vector3d& operator+=(double d) { return Vector3d::operator+= (d); }
-  virtual Vector3d operator+(double d) const { return Vector3d::operator+ (d); }
-  virtual Vector3d operator-(double d) const { return Vector3d::operator- (d); }
-
-  void Print(ostream& os) const;
-  
-};
-
-#endif
-
-
Index: trunk/SophyaLib/Samba/utilgeom.cc
===================================================================
--- trunk/SophyaLib/Samba/utilgeom.cc	(revision 757)
+++ 	(revision )
@@ -1,125 +1,0 @@
-#include <math.h>
-#include <iostream.h>
-#include "utilgeom.h"
-//++
-// Titre	Some utilitary functions for geometry (utilgeom.h)...
-//--
-//++
-int sign(double d) 
-//
-//--
-{
-  return (d >= 0) - (d < 0);
-}
-//++
-double absd(double d) 
-//
-//--
-{
-  return sqrt(d*d);
-}
-//++
-double mod(double d, double periode) 
-//
-//--
-{
-  if( d >= 0 ) return d-floor(d/periode)*periode;
-  else return d-ceil(d/periode)*periode+periode;
-}
-//++
-void swap(double& d1, double& d2) 
-//
-//--
-{
-  double temp;
-  temp=d2;
-  d2=d1;
-  d1=temp;
-}
-//++
-double min(double d1, double d2) 
-//
-//--
-{
-  if( d1 >= d2 ) return d2;
-  else return d1;
-}
-//++
-double max(double d1, double d2) 
-//
-//--
-{
-  return -min(-d1,-d2);
-}
-//++
-int arrondi(double d) 
-//
-//--
-{
-  return (int)(((d-floor(d)) >= 0.5)*ceil(d)+((d-floor(d)) < 0.5)*floor(d));
-}
-//++
-long rangijd(int nc, int i, int j, int d) 
-//
-//--
-{
-  if( i < j ) return 2*i*(2*nc-i-1)+4*(j-i-1)+d;
-  if( i > j ) return 2*j*(2*nc-j-1)+4*(i-j-1)+d+2;
-  if( i == j ) 
-    {
-      cout << "pas d'intersection entre le cercle " << i 
-	   << " et le cercle " << j << "." << endl;
-      return -1;
-    }
-  else return -1;
-}
-//++
-long rangdiff(int nc, int i, int j, int d) 
-//
-//--
-{
-  if( i == j ) 
-    {
-      cout << "diff n'est pas defini entre le cercle " << i 
-	   << " et le cercle " << j << "." << endl;
-      return 0;
-    }
-  int indm=(int)min(i,j);
-  int indM=(int)max(i,j);
-  return indm*(2*nc-indm-1)+2*(indM-indm-1)+d;
-}
-//++
-long rangik(int NtotEch, int i, int k)
-//
-//--
-{
-  return NtotEch*i+k;
-}
-//++
-long ranghk(int NtotEch, int h, int k) 
-//
-//--
-{
-  return NtotEch*h+k;
-}
-//++
-double scangle(double sinus, double cosinus) 
-//
-//--
-{	// Modifie par D. Y. pour eviter des plantages aux petis angles
-	double eps=1.e-10;
-	if( fabs(1.-sinus*sinus-cosinus*cosinus) > eps) 
-    {	   // On renormalise D.Y.
-	  // cerr << "angle non valide: "<< "sinus = " << sinus << " cosinus = " << cosinus << '\n'; // BUGGG
-      double norm= cosinus*cosinus+sinus*sinus;
-      cosinus=cosinus/norm;
-      sinus=sinus/norm;     
-    }
-  if( fabs(cosinus) >1. ) cosinus=1.*copysign(1.,cosinus);
-  if( fabs(sinus) >1. ) sinus=1.*copysign(1.,sinus);
-/*  if( fabs(sinus) == 0. ) sinus=0.;
-    if( fabs(cosinus) == 0. ) cosinus=0.;
-*/
-  return acos(cosinus)*copysign(1,sinus);
-}
-
Index: trunk/SophyaLib/Samba/utilgeom.h
===================================================================
--- trunk/SophyaLib/Samba/utilgeom.h	(revision 757)
+++ 	(revision )
@@ -1,23 +1,0 @@
-#ifndef MYFUNC_H_SEEN
-#define MYFUNC_H_SEEN
-#include <math.h>
-#ifdef __MWERKS__
- #include "mwerksmath.h"		// D.Y. Portage mac/ Math.h Unix
-#include "unixmac.h"
-#endif
-#include <vector>
-
-const double pi2=2.*M_PI;
-const double dtor=M_PI/180.;
-const double pi_over_2=M_PI/2.;
-
-int sign(double);
-double absd(double);
-double mod(double, double=0);
-void swap(double&, double&);
-double min(double, double);
-double max(double, double);
-int arrondi(double);
-double scangle(double sinus, double cosinus);
-
-#endif
Index: trunk/SophyaLib/Samba/vector3d.cc
===================================================================
--- trunk/SophyaLib/Samba/vector3d.cc	(revision 757)
+++ 	(revision )
@@ -1,411 +1,0 @@
-#include "machdefs.h"
-#include <math.h>
-#include "pexceptions.h"
-#include "vector3d.h"
-#include "utilgeom.h"
-//++ 
-// Class	Vector3d
-//
-// include      vector3d.h utilgeom.h  longlat.h math.h
-//
-//
-//    3-D geometry. 
-//    All computations are made with angles in radians and with spherical
-//    coordinates theta, phi.
-//
-//    Concerning Euler's angles, the reference is :
-// 
-//    "Classical Mechanics" 2nd edition, H. Goldstein, Addison Wesley
-//--
-//++
-// Titre	Constructors
-//--
-//++
-Vector3d::Vector3d()
-//
-//-- 
-{
-  Setxyz(1.,0.,0.);
-}
-//++
-Vector3d::Vector3d(double x, double y, double z) 
-//
-//-- 
-{
-  _x=x;
-  _y=y;
-  _z=z;
-  xyz2ThetaPhi();
-}
-//++
-Vector3d::Vector3d(double theta, double phi) 
-//
-//-- 
-{
-  _theta=mod(theta,M_PI); // dans [0;pi]
-  _phi=mod(phi,pi2); // dans [0;2pi]
-  ThetaPhi2xyz();
-}
-//++
-Vector3d::Vector3d(const LongLat& ll)
-//
-//-- 
-{
-  _theta=ll.Theta(); // dans [0;pi]
-  _phi=ll.Phi(); // dans [0;2pi]
-  ThetaPhi2xyz();
-}
-//++
-Vector3d::Vector3d(const Vector3d& v) 
-//
-//-- 
-{
-  _x=v._x;
-  _y=v._y;
-  _z=v._z;
-  _theta=v._theta;
-  _phi=v._phi;
-}
-//++
-// Titre	Public methods
-//--
-//++
-void Vector3d::SetThetaPhi(double theta, double phi)
-//
-//-- 
-{
-  _theta=mod(theta,M_PI);
-  _phi=mod(phi,pi2);
-  ThetaPhi2xyz();
-}
-//++
-void Vector3d::Setxyz(double x, double y, double z) 
-//
-//-- 
-{
-  _x=x;
-  _y=y;
-  _z=z;
-  xyz2ThetaPhi();
-}
-//++
-void Vector3d::ThetaPhi2xyz()
-//
-//-- 
-{
-  _x=sin(_theta)*cos(_phi);
-  _y=sin(_theta)*sin(_phi);
-  _z=cos(_theta);
-}
-//++
-void Vector3d::xyz2ThetaPhi()
-//
-//-- 
-{
-  double norm=this->Norm();
-  if( norm != 0. ) 
-    {
-      _theta=acos(_z/norm); // dans [0,Pi]
-      if( mod(_theta,M_PI) == 0. ) _phi=0.; // on est sur +-Oz, le vecteur z est en phi=0
-      //      else _phi=acos(_x/sin(_theta)/norm)+M_PI*(_y<0);
-      else _phi=scangle(_y/sin(_theta)/norm,_x/sin(_theta)/norm);
-    }
-  else // vecteur nul
-    {
-      _theta=0.;
-      _phi=0.;
-    }
-}
-//++
-Vector3d& Vector3d::Normalize() 
-//
-//-- 
-{
-  double norm=this->Norm();
-  if( norm != 0. )  (*this)/=norm;
-  else cerr << "Division par zero" << endl;
-  return *this;
-} 
-//++
-double Vector3d::Norm() const 
-//
-//-- 
-{
-  return sqrt(_x*_x+_y*_y+_z*_z);
-}
-//++
-double Vector3d::Psc(const Vector3d& v) const 
-//
-//    dot product
-//-- 
-{
-  return _x*v._x+_y*v._y+_z*v._z;
-}
-//++
-double Vector3d::SepAngle(const Vector3d& v) const 
-//
-//    angular gap between 2 vectors in [0,Pi]
-//-- 
-{
-  double n1=this->Norm();
-  double n2=v.Norm();
-  double ret;
-  if( n1!=0. && n2!=0. ) ret=acos((this->Psc(v))/n1/n2);
-  else 
-    {
-      cerr << "Division par zero" << endl;
-      ret=0.;
-    }
-  return ret;
-}
-//++
-Vector3d Vector3d::Vect(const Vector3d& v) const 
-//
-//    vector product
-//-- 
-{
-  double xo=_y*v._z-_z*v._y;
-  double yo=_z*v._x-_x*v._z;
-  double zo=_x*v._y-_y*v._x;
-  return Vector3d(xo,yo,zo);
-}
-//++
-Vector3d Vector3d::VperpPhi() const 
-//
-//    perpendicular vector, with equal phi
-//-- 
-{
-  double theta;
-  if( _theta != pi_over_2 ) theta=_theta+(0.5-(_theta>pi_over_2))*M_PI; // on tourne theta de +-pi/2
-  else theta=0.;
-  return Vector3d(theta,_phi);
-}
-//++
-Vector3d Vector3d::VperpTheta() const 
-//
-//    perpendicular vector with equal theta
-//-- 
-{ 
-  throw PError("Vector3d::VperpTheta() - Logic Error DY/Reza 20/02/2000"); 
-  // Bug ??? (D. Yvon, Fevrier 2000)
-  //  double phi=mod(_phi+pi_over_2,pi2); // on tourne phi de pi/2
-  //  return Vector3d(_theta,phi);
-}
-
-Vector3d Vector3d::EPhi() const 
-{
-  Vector3d temp;
-  if ( fabs(_z) == 1. ) // si on est en +- Oz
-    {
-      temp=Vector3d(1.,0.,0.);
-    }
-  else 
-    {
-      Vector3d k(0,0,-1);
-      temp=this->Vect(k);
-      temp.Normalize();
-    }
-  return temp;
-}
-//++
-Vector3d Vector3d::ETheta() const 
-//
-//-- 
-{
-  Vector3d temp=this->Vect(EPhi());
-  temp.Normalize();
-  return temp;
-}
-
-//++
-Vector3d Vector3d::Euler(double phi, double theta, double psi) const 
-//
-//    Euler's rotations
-//-- 
-{
-  double cpsi=cos(psi);
-  double ctheta=cos(theta);
-  double cphi=cos(phi);
-  double spsi=sin(psi);
-  double stheta=sin(theta);
-  double sphi=sin(phi);
-  double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
-    +(cpsi*sphi+ctheta*cphi*spsi)*_y
-    +spsi*stheta*_z;
-  double ynew=(-spsi*cphi-ctheta*sphi*cpsi)*_x
-    +(-spsi*sphi+ctheta*cphi*cpsi)*_y
-    +cpsi*stheta*_z;
-  double znew=stheta*sphi*_x-stheta*cphi*_y+ctheta*_z;
-  return Vector3d(xnew,ynew,znew);
-}
-//++
-Vector3d Vector3d::InvEuler(double phi, double theta, double psi) const 
-//
-//    inverse rotation
-//-- 
-{
-  double cpsi=cos(psi);
-  double ctheta=cos(theta);
-  double cphi=cos(phi);
-  double spsi=sin(psi);
-  double stheta=sin(theta);
-  double sphi=sin(phi);
-  double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
-    -(spsi*cphi+ctheta*sphi*cpsi)*_y
-    +sphi*stheta*_z;
-  double ynew=(cpsi*sphi+ctheta*cphi*spsi)*_x
-    +(-spsi*sphi+ctheta*cphi*cpsi)*_y
-    -cphi*stheta*_z;
-  double znew=stheta*spsi*_x+stheta*cpsi*_y+ctheta*_z;
-  return Vector3d(xnew,ynew,znew);
-}
-//++
-Vector3d Vector3d::Rotate(const Vector3d& omega, double phi)
-//
-//    rotation of angle phi around an axis omega (Maxwell's rule)
-//-- 
-{
-  Vector3d rotationaxis=omega;
-  rotationaxis.Normalize();
-  double n=this->Psc(rotationaxis);
-  Vector3d myself=*this;
-  Vector3d rotate=n*rotationaxis+(myself-n*rotationaxis)*cos(phi)-(myself^rotationaxis)*sin(phi);
-  return rotate;
-} 
-//++  
-void Vector3d::Print(ostream& os) const 
-//
-//--
-{
-  os << "Vector3: (X,Y,Z)= (" << _x << ", " << _y << ", " << _z 
-     << ") --- Theta,Phi= " << _theta << ", " << _phi << "\n"
-     << "Norme = " << this->Norm() << endl;
-}
-//++
-// Titre	Operators
-//--
-//++
-Vector3d& Vector3d::operator += (const Vector3d& v) 
-//
-//-- 
-{
-  *this=*this+v;
-  return *this;
-}
-//++
-Vector3d& Vector3d::operator -= (const Vector3d& v) 
-//
-//-- 
-{
-  *this=*this-v;
-  return *this;
-}
-//++
-Vector3d& Vector3d::operator += (double d) 
-//
-//-- 
-{
-  Setxyz(_x+d,_y+d,_z+d);
-  return *this;
-}
-//++
-Vector3d& Vector3d::operator /= (double d) 
-//
-//-- 
-{
-  if( d != 0. ) Setxyz(_x/d,_y/d,_z/d);
-  else cerr << "Division par zero." << endl;
-  return *this;
-}
-//++
-Vector3d& Vector3d::operator *= (double d) 
-//
-//-- 
-{
-  Setxyz(_x*d,_y*d,_z*d);
-  return *this;
-}
-//++
-Vector3d Vector3d::operator ^ (const Vector3d& v) const 
-//
-//    vector product
-//-- 
-{
-  return this->Vect(v);
-}
-//++
-Vector3d Vector3d::operator + (double d) const 
-//
-//-- 
-{
-  return Vector3d(_x+d,_y+d,_z+d);
-}
-//++
-Vector3d Vector3d::operator + (const Vector3d& v) const 
-//
-//-- 
-{
-  return Vector3d(_x+v._x,_y+v._y,_z+v._z);
-}
-//++
-Vector3d Vector3d::operator - (double d) const 
-//
-//-- 
-{
-  return *this+(-d);
-}
-//++
-Vector3d Vector3d::operator - (const Vector3d& v) const 
-//
-//-- 
-{
-  return *this+(v*(-1.));
-}
-//++
-Vector3d Vector3d::operator * (double d) const 
-//
-//-- 
-{
-  return Vector3d(d*_x,d*_y,d*_z);
-}
-//++
-double Vector3d::operator * (const Vector3d& v) const 
-//
-//    dot product
-//-- 
-{
-  return this->Psc(v);
-}
-//++
-Vector3d Vector3d::operator / (double d) const 
-//
-//-- 
-{
-  Vector3d ret=*this;
-  if( d != 0. ) ret/=d;
-  else  cerr << "Division par zero." << endl;
-  return ret;
-}
-//++
-Vector3d& Vector3d::operator = (const Vector3d& v) 
-//
-//-- 
-{
-  if( this != &v ) 
-    {
-      _x=v._x; 
-      _y=v._y; 
-      _z=v._z; 
-      _theta=v._theta; 
-      _phi=v._phi;
-    }
-  return *this;
-}
-//++
-bool Vector3d::operator == (const Vector3d& v) 
-//
-//-- 
-{
-  return (this==&v);
-}
-
Index: trunk/SophyaLib/Samba/vector3d.h
===================================================================
--- trunk/SophyaLib/Samba/vector3d.h	(revision 757)
+++ 	(revision )
@@ -1,146 +1,0 @@
-#ifndef VECTOR3D_H_SEEN
-#define VECTOR3D_H_SEEN
-
-#include <math.h>
-#include <iostream.h>
-#include <stdio.h>
-#include <string.h>
-#ifdef __MWERKS__
-#include "unixmac.h"
-#endif
-#include "longlat.h"
-
-/*
-  Geometrie en dimension 3. 
-  Tous les calculs sont faits en radians 
-  et en coordonnees spheriques theta,phi
-  pour les rotations (angles d'Euler) ma source est 
-  "Classical Mechanics" 2nd edition, H. Goldstein, Addison Wesley
-*/
-/*!    3-D geometry. 
-
-    All computations are made with angles in radians and with spherical
-    coordinates theta, phi.
-
-    Concerning Euler's angles, the reference is :
- 
-    "Classical Mechanics" 2nd edition, H. Goldstein, Addison Wesley
-*/
-
-class Vector3d 
-{
-
- public:
-  
-  Vector3d();
-  Vector3d(double x, double y, double z);
-  Vector3d(double theta, double phi);
-  Vector3d(const LongLat&);
-  Vector3d(const Vector3d&);
-
-//   To manipulate the vector
-  virtual void Setxyz(double x, double y, double z);
-  virtual void SetThetaPhi(double theta,  double phi);
-  virtual void ThetaPhi2xyz();
-  virtual void xyz2ThetaPhi();
-
-// Acces to coordinates
-  inline double Theta() const {return _theta;}
-  inline double Phi() const {return _phi;}
-  inline double X() const {return _x;}
-  inline double Y() const {return _y;}
-  inline double Z() const {return _z;}
-
-  virtual Vector3d& Normalize();
-  virtual double Norm() const;
-
-  // produit scalaire
-  virtual double Psc(const Vector3d&) const;
-
-  // ecart angulaire entre 2 vecteurs dans [0,Pi]
-  /*!   angular gap between 2 vectors in [0,Pi] */
-  virtual double SepAngle(const Vector3d&) const;
-
-  // produit vectoriel
-  /*!    vector product */
-  virtual Vector3d Vect(const Vector3d&) const;
-
-  // vecteur perpendiculaire de meme phi
-  /*!    perpendicular vector, with equal phi */
-  virtual Vector3d VperpPhi() const;
-
-  // vecteur perpendiculaire de meme theta
-  /*!    perpendicular vector with equal theta */
-  virtual Vector3d VperpTheta() const;
-
-  virtual Vector3d ETheta() const;
-  virtual Vector3d EPhi() const;
-
-  // rotations d'Euler
-  /*!    Euler's rotations */
-  // rotations d Euler
-  virtual Vector3d Euler(double, double, double) const;
-
-  // rotation inverse
-  /*!    inverse rotation */
-  Vector3d InvEuler(double, double, double) const;
-
-  // rotation d'angle phi autour d'un axe omega (regle du tire-bouchon)
-  /*!    rotation of angle phi around an axis omega (Maxwell's rule) */
-  Vector3d Rotate(const Vector3d& omega,double phi);
-
-  /*virtual*/ Vector3d& operator=(const Vector3d&); // $CHECK$ EA 101299
-  virtual Vector3d& operator+=(const Vector3d&);
-  virtual Vector3d& operator-=(const Vector3d&);
-  virtual Vector3d operator+(const Vector3d&) const;
-  virtual Vector3d operator-(const Vector3d&) const;
-
-  virtual Vector3d& operator+=(double);
-  virtual Vector3d& operator/=(double);
-  virtual Vector3d& operator*=(double);
-
-  virtual Vector3d operator+(double) const;
-  virtual Vector3d operator-(double) const;
-  virtual Vector3d operator*(double) const;
-  virtual Vector3d operator/(double) const;
-
-  /*!    vector product */
-  virtual Vector3d operator^(const Vector3d&) const; // produit vectoriel
-  /*!    dot product */
-  virtual double operator*(const Vector3d&) const; // produit scalaire
-
-  bool operator==(const Vector3d&);
-  
-  virtual void Print(ostream& os) const;
-
- protected:
-
-  double _x;
-  double _y;
-  double _z;
-  double _theta;
-  double _phi;
-
-};
-
-inline ostream& operator<<(ostream& s, const Vector3d& v) 
-{  
-  v.Print(s);  
-  return s;  
-}
-
-// fonctions globales
-
-inline Vector3d operator*(double d, const Vector3d& v) 
-{
-  return v*d;
-}
-
-inline Vector3d operator+(double d, const Vector3d& v) 
-{
-  return v+d;
-}
-
-#endif
-
-
