Index: trunk/SophyaLib/Samba/Makefile
===================================================================
--- trunk/SophyaLib/Samba/Makefile	(revision 228)
+++ trunk/SophyaLib/Samba/Makefile	(revision 228)
@@ -0,0 +1,93 @@
+include ../Mgr/Makefile.h
+ 
+all: $(LIB)libSamba.a
+clean:
+	rm -f $(OBJ)bruit.o $(OBJ)gorfph.o $(OBJ)localmap.o $(OBJ)mlobe.o $(OBJ)scan.o $(OBJ)spheregorski.o $(OBJ)spherethetaphi.o
+	rm $(LIB)libSamba.a
+$(LIB)libSamba.a : $(OBJ)bruit.o $(OBJ)gorfph.o $(OBJ)localmap.o $(OBJ)mlobe.o $(OBJ)scan.o $(OBJ)spheregorski.o $(OBJ)spherethetaphi.o
+	$(AR) $(ARFLAGS) $@ $?
+ifeq ($(CXX),cxx)
+	$(CXX) $? $(CPPFLAGS) $(CXXFLAGS) -L$(LIB) -lSamba -o $(OBJ)xx.x -Hf
+endif 
+ 
+$(OBJ)bruit.o: bruit.cc $(INC)fmath.h \
+ $(INC)cimage.h \
+ $(INC)rzimage.h \
+ $(INC)rzvect.h \
+ $(INC)machine.h \
+ $(INC)defs.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)perrors.h \
+ $(INC)datatypes.h \
+ $(INC)ppersist.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h bruit.h \
+ $(INC)nbrandom.h \
+ $(INC)hbook.h
+$(OBJ)gorfph.o: gorfph.cc
+$(OBJ)localmap.o: localmap.cc localmap.h pixelmap.h \
+ $(INC)ppersist.h \
+ $(INC)defs.h \
+ $(INC)machine.h \
+ $(INC)perrors.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h sphericalmap.h \
+ $(INC)nbmath.h \
+ $(INC)cvector.h \
+ $(INC)matrix.h \
+ $(INC)peida.h \
+ $(INC)utils.h \
+ $(INC)fmath.h
+$(OBJ)mlobe.o: mlobe.cc mlobe.h sphericalmap.h \
+ $(INC)nbmath.h pixelmap.h \
+ $(INC)ppersist.h \
+ $(INC)perrors.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h \
+ $(INC)cvector.h \
+ $(INC)matrix.h \
+ $(INC)peida.h \
+ $(INC)utils.h \
+ $(INC)fmath.h unitvec.h \
+ $(INC)timing.h
+$(OBJ)scan.o: scan.cc scan.h $(INC)ppersist.h \
+ $(INC)defs.h \
+ $(INC)machine.h \
+ $(INC)perrors.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h \
+ $(INC)nbmath.h
+$(OBJ)spheregorski.o: spheregorski.cc spheregorski.h sphericalmap.h \
+ $(INC)nbmath.h pixelmap.h \
+ $(INC)ppersist.h \
+ $(INC)perrors.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h \
+ $(INC)cvector.h \
+ $(INC)matrix.h \
+ $(INC)peida.h \
+ $(INC)utils.h \
+ $(INC)fmath.h \
+ $(INC)strutil.h
+$(OBJ)spherethetaphi.o: spherethetaphi.cc spherethetaphi.h sphericalmap.h \
+ $(INC)nbmath.h pixelmap.h \
+ $(INC)ppersist.h \
+ $(INC)perrors.h \
+ $(INC)exceptions.h \
+ $(INC)peidainit.h \
+ $(INC)pclassids.h \
+ $(INC)dvlist.h \
+ $(INC)cvector.h \
+ $(INC)matrix.h \
+ $(INC)peida.h \
+ $(INC)utils.h \
+ $(INC)fmath.h pdataarray.h prefcount.h
Index: trunk/SophyaLib/Samba/Samba.o.list
===================================================================
--- trunk/SophyaLib/Samba/Samba.o.list	(revision 228)
+++ trunk/SophyaLib/Samba/Samba.o.list	(revision 228)
@@ -0,0 +1,7 @@
+bruit.o
+gorfph.o
+localmap.o
+mlobe.o
+scan.o
+spheregorski.o
+spherethetaphi.o
Index: trunk/SophyaLib/Samba/bruit.cc
===================================================================
--- trunk/SophyaLib/Samba/bruit.cc	(revision 228)
+++ trunk/SophyaLib/Samba/bruit.cc	(revision 228)
@@ -0,0 +1,457 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "fmath.h"
+#include "cimage.h"
+
+#include "bruit.h"
+#include "nbrandom.h"
+// #include "rancern.h"
+#include "hbook.h"
+
+//  Le code des classes   NoiseGenerator  RWalkNoise
+
+
+/* --Methode-- */
+NoiseGenerator::NoiseGenerator(float sigma)
+{
+if (sigma < 0.)  sigma = 1.;
+mNCoups = 0;
+mSigma = sigma;
+//printf("-- NoiseGenerator::NoiseGenerator(%g) (Constructeur) ---\n", sigma);
+}
+
+/* --Methode-- */
+NoiseGenerator::~NoiseGenerator()
+{
+    //printf("-- NoiseGenerator::~NoiseGenerator() (Destructeur) --- \n");
+}
+
+/* --Methode-- */
+float NoiseGenerator::Noise()
+{
+mNCoups++;
+return(NorRand()*mSigma);
+}
+
+/* --Methode-- */
+RWalkNoise::RWalkNoise(float sigma)
+    : NoiseGenerator(sigma)
+{
+mState = 0.;
+//printf("--  RWalkNoise::RWalkNoise(%g) (Constructeur) ---\n", sigma);
+}
+
+
+/* --Methode-- */
+RWalkNoise::~RWalkNoise()
+{
+//printf("--  RWalkNoise::~RWalkNoise (Destructeur) ---\n");
+
+}
+
+/* --Methode-- */
+float RWalkNoise::Noise()
+{
+mState += NoiseGenerator::Noise();
+return(mState);
+}
+
+
+/* --Methode-- */
+OOFNoise::OOFNoise(float sigma, int typacf, int mem, float tau)
+    : NoiseGenerator(sigma)
+{
+if (typacf != ACF_Exp)   typacf = ACF_Exp;
+// Fonction d'autocorrelation donnant un bruit en 1/f
+//                  
+mTypACF = typacf;
+if (mem < 1)  mem = 1;
+mMemL = mem;
+if (tau < 1.e-6)  tau = 1.e-6;
+mTau = tau;
+mState = new double[mem];
+mDyn = new double[mem];
+int i;
+mDyn[0]=1.;
+mState[0]=0.;
+for(i=1; i<mem; i++) {
+  mState[i] = 0.;   mDyn[i] = 1./sqrt((double)i);
+}
+//printf("--  OOFNoise::OOFNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
+//       sigma, typacf, mem, tau);
+}
+
+
+/* --Methode-- */
+OOFNoise::~OOFNoise()
+{
+delete[] mState;
+delete[] mDyn;
+// printf("--  OOFNoise::~OOFNoise (Destructeur) ---\n");
+}
+
+
+/* --Methode-- */
+float OOFNoise::Noise()
+{
+int i;
+double rn;
+
+for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
+mState[0] = NoiseGenerator::Noise();
+rn = 0.;
+for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
+return(rn);
+}
+
+/* --Methode-- */
+void OOFNoise::Print()
+{
+int i,j;
+printf("OOFNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
+printf("OOFNoise::Print() Vecteur de Dynamique / Etat : \n");
+for(i=0; i<mMemL-1; i+=2) 
+  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
+         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
+return;
+}
+
+/* --Methode-- */
+EXPNoise::EXPNoise(float sigma, int typacf, int mem, float tau)
+    : NoiseGenerator(sigma)
+{
+if (typacf != ACF_Exp)   typacf = ACF_Exp;
+// Fonction d'autocorrelation exponentiel disponible uniquement
+//                  Reza , Juillet 97
+mTypACF = typacf;
+if (mem < 1)  mem = 1;
+mMemL = mem;
+if (tau < 1.e-6)  tau = 1.e-6;
+mTau = tau;
+mState = new double[mem];
+mDyn = new double[mem];
+int i;
+for(i=0; i<mem; i++) {
+  mState[i] = 0.;   mDyn[i] = exp(-(double)i/tau);
+}
+//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
+//       sigma, typacf, mem, tau);
+}
+
+
+/* --Methode-- */
+EXPNoise::~EXPNoise()
+{
+delete[] mState;
+delete[] mDyn;
+//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
+}
+
+
+/* --Methode-- */
+float EXPNoise::Noise()
+{
+int i;
+double rn;
+
+for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
+mState[0] = NoiseGenerator::Noise();
+rn = 0.;
+for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
+return(rn);
+}
+
+/* --Methode-- */
+void EXPNoise::Print()
+{
+int i,j;
+printf("EXPNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
+printf("EXPNoise::Print() Vecteur de Dynamique / Etat : \n");
+for(i=0; i<mMemL-1; i+=2) 
+  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
+         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
+
+return;
+}
+/* --Methode-- */
+MemNoise::MemNoise(float sigma, int mem, float tau, int ava)
+    : NoiseGenerator(sqrtf(fabsf(tau))*sigma)
+{
+/* on tire les instants des impulsions successives selon une
+loi de distance exponentielle de duree tau
+*/
+const int place = 1000;
+mTypACF = ACF_Exp;
+if (mem < 1)  mem = 1;
+mMemL = place;
+mCoupHaut=mem;
+if (tau < 1.)  tau = 1.;
+mTau = tau;
+mStPos = new float[place];
+mStNeg = new float[place];
+mTePos = new float[place];
+mTeNeg = new float[place];
+// printf("apres reservation\n");
+int i;
+for(i=0; i<place; i++)
+{
+mStNeg[i] = 0.; mStPos[i] = 0.;
+}
+// printf("apres initialisation tableaux \n");
+
+mNappel=0;
+mNtirage=1;
+mNappLast=0;
+mMemPos=1;
+mMemNeg=1;
+mduree=0.;
+mTePos[0]=0.;
+mTeNeg[0]=0.;
+mTdernier=0.;
+
+Avance(ava);
+// printf("fin du createur MemNoise \n");
+//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
+//       sigma, typacf, mem, tau);
+}
+
+
+/* --Methode-- */
+MemNoise::~MemNoise()
+{
+delete[] mStPos;
+delete[] mStNeg;
+delete[] mTePos;
+delete[] mTeNeg;
+//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
+}
+
+/* --Methode-- */
+float MemNoise::Noise()
+{
+return(Avance(0));
+}
+
+/* --Methode-- */
+float MemNoise::Avance(long Asauter)
+{
+int i,j;
+float boum;
+float rn;
+float Somme,SbaryT,Poids,BaryT,Approx,Erreur;
+int nReduit;
+const float epsilon=.05;
+// printf("entree dans methode Noise de MemNoise \n");
+long mNsaut;
+for (mNsaut = Asauter ; mNsaut >= 0 ; mNsaut--)
+    {  // debut mNsaut
+	// printf("mNsaut = %ld\n",mNsaut);
+mNappel++;
+mTdernier += 1.;
+for(i=0; i<mMemPos; i++) mTePos[i] += 1.;
+for(i=0; i<mMemNeg; i++) mTeNeg[i] += 1.;
+
+if (mduree < mTdernier)
+    {
+	/* la duree depuis le dernier tirage est maintenant ecoulee
+	   on commence par remettre a jour le temps ecoule entre la derniere
+	   impulsion (celle qu'on va tirer maintenant) et l'instant present  */
+	   
+	   mTdernier = mTdernier - mduree;
+        
+	/* puis on tire une nouvelle impulsion de bruit */
+
+	   boum =  NoiseGenerator::Noise();
+
+        /* on decale la memoire d'une case vers le passe 
+	   puis on met a jour la case [0] */
+
+if(boum>0)
+    {
+	for(i=mMemPos; i>0; i--)
+	    {
+		mStPos[i] = mStPos[i-1];
+		mTePos[i] = mTePos[i-1];
+	    }
+	mTePos[0] = mTdernier;
+	mStPos[0] = boum;
+	if (mMemPos<mMemL) mMemPos++;
+/* a chaque tirage, on regroupe les impulsions precedentes */
+for(i=mMemPos-2 ; i>0 ; i--)
+    {
+	if(mStPos[i]==0.)continue;
+	//	printf("i,temps,bruit %d,%g,%g \n",i,mTePos[i],mStPos[i]);
+	Somme = mStPos[i]/sqrtf(mTePos[i]);
+	    SbaryT = mStPos[i] * mTePos[i];
+	    Poids = mStPos[i];
+	    for(j=i-1 ; j>=0 ; j--)
+		{
+		    //		printf("  j,temps,bruit %d,%g,%g \n",j,mTePos[j],mStPos[j]);
+		Somme += mStPos[j]/sqrtf(mTePos[j]);
+		    Poids += mStPos[j];
+		    SbaryT += mStPos[j] * mTePos[j];
+		    BaryT = SbaryT / Poids;
+		    Approx = Poids / sqrtf(BaryT);
+		    Erreur = (Approx - Somme) / Somme;
+		    //		    		    printf("i,j = %d %d Approx = %g \n",i,j,fabsf(Erreur));
+		    //		    printf("Temps,Bruit = %g %g \n",BaryT,Poids);
+		    if(fabsf(Erreur)<epsilon)
+			{
+			    mStPos[i]=Poids;
+			    mTePos[i]=BaryT;
+			    mStPos[j]=0.;
+			}
+		    else break;
+		}
+    }
+nReduit=0;
+for(i=0; i<mMemPos-1 ; i++)
+    {
+	if(mStPos[i]!=0)
+	    {
+		mStPos[nReduit]=mStPos[i];
+		mTePos[nReduit]=mTePos[i];
+		nReduit++;
+	    }
+    }
+mMemPos=nReduit+1;
+
+    }
+else
+    {
+	for(i=mMemNeg; i>0; i--)
+	    {
+		mStNeg[i] = mStNeg[i-1];
+		mTeNeg[i] = mTeNeg[i-1];
+	    }
+	mTeNeg[0] = mTdernier;
+	mStNeg[0] = boum;
+	if (mMemNeg<mMemL) mMemNeg++;
+    
+
+/* idem pour negatifs */
+for(i=mMemNeg-2 ; i>0 ; i--)
+    {
+	if(mStNeg[i]==0.)continue;
+	//	printf("i,temps,bruit %d,%g,%g \n",i,mTeNeg[i],mStNeg[i]);
+	Somme = mStNeg[i]/sqrtf(mTeNeg[i]);
+	    SbaryT = mStNeg[i] * mTeNeg[i];
+	    Poids = mStNeg[i];
+	    for(j=i-1 ; j>=0 ; j--)
+		{
+		    //		printf("  j,temps,bruit %d,%g,%g \n",j,mTeNeg[j],mStNeg[j]);
+		Somme += mStNeg[j]/sqrtf(mTeNeg[j]);
+		    Poids += mStNeg[j];
+		    SbaryT += mStNeg[j] * mTeNeg[j];
+		    BaryT = SbaryT / Poids;
+		    Approx = Poids / sqrtf(BaryT);
+		    Erreur = (Approx - Somme) / Somme;
+		    //		    printf("i,j = %d %d Approx = %g \n",i,j,Erreur);
+		    //		    printf("Temps,Bruit = %g %g \n",BaryT,Poids);
+		    if(fabsf(Erreur)<epsilon)
+			{
+			    mStNeg[i]=Poids;
+			    mTeNeg[i]=BaryT;
+			    mStNeg[j]=0.;
+			}
+		    else break;
+		}
+    }
+nReduit=0;
+for(i=0; i<mMemNeg-1 ; i++)
+    {
+	if(mStNeg[i]!=0)
+	    {
+		mStNeg[nReduit]=mStNeg[i];
+		mTeNeg[nReduit]=mTeNeg[i];
+		nReduit++;
+	    }
+    }
+mMemNeg=nReduit+1;
+    }
+/* on tire la duree a attendre avant la prochaine impulsion de bruit */
+// $CHECK$ avec Francois - rndm() remplace par Reza 10/01/99
+mduree = -mTau*log(frand01());  
+
+// if(mNtirage<mMemL) mNtirage++;
+mNtirage++;
+/*
+int idn=1;
+float pds=1.;
+float nul=0.;
+hfill_(&idn,&mduree,&nul,&pds);
+*/
+// printf("Ntirage %d,mMemPos,Neg %d,%d, Suppress = %d \n",mNtirage,mMemPos,mMemNeg, Suppress);
+    }
+/* On calcule le bruit total */
+
+//printf("mMemPos,Neg,%d %d\n",mMemPos,mMemNeg);
+//for(i=0;i<mMemPos;i++) printf("mTePos %g \n",mTePos);
+//for(i=0;i<mMemNeg;i++) printf("mTeNeg %g \n",mTeNeg);
+if(mNsaut == 0)
+    {
+	// printf("coucou\n");
+rn = 0.;
+int indlim=0;
+for(i=0; i < mMemPos-1; i++) 
+{
+if(mTePos[i] > .1)break; 
+rn += mStPos[i];
+indlim++;
+}
+// float ajout=0.;
+for(i=indlim; i < mMemPos-1; i++)  rn += mStPos[i]/sqrtf(mTePos[i]);
+indlim=0;
+for(i=0; i < mMemNeg-1; i++) 
+{
+if(mTeNeg[i] > .1) break;
+rn += mStNeg[i];
+indlim++;
+}
+for(i=indlim; i < mMemNeg-1; i++) rn += mStNeg[i]/sqrtf(mTeNeg[i]);
+}
+}
+return(rn);
+}
+/* --Methode-- */
+
+int MemNoise::Print()
+{
+int i,j,rc=0;
+printf("MemNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
+printf("mNtirage= %d , mMemPos= %d , mMemNeg= %d \n",mNtirage,mMemPos,mMemNeg);
+// printf("MemNoise::Print() Vecteur de Dynamique / Etat : \n");
+/* for(i=0; i<mMemL-1; i+=2) 
+  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
+         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
+	 */
+return(rc);
+}
+
+
+// Generateur de bruit blanc + 1/f avec fknee
+
+
+/* --Methode-- */
+SumNoise::SumNoise(float fknee, float sig)
+    : NoiseGenerator(sig)
+{
+// Reza 27/01/98 :
+// Je calcule les sigma du MemNoise d'apres la formule de Francois 
+mn = new MemNoise(sig*fknee);  
+}
+
+
+/* --Methode-- */
+SumNoise::~SumNoise()
+{
+delete mn;
+}
+
+/* --Methode-- */
+float SumNoise::Noise()
+{
+return(NoiseGenerator::Noise() + mn->Noise());
+}
Index: trunk/SophyaLib/Samba/bruit.h
===================================================================
--- trunk/SophyaLib/Samba/bruit.h	(revision 228)
+++ trunk/SophyaLib/Samba/bruit.h	(revision 228)
@@ -0,0 +1,108 @@
+#ifndef BRUIT_H_SEEN
+#define BRUIT_H_SEEN
+
+class NoiseGenerator {
+public :
+			NoiseGenerator(float sigma=1.); 
+  virtual		~NoiseGenerator();
+  virtual float		Noise();
+  inline  unsigned long int           NCoups() { return(mNCoups); }
+
+protected :
+unsigned long  int mNCoups;
+  float mSigma;
+};
+
+class  RWalkNoise : public NoiseGenerator {
+public :
+			RWalkNoise(float sigma=1.);
+  virtual		~RWalkNoise();
+  virtual float		Noise();
+
+protected:
+  float mState;
+};
+
+class OOFNoise : public NoiseGenerator {
+public :
+  enum {
+     ACF_Def = 0,
+     ACF_Exp = 1
+  };
+			OOFNoise(float sigma=1., int typacf=ACF_Exp, 
+                                 int mem=128, float tau=1.);
+  virtual		~OOFNoise();
+  virtual float         Noise();
+  void                  Print();
+
+protected:
+  int mTypACF;
+  int mMemL;
+  float mTau; 
+  double* mState;
+  double* mDyn;
+};
+
+class EXPNoise : public NoiseGenerator {
+public :
+  enum {
+     ACF_Def = 0,
+     ACF_Exp = 1
+  };
+			EXPNoise(float sigma=1., int typacf=ACF_Exp, 
+                                 int mem=128, float tau=1.);
+  virtual		~EXPNoise();
+  virtual float         Noise();
+  void                  Print();
+
+protected:
+  int mTypACF;
+long  int mMemL;
+  float mTau; 
+  double* mState;
+  double* mDyn;
+};
+
+class MemNoise : public NoiseGenerator {
+public :
+  enum {
+     ACF_Def = 0,
+     ACF_Exp = 1
+  };
+			MemNoise(float sigma=1., int mem=128, float tau=1., int ava=4096);
+  virtual		~MemNoise();
+  virtual float         Noise();
+  virtual float         Avance(long Asauter);
+  int                   Print();
+ 
+protected:
+  int mTypACF;
+  int mMemL; int mCoupHaut;
+unsigned long   mNappel;
+unsigned long   mNappLast;
+  float mTeappel;
+  float mTau;
+  float mduree;
+  float mTdernier;
+  int mNtirage;
+  int mMemPos;
+  int mMemNeg;
+  float* mStPos;
+  float* mTePos;
+  float* mStNeg;
+  float* mTeNeg;
+};
+
+class  SumNoise : public NoiseGenerator {
+public :
+			SumNoise(float fknee, float sig=1.);
+  virtual		~SumNoise();
+  virtual float		Noise();
+
+protected:
+  MemNoise* mn;
+};
+
+
+#endif
+
Index: trunk/SophyaLib/Samba/gorfph.cc
===================================================================
--- trunk/SophyaLib/Samba/gorfph.cc	(revision 228)
+++ trunk/SophyaLib/Samba/gorfph.cc	(revision 228)
@@ -0,0 +1,39 @@
+#include <iostream.h>
+
+extern "C" {
+void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,
+	      float*,float*,float*);
+void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,
+	      double*, double*,double*,double*,double*,float*);
+void  input_map_(char*,float*,int&);
+void  ecrire_fits_(int&,int&,int&,float*,char*,char*,int&,float&,float&);
+
+}
+
+void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,
+	      float*,float*,float*)
+{
+cerr << "anafast_() bidon ... - Fonction vide - Reza 04/99 " << endl;
+throw " Unsupported function - anafast_() bidon ";
+}
+
+void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,
+	      double*, double*,double*,double*,double*,float*)
+{
+cerr << "synfast_() bidon ... - Fonction vide - Reza 04/99 " << endl;
+throw " Unsupported function - synfast_() bidon ";
+}
+
+void  input_map_(char*,float*,int&)
+{
+cerr << "input_map_() bidon ... - Fonction vide - Reza 04/99 " << endl;
+throw " Unsupported function - input_map_() bidon " ;
+}
+
+void  ecrire_fits_(int&,int&,int&,float*,char*,char*,int&,float&,float&)
+{
+cerr << "ecrire_fits_() bidon ... - Fonction vide - Reza 04/99 " << endl;
+throw " Unsupported function - ecrire_fits_() bidon " ;
+}
+
+
Index: trunk/SophyaLib/Samba/localmap.cc
===================================================================
--- trunk/SophyaLib/Samba/localmap.cc	(revision 228)
+++ trunk/SophyaLib/Samba/localmap.cc	(revision 228)
@@ -0,0 +1,453 @@
+#include "localmap.h"
+#include <iostream.h>
+#include "nbmath.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)
+//
+//    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 
+//    l'equateur, 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(...))
+
+//
+//-- 
+//++
+//
+// Links	Parents
+//
+//    PixelMap
+//    
+//-- 
+//++
+//
+// Links	Descendants
+//
+//     
+//-- 
+//++
+// Titre	Constructeurs
+//--
+//++
+
+LocalMap::LocalMap()
+
+//--
+
+{InitNul();}
+
+//++
+
+LocalMap::LocalMap(int_4 nx, int_4 ny) : mSzX_(nx), mSzY_(ny)
+
+//--
+{
+  InitNul();
+  mNPix_=nx*ny;
+  mPix_ = new r_8[mNPix_];
+}
+
+
+//++
+// Titre	Destructeur
+//--
+//++
+LocalMap::~LocalMap()
+
+//--
+{
+  Clear();
+}
+
+
+
+void        LocalMap::WriteSelf(POutPersist& s) const
+{
+  char strg[256];
+  if (mInfo_) {sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d HasInfo",  mSzX_, mSzY_);} else  {
+    sprintf(strg, "LocalMap: NPixX=%6d  NPixY=%9d ",  mSzX_, mSzY_);
+  }	    
+  s.PutLine(strg);
+  if (mInfo_)  mInfo_->Write(s);
+  s.PutI4(mSzX_);
+  s.PutI4(mSzY_);
+  s.PutI4(mNPix_);
+  s.PutI4(x0_);
+  s.PutI4(y0_);
+  s.PutI4(originFlag_);
+  s.PutI4(SzFlag_);
+  s.PutR4(cos_angle_);
+  s.PutR4(sin_angle_);
+  s.PutR8(theta0_);
+  s.PutR8( phi0_);
+  s.PutR8( tgAngleX_);
+  s.PutR8(tgAngleY_);
+  s.PutR8s(mPix_, mNPix_);
+}
+void        LocalMap::ReadSelf(PInPersist& s)
+{
+  Clear();
+  char strg[256];
+  s.GetLine(strg, 255);
+// Pour savoir s'il y avait un DVList Info associe
+  bool hadinfo = false;
+  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
+  if (hadinfo) {    // Lecture eventuelle du DVList Info
+    if (mInfo_ == NULL)  mInfo_ = new DVList;
+    mInfo_->Read(s);
+  }
+  s.GetI4(mSzX_);
+  s.GetI4(mSzY_);
+  s.GetI4(mNPix_);
+  s.GetI4(x0_);
+  s.GetI4(y0_);
+  s.GetI4(originFlag_);
+  s.GetI4(SzFlag_);
+  s.GetR4(cos_angle_);
+  s.GetR4(sin_angle_);
+  s.GetR8(theta0_);
+  s.GetR8( phi0_);
+  s.GetR8( tgAngleX_);
+  s.GetR8(tgAngleY_);
+  mPix_ = new r_8[mNPix_+1];
+  s.GetR8s(mPix_, mNPix_);
+}
+
+//++
+// Titre	Methodes
+//--
+
+
+//++
+int_4 LocalMap::NbPixels() const
+
+//    Retourne le nombre de pixels du découpage
+//--
+{
+  return(mNPix_);
+}
+
+//++
+r_8& LocalMap::PixVal(int_4 k)
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if ( (k<0) || (k >= mNPix_) ) {
+    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+    //  THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+  THROW(rangeCheckErr);
+
+    //throw "LocalMap::PIxVal Pixel index out of range";
+  }
+  //return (*(pixel_->Data()+k));
+  return(mPix_[k]);
+}
+
+//++
+
+r_8 const& LocalMap::PixVal(int_4 k) const
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if ( (k<0) || (k >= mNPix_) ) {
+    cout << " LocalMap::PIxVal : exceptions  a mettre en place" <<endl;
+  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+
+    throw "LocalMap::PIxVal Pixel index out of range"; 
+  } 
+  //return *(pixel_->Data()+k);
+  return(mPix_[k]);
+}
+
+//++
+int_4 LocalMap::PixIndexSph(float theta, float phi) const
+
+//    Retourne l'indice du pixel vers lequel pointe une direction définie par 
+//    ses coordonnées sphériques
+//--
+{
+  int i,j;
+  if (!(originFlag_==1) || !(SzFlag_==1) ) {
+    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
+    exit(0);
+  }
+  // theta et phi en coordonnees relatives (on se ramene a une situation 
+  // par rapport au plan de reference)
+  float theta_aux=theta;
+  float phi_aux=phi;
+  UserToReference(theta_aux, phi_aux);
+  // coordonnees dans le plan local en unites de pixels
+  float x,y;
+  AngleProjToPix(theta_aux, phi_aux, x,y);
+
+  float xmin=-x0_-0.5;
+  float xmax=xmin+mSzX_;
+  if( (x >  xmax) || (x < xmin ) ) return(-1);
+  float xcurrent=xmin;
+  for( i=0; i < mSzX_; i++ ) {
+    xcurrent += 1.;
+    if( x < xcurrent ) break;
+  }
+  float ymin=-y0_-0.5;
+  float ymax=ymin+mSzY_;
+  if( (y >  ymax) || (y < ymin ) ) return(-1);
+  float ycurrent=ymin;
+  for( j=0; j < mSzY_; j++ ) {
+    ycurrent += 1.;
+    if( y < ycurrent ) break;
+  }
+
+
+  return (j*mSzX_+i);
+}
+
+//++
+void  LocalMap::PixThetaPhi(int_4 k, float& theta, float& phi) const
+
+//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
+//--
+{
+  if (!(originFlag_==1) || !(SzFlag_==1) ) {
+    cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
+    exit(0);
+  }
+  int i,j;
+  Getij(k,i,j);
+  float X=float(i-x0_);
+  float Y=float(j-y0_);
+  // situation de ce pixel dans le plan de reference
+  float x=X*cos_angle_-Y*sin_angle_;
+  float y=X*sin_angle_+Y* cos_angle_;
+  // projection sur la sphere
+  PixProjToAngle(x, y,theta, phi);
+  // retour au plan utilisateur
+  ReferenceToUser(theta, phi);
+}
+ 
+//++
+r_8    LocalMap::PixSolAngle(int_4 k) const
+//    Pixel Solid angle  (steradians)
+//    All the pixels have not necessarly the same size in (theta, phi)
+//    because of the projection scheme which is not yet fixed. 
+//    
+//--
+{
+  int i,j;
+  Getij(k,i,j);
+  float X=float(i-x0_);
+  float Y=float(j-y0_);
+  float XR=X+float(i)*0.5;
+  float XL=X-float(i)*0.5;
+  float YU=Y+float(j)*0.5;
+  float YL=Y-float(j)*0.5;
+  // situation  dans le plan de reference
+  float x0=XL*cos_angle_-YL*sin_angle_;
+  float y0=XL*sin_angle_+YL* cos_angle_;
+  float xa=XR*cos_angle_-YL*sin_angle_;
+  float ya=XR*sin_angle_+YL* cos_angle_;
+  float xb=XL*cos_angle_-YU*sin_angle_;
+  float yb=XL*sin_angle_+YU* cos_angle_;
+  // projection sur la sphere
+  float tet0,phi0,teta,phia,tetb,phib;
+  PixProjToAngle(x0, y0,tet0, phi0);
+  PixProjToAngle(xa, ya,teta, phia);
+  PixProjToAngle(xb, yb,tetb, phib);
+  // angle solide
+  float sol=fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
+  return r_8(sol);
+}
+
+
+//++
+void  LocalMap::SetOrigin(float theta0, float phi0, float angle)
+
+//    
+//--
+{
+  theta0_=theta0;
+  phi0_=phi0;
+  x0_=mSzX_/2;
+  y0_=mSzY_/2;
+  angle=angle*Pi/180.;
+  cos_angle_=cos(angle);
+  sin_angle_=sin(angle);
+  originFlag_ =1; 
+}
+
+ 
+//++
+void  LocalMap::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
+
+//    
+//--
+{
+  theta0_=theta0;
+  phi0_=phi0;
+  x0_=x0;
+  y0_=y0;
+  angle=angle*Pi/180.;
+  cos_angle_=cos(angle);
+  sin_angle_=sin(angle);
+  originFlag_ =1; 
+} 
+
+//++
+void  LocalMap::SetSize(float angleX, float angleY)
+
+//    
+//--
+{
+ tgAngleX_=r_8(tan(angleX*Pi/360.));
+ tgAngleY_=r_8(tan(angleY*Pi/360.));
+ SzFlag_=1;
+}
+//++
+void  LocalMap::Project(SphericalMap& sphere) const
+
+//    Projection to spherical map   
+//--
+{
+  for (int m=0; m<mNPix_; m++) {
+    float theta,phi;
+    PixThetaPhi(m,theta,phi);
+    sphere(theta,phi)=mPix_[m];
+    //    cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
+  }
+}
+
+void LocalMap::InitNul()
+//
+//    initialise à zéro certaines variables internes
+{ 
+originFlag_=0;
+SzFlag_=0;
+cos_angle_=1.;
+sin_angle_=0.;
+mPix_=NULL;
+}
+
+void LocalMap::Clear()
+//
+//    
+{ 
+  if (mPix_)   delete[] mPix_;
+}
+
+void    LocalMap::Getij(int k, int& i, int& j) const
+{
+  i=(k+1)%mSzX_-1;
+  if (i==-1) i=mSzX_-1;
+  j=(k-i+2)/mSzX_;
+}
+
+void  LocalMap::ReferenceToUser(float &theta, float &phi) const
+
+//       
+{
+  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
+    //cout << " LocalMap::ReferenceToUser : exceptions  a mettre en place" <<endl;
+  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+
+    throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 
+  } 
+  //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
+  theta=(r_4)theta0_+theta-(r_4)Pi*0.5;
+  if (theta<0.) {
+    theta=-theta;
+    phi+=(r_4)Pi;
+  } else 
+    if (theta>(r_4)Pi) {
+      theta=2.*(r_4)Pi-theta;
+      phi+=(r_4)Pi;
+    } 
+  
+  phi=(r_4)phi0_+phi;
+  while (phi>=2.*(r_4)Pi) phi-=2.*(r_4)Pi;
+  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
+    cout <<  " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
+    cout << " theta= " << theta << " phi= " << phi << endl;
+    exit(0);
+  }
+  //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
+}
+
+void  LocalMap::UserToReference(float &theta, float &phi) const
+
+//       
+{
+  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
+    cout << " LocalMap::UserToReference : exceptions  a mettre en place" <<endl;
+  //    THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
+
+    throw "LocalMap::UserToReference (theta,phi) out of range"; 
+  } 
+  float phi1=phi-(r_4)phi0_;
+  if (phi1<0.) phi1+=2.*(r_4)Pi;
+
+  float theta1=theta-(r_4)theta0_+(r_4)Pi*0.5;
+  if (theta1<0.) {
+    theta=-theta1;
+    phi1+=(r_4)Pi;
+  } else 
+    if (theta1>(r_4)Pi) {
+      theta=2.*(r_4)Pi-theta1;
+      phi1+=(r_4)Pi;
+    }
+  while (phi1>=2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
+  phi=phi1;
+  if (theta > (r_4)Pi || theta < 0.  || phi<0. || phi >=2* (r_4)Pi ) {
+    cout <<  " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
+    cout << " theta= " << theta << " phi= " << phi << endl;
+    exit(0);
+  }
+}
+
+void LocalMap::PixProjToAngle(float x, float y,float &theta, float &phi) const {
+
+// (x,y) representent les coordonnees en unites de pixels d'un point DANS 
+// LE PLAN DE REFERENCE. On recupere (theta,phi) par rapport au repere "absolu"
+// theta=pi/2,phi=0.
+
+theta=(r_4)Pi*0.5-atan(2*y* tgAngleY_/(float) mSzY_);
+phi=atan2(2*x* tgAngleX_,(float) mSzX_);
+if( phi< 0. ) phi+= DeuxPi;
+}
+
+void  LocalMap::AngleProjToPix(float theta, float phi, float& x, float& y) const {
+
+// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere 
+// (i,j)  DANS LE PLAN DE REFERENCE. 
+  if (phi>=(r_4)Pi) phi-= DeuxPi;
+  //  y=0.5*mSzY_*cot(theta)/tgAngleY_;  $CHECK-REZA-04/99$
+  y=0.5*mSzY_/tan(theta)/tgAngleY_;  // ? cot = 1/tan ?  
+  x=0.5*mSzX_*tan(phi)/tgAngleX_;
+}
+
+
+
+
+
Index: trunk/SophyaLib/Samba/localmap.h
===================================================================
--- trunk/SophyaLib/Samba/localmap.h	(revision 228)
+++ trunk/SophyaLib/Samba/localmap.h	(revision 228)
@@ -0,0 +1,123 @@
+#ifndef LOCALMAP_SEEN
+#define LOCALMAP_SEEN
+
+#include "pixelmap.h"
+#include "sphericalmap.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)
+//
+//    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 
+//    l'equateur, 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 : public PixelMap {
+public:
+   LocalMap();
+   LocalMap(int_4 nx, int_4 ny);
+   virtual ~LocalMap();
+   // Overloading of () to access pixel number k.
+inline  r_8&        operator()(int_4 k)
+                          {return(PixVal(k));}
+inline  r_8 const&  operator()(int_4 k) const
+                          {return(PixVal(k));}
+inline  r_8&        operator()(int ix, int iy)
+                { return PixVal(iy*mSzX_+ix) ; };
+inline  r_8 const&        operator()(int ix, int iy) const
+                { return PixVal(iy*mSzX_+ix) ; };
+   
+// ---------- Persistence handling
+   
+   enum {classId = 0xF003 };
+int_4               ClassId() const        { return classId; }
+
+virtual void        WriteSelf(POutPersist&) const;
+virtual void        ReadSelf(PInPersist&);
+
+// ---------- Definition of PixelMap abstract methods
+
+   // Number of pixels
+   virtual int_4       NbPixels() const;
+   
+   // Value of pixel number k
+   virtual r_8&        PixVal(int_4 k);
+   virtual r_8 const&  PixVal(int_4 k) const;
+
+   // Index of pixel at (theta,phi)
+   virtual int_4       PixIndexSph(float theta, float phi) const;
+   
+   // Spherical coordinates of center of pixel number k
+   virtual void        PixThetaPhi(int_4 k, float& theta, float& phi) const;
+   // Pixel Solid angle  (steradians)
+   virtual r_8         PixSolAngle(int_4 k) const; 
+
+// ---------- Specific methods
+
+   // Origin (with angle between x axis and phi axis, in degrees)
+    virtual void  SetOrigin(float theta0, float phi0, float angle=0.); // x0,y0 default: middle of map
+  virtual void        SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle=0.);
+   // Pixel size (degres) 
+   virtual void        SetSize(float angleX, float angleY);
+
+   // Projection to/from spherical map
+   //virtual void        Extract(SphericalMap const& sphere);
+   virtual void        Project(SphericalMap& sphere) const;
+  
+   // There should be a more complex algorithm somewhere to combine *several*
+   // local maps to a full sphere.
+   //   -> static method, or separate class
+
+   void Pixelize(int_4,int_4); // Allocate pixel array
+
+// ------------- méthodes internes ----------------------
+          
+private :
+
+   void	     InitNul();
+   void      Clear();
+   void      Getij(int k, int& i, int& j) const;
+   void      ReferenceToUser(float &theta, float &phi) const; 
+   void      UserToReference(float &theta, float &phi) const;
+   void      PixProjToAngle(float x, float y,float &theta, float &phi) const; 
+   void      AngleProjToPix(float theta, float phi, float& x, float& y) const; 
+// ------------- variables internes -----------------------
+
+   int_4   mSzX_, mSzY_;
+   int_4   mNPix_;
+   int_4   x0_;
+   int_4   y0_;
+   int_4    originFlag_;
+   int_4    SzFlag_;
+   r_4   cos_angle_;
+   r_4   sin_angle_;
+   r_8     theta0_;
+   r_8     phi0_;
+   r_8     omeg_;
+   r_8     tgAngleX_;
+   r_8     tgAngleY_;
+   r_8*    mPix_;
+};
+
+#endif
Index: trunk/SophyaLib/Samba/mlobe.cc
===================================================================
--- trunk/SophyaLib/Samba/mlobe.cc	(revision 228)
+++ trunk/SophyaLib/Samba/mlobe.cc	(revision 228)
@@ -0,0 +1,172 @@
+#include "mlobe.h"
+#include "unitvec.h"
+
+#include "timing.h"
+// #include "nbrandom.h"
+
+/* --Methode-- */
+MainLobe::MainLobe(float sig, float dels, int nc)
+{
+if (sig < 1.e-9)  sig = M_PI/100.;   // Sigma du lobe gaussien (radians)
+mSigma = sig;
+mDels = (dels < 1.e-5) ?  0.5 : dels;   // Epaisseur de chaque couche en unite de sigma
+mNC =  (nc < 1) ? 1 : nc;
+mNpix = 1;    // Nb total de pixels du lobe
+
+int i,j,k;
+for(i=1; i<=mNC; i++)  mNpix += 6*i;
+//  mT0 , mP0  Teta et Phi des pixels, Pixel 0 oriente vers (0, 0)
+mT0 = new float[mNpix];
+mP0 = new float[mNpix];
+//  mTC , mPC  Teta et Phi des pixels, Pixel 0 oriente vers (Teta, phi)
+mTC = new float[mNpix];
+mPC = new float[mNpix];
+mDx = new float[mNpix];
+mDy = new float[mNpix];
+mDr = new float[mNpix];
+// Valeurs des pixels 
+mVal = new double[mNpix];
+
+
+double teta,phi,dphi;
+double val;
+mT0[0] = mP0[0] = mTC[0] = mPC[0] = 0.;
+mVal[0] = 1.;  mDx[0] = mDy[0] = 0.;
+k = 1;  mTotVal = 1.;
+for(i=1; i<=mNC; i++) {
+  teta = (double)i*mDels;
+  val = exp(-teta*teta);
+  teta *= mSigma;
+  if (i%2)  phi = M_PI/6.;
+  else phi = 0.;
+  dphi = M_PI/(3.*i);
+  for(j=0; j<6*i; j++) {
+    mT0[k] = mTC[k] = teta;
+    phi += dphi;
+    mP0[k] = mPC[k] = phi;
+    mVal[k] = val;
+    mTotVal += val;
+    mDx[k] = teta*cos(phi);
+    mDy[k] = teta*sin(phi);
+    mDr[k] = teta;
+    k++;
+    }
+  }
+mDirTet = mDirPhi = 0.;
+return;
+}
+
+
+/* --Methode-- */
+MainLobe::~MainLobe()
+{
+delete mT0;
+delete mP0;
+delete mTC;
+delete mPC;
+delete mDx;
+delete mDy;
+delete mDr;
+delete mVal;
+}
+
+/* --Methode-- */
+void MainLobe::SetDirection(float teta, float phi)
+{
+int k;
+
+mDirTet = teta;
+mDirPhi = phi;
+
+if (teta < 0.1*mDels*mSigma) {
+  for(k=0; k<mNpix; k++) {
+    mTC[k] = mT0[k]; mPC[k] = mP0[k]; 
+    }
+  return;
+  }
+
+
+/*
+if (teta > 10.*mSigma*mDels*mNC) {   // On utilise des formules approchees
+  float stet = sin((double)teta);
+  float ctet = cos((double)teta);
+  for(k=0; k<mNpix; k++) {
+    mTC[k] = teta - mDy[k]*stet - mDr[k]*ctet;
+    mPC[k] = phi + mDx[k];
+    }
+  return;
+  }
+  */
+
+// Formules compliquees 
+
+double dx,dy,dz;
+double sphi = sin((double)phi);
+double cphi = cos((double)phi);
+double stet = sin((double)teta);
+double ctet = cos((double)teta);
+
+UnitVector pd((double)teta, (double)phi);
+UnitVector vc;
+for(k=0; k<mNpix; k++) {
+  dz = mDy[k]*stet;
+  dx = -(mDx[k]*sphi+mDy[k]*ctet*cphi);
+  dy = mDx[k]*cphi-mDy[k]*ctet*sphi;
+  vc.setxyz(pd.x+dx, pd.y+dy, pd.z+dz);
+  mTC[k] = vc.Teta();
+  mPC[k] = vc.Phi();
+  if (mTC[k] < 0.)  printf("ERREUR !!!! mTC[%d] = %g \n", k, mTC[k]);
+  if (mPC[k] < 0.) mPC[k] += M_PI*2.;
+  if (mPC[k] >= M_PI*2.) mPC[k] -= M_PI*2.;
+  }
+return;  
+}
+
+/* --Methode-- */
+double MainLobe::Value(int kpx, float& teta, float& phi)
+{
+if ( (kpx < 0) || (kpx >= mNpix) ) { teta = phi = 0.; return(0.); }
+teta = mTC[kpx];
+phi = mPC[kpx];
+return(mVal[kpx]);
+}
+
+/* --Methode-- */
+double MainLobe::Convol(SphericalMap& sph)
+{
+double ret=0.;
+int k;
+for(k=0; k<mNpix; k++)  ret += mVal[k]*sph(mTC[k], mPC[k]);
+/*
+if( (dbg > 0) && (k%dbg) ) 
+  printf("MainLobe::Convol k= %d T,P= %g %g , Lobe=%g Sph= %g  Ret=%g\n", 
+         k, mTC[k], mPC[k], mVal[k], sph(mTC[k], mPC[k]), ret );
+*/
+return(ret/mTotVal);
+}
+
+/* --Methode-- */
+ostream& operator << (ostream& s, MainLobe const& lob)
+{
+s << "MainLobe_Info() : Sigma= " << lob.mSigma << " Dels= " << lob.mDels << " NC= " 
+  << lob.mNC << " NTotPix= " << lob.mNpix << " Integ= " << lob.mTotVal << endl;
+
+int i,j,k;
+cout << " Direction Teta,Phi = " << lob.mDirTet << " , " << lob.mDirPhi << endl;      
+cout << " Pixel 0 , Teta= " << lob.mT0[0] << " Phi= " << lob.mP0[0] 
+     << " Val= " << lob.mVal[0] << endl;
+k = 1;
+for(i=1; i<lob.mNC; i++) {
+  cout << " Couche # " << i << " NPix= " << i*6 << endl;
+  for(j=0;j<6*i;j++) {
+    cout << "K= " << k << "Teta,Phi= "  << lob.mTC[k] << ", " << lob.mPC[k] 
+         << " (" << lob.mT0[k] << ", " << lob.mP0[k] << ") "
+         << " Val= " << lob.mVal[k] << endl; 
+    k++;
+    }
+  }
+
+
+return s;
+}
+
Index: trunk/SophyaLib/Samba/mlobe.h
===================================================================
--- trunk/SophyaLib/Samba/mlobe.h	(revision 228)
+++ trunk/SophyaLib/Samba/mlobe.h	(revision 228)
@@ -0,0 +1,50 @@
+#ifndef SCSKYU_H_SEEN
+#define SCSKYU_H_SEEN
+
+#include <math.h>
+#include <iostream.h>
+
+#include "sphericalmap.h"
+
+
+
+// Classe MainLobe
+// Classe pour calcul du lobe principal
+// Les valeurs du lobe sont calcules dans des pixels
+// disposes sur "nc" couches avec pavage hexagonal autour de Teta=0, Phi=0
+// Le lobe est gaussien, de sigma = "sig" - Les couches sont construites 
+// toutes les "dels" * sigma 
+// NPix() renvoie le nombre total de pixels du lobe 
+// SetDirection() permet d orienter le lobe dans la direction "teta, phi"
+// Value() renvoie la valeur du lobe, ainsi que les teta,phi du pixel "kpx"
+// Convol() renvoie la valeur lobe convolue avec le ciel (sphere) "sph"
+
+class MainLobe {
+public:
+			MainLobe(float sig, float dels=0.5, int nc=3);
+			~MainLobe();
+        inline int      NPix() {return mNpix; }            
+        inline float    Sigma() {return mSigma; }            
+        inline double   Integral() {return mTotVal; }            
+               void     SetDirection(float teta, float phi);
+	       double   Value(int kpx, float& teta, float& phi);
+               double   Convol(SphericalMap& sph);
+        friend ostream& operator << (ostream& s, MainLobe const& lob);
+	inline void     Print() { cout << (*this); } 
+
+float mSigma, mDels;
+int mNpix, mNC;
+double mTotVal;
+float mDirTet, mDirPhi;
+float* mT0; 
+float* mP0; 
+float* mTC; 
+float* mPC; 
+float* mDx;
+float* mDy;
+float* mDr;
+double* mVal;
+};
+
+
+#endif
Index: trunk/SophyaLib/Samba/pdataarray.h
===================================================================
--- trunk/SophyaLib/Samba/pdataarray.h	(revision 228)
+++ trunk/SophyaLib/Samba/pdataarray.h	(revision 228)
@@ -0,0 +1,46 @@
+#ifndef PDATAARRAY_H
+#define PDATAARRAY_H
+#include "prefcount.h"
+
+template <class T>
+class PDataArray : public PRefCounted {
+public:
+                            PDataArray(long n, 
+                                       bool isTmp=false)
+                                                    :PRefCounted(isTmp) 
+                                                    {data = new T[n]; sz=n;}
+                            ~PDataArray()           {delete[] data;}
+  
+  T*                        Data()                  {return data;}
+  T const*                  Data() const            {return data;}
+  long                      Size() const            {return sz;}
+  
+  // Pointeurs/iterateurs
+  T*                        begin()                 {return data;}
+  T const*                  begin() const           {return data;}
+  T*                        end()                   {return data+sz;}
+  T const*                  end() const             {return data+sz;}
+
+  
+  PDataArray<T>*            Clone(bool isTmp=false) const;
+  
+private:
+  T*                        data;
+  int                       sz;
+};
+
+
+template <class T>
+PDataArray<T>* PDataArray<T>::Clone(bool isTmp) const
+{
+#ifdef TRACE_COPY
+  cerr << "PDataArray::Clone" << endl;
+#endif
+   PDataArray<T>* p = new PDataArray<T>(sz, isTmp);
+   memcpy(p->data, data, sz*sizeof(T));
+   return p;
+}
+
+
+#endif
+
Index: trunk/SophyaLib/Samba/pixelmap.h
===================================================================
--- trunk/SophyaLib/Samba/pixelmap.h	(revision 228)
+++ trunk/SophyaLib/Samba/pixelmap.h	(revision 228)
@@ -0,0 +1,70 @@
+#ifndef PIXELMAP_SEEN
+#define PIXELMAP_SEEN
+
+#include "ppersist.h" 
+#include "dvlist.h"
+#include <iostream.h>
+
+// General map of pixels on part of sphere or whole sphere
+// Class hierarchy :
+//   PixelMap
+//      SphericalMap
+//         SphereThetaPhi
+//         SphereGorski
+//         SphereIco
+//      LocalMap
+
+
+class PixelMap : public PPersist {
+public:
+   PixelMap():mInfo_(NULL) {}
+   virtual ~PixelMap() {}
+   
+   // Number of pixels
+   virtual int_4       NbPixels() const=0;
+   
+   // Value of pixel number k
+   virtual r_8&        PixVal(int_4 k)=0;
+   virtual r_8 const&  PixVal(int_4 k) const=0;
+   
+   // Index of pixel at (theta,phi)
+   virtual int_4       PixIndexSph(float theta, float phi) const=0;
+
+   // Value of pixel number at (theta,phi)
+   virtual r_8&        PixValSph(float theta, float phi)
+                           {return PixVal(PixIndexSph(theta,phi));}
+   virtual r_8 const&  PixValSph(float theta, float phi) const
+                           {return PixVal(PixIndexSph(theta,phi));}                           
+   
+   // Spherical coordinates of center of pixel number k
+   virtual void        PixThetaPhi(int_4 k, float& theta, float& phi) const=0;
+
+   // Pixel  (steradians)
+   virtual r_8         PixSolAngle(int_4 k) const =0; 
+
+   // Overloading of () to access pixel number k.
+   inline  r_8&        operator()(int_4 k)
+                          {return(PixVal(k));}
+   inline  r_8 const&  operator()(int_4 k) const
+                          {return(PixVal(k));}
+
+   // Note : no overloading of (float,float) to access pixel (theta,phi).
+   // overloading of (float,float) in SphericalMap
+   // overloading of (int,int)     in CartesianMap
+
+//++
+DVList&  Info()
+//
+//      Renvoie une rM-ifM-irence sur l'objet DVList AssociM-i
+//--
+{
+if (mInfo_ == NULL)  mInfo_ = new DVList;
+return(*mInfo_);
+}
+
+   protected :
+
+  DVList* mInfo_;        // Infos (variables) attachees 
+};
+
+#endif
Index: trunk/SophyaLib/Samba/prefcount.h
===================================================================
--- trunk/SophyaLib/Samba/prefcount.h	(revision 228)
+++ trunk/SophyaLib/Samba/prefcount.h	(revision 228)
@@ -0,0 +1,65 @@
+#ifndef PREFCOUNT_H
+#define PREFCOUNT_H
+
+// Classe d'objet avec comptage de référence.
+// Le nombre de référence est incrémenté lors de la création et lors d'un Attach
+// (inutile d'appeler Attach si on crée l'objet).
+// Si on cree l'objet avec isTmp = true, alors le possesseur est temporaire,
+// et en tout cas, possede le PRefCount en lecture seule.
+// Par exemple, une fonction qui retourne un objet de type PRefCounted le créera
+// en mode temporaire. Un constructeur par copie appelé derrière ne dupliquera pas
+// les données.
+// HasTmpOnly retourne true si un objet n'est référencé que de façon temporaire.
+
+#ifdef TRACE_REFCOUNT
+#define DBG_DESTROY cerr << "PRefCounted destroyed" << endl;
+#define DBG_ATTACH  cerr << "PRefCounted attached" << endl;
+#define DBG_DETACH  cerr << "PRefCounted detached" << endl;
+#else
+#define DBG_DESTROY 
+#define DBG_ATTACH
+#define DBG_DETACH
+#endif
+
+
+class PRefCounted {
+public:
+                               PRefCounted(bool isTmp=false) 
+                                                        :refcount(1), tmpcount(isTmp?1:0) {}
+  virtual                      ~PRefCounted()           {DBG_DESTROY}
+  
+  void                         Attach(bool isTemp=false) 
+                                                        {if (isTemp && HasTmpOnly()) tmpcount++;
+                                                         DBG_ATTACH
+                                                         refcount++;}
+  void                         Detach()                 {refcount--;
+                                                         DBG_DETACH
+                                                         if (!refcount) delete this;}
+  
+  int                          RefCount() const         {return refcount;}
+  bool                         HasTmpOnly() const       {return tmpcount == refcount;}
+  
+private:
+  int refcount;
+  int tmpcount;
+};
+
+// Macro pour les constructeurs par copie d'objets référençant des PRefCounted
+// Il faut que le PRefCounted ait une méthode T* Clone(bool isTmp);
+
+template <class refcount>
+void PRefCountCopy(refcount*& dest, refcount const* src, bool isTmp=false) {
+  if (dest) dest->Detach();
+  if (!src) {
+    dest=0;
+  } else if (src->HasTmpOnly()) {
+    dest = const_cast<refcount*>(src);
+    dest->Attach(isTmp);
+  } else {
+    dest = src->Clone(isTmp);
+  }
+}
+
+
+#endif
+
Index: trunk/SophyaLib/Samba/scan.cc
===================================================================
--- trunk/SophyaLib/Samba/scan.cc	(revision 228)
+++ trunk/SophyaLib/Samba/scan.cc	(revision 228)
@@ -0,0 +1,331 @@
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "scan.h"
+
+// valeurs de Pi, 2*Pi, etc                                              
+#include "nbmath.h"
+
+
+//
+//++
+// Class	Scan
+//
+// include      scan.h  nbmath.h 
+//
+//    Cette classe permet de stocker et traiter les données pour un balayage 
+//    d'une portion de ciel pour un jeu de valeurs fixé des paramètres précisé
+//    dans le constructeur.    
+//--
+//++
+//
+// Links	Parents
+//
+//    
+//-- 
+//++
+//
+// Links	Descendants
+//
+//     
+//-- 
+//++
+// Titre	Constructeurs
+//--
+//++
+Scan::Scan(float Ouv,float* Omega,float Fech,float T, 
+	   float t0=0., float phi0=0.)
+
+//  
+//    * Ouv      = angle d'ouverture (rad)
+//    * Omega[3] = direction de l'axe de rotation du satellite (teta,phi)
+//    et vitesse de rotation (rad/s)
+//    * Fech     = fréquence d'échantillonnage (Hz)
+//    * T        = temps total de prise des données (s)
+//    * t0       = instant de départ (s)
+//    * phi0     = offset antenne (rad)
+//--
+{
+  mInfo_=NULL;
+  // On memorise les arguments d'appel
+  Ouverture_= Ouv;
+  OmegaTeta_= Omega[0];
+  //
+  // je ne comprends pas ce qui se passe ci-apres (GLM)
+  if( (Omega[0]== 0.0) || (Omega[0]==(float)Pi) ) 
+    OmegaPhi_ = Omega[1];
+  else
+    OmegaPhi_ = Omega[1]+(float)Pi/2.;
+  OmegaRad_ = Omega[2];
+  FrequenceEch_  = Fech;
+  TempsFinal_    = T;
+  TempsInitial_  = t0;
+  PhiZero_ = phi0;
+
+  float aux= 0.0;
+  // Nombre total de points
+  aux= (TempsFinal_-TempsInitial_)*FrequenceEch_;
+  NmaxPts_= (int_4)(aux);
+
+  // Nombre total de tours
+  aux= OmegaRad_*(TempsFinal_-TempsInitial_)/(2.*Pi);
+  NmaxTrs_= (int_4)(aux);
+
+  // Nombre de points par tour
+  aux= 2.*Pi*FrequenceEch_/OmegaRad_;
+  NPts1Tr_= (int_4)(aux);
+
+  // Creation et initialisation du vecteur des mesures aux points
+  sPix_ = new r_8[NmaxPts_];
+  for( int_4 i=0; i<NmaxPts_; i++ )  sPix_[i]= 0.;
+
+  // matrice de passage du systeme lie au satellite au systeme fixe
+
+  Rota_[0]=  cos((double)OmegaPhi_);
+  Rota_[1]= -sin((double)OmegaPhi_)*cos((double)OmegaTeta_);
+  Rota_[2]=  sin((double)OmegaPhi_)*sin((double)OmegaTeta_);
+  Rota_[3]=  sin((double)OmegaPhi_);
+  Rota_[4]=  cos((double)OmegaPhi_)*cos((double)OmegaTeta_);
+  Rota_[5]= -cos((double)OmegaPhi_)*sin((double)OmegaTeta_);
+  Rota_[6]=  0.0;
+  Rota_[7]=  sin((double)OmegaTeta_);
+  Rota_[8]=  cos((double)OmegaTeta_);
+
+
+  // printf("%f  %f  %f \n",Rota_[0],Rota_[1],Rota_[2]);
+  // printf("%f  %f  %f \n",Rota_[3],Rota_[4],Rota_[5]);
+  // printf("%f  %f  %f \n",Rota_[6],Rota_[7],Rota_[8]);
+
+}
+//++
+Scan::Scan(const Scan& s)
+
+//    Constructeur de copie
+//--
+{
+  NmaxPts_      = s. NmaxPts_;
+  Ouverture_    = s.Ouverture_;
+  OmegaTeta_    = s.OmegaTeta_;
+  OmegaPhi_     = s.OmegaPhi_;
+  OmegaRad_     = s.OmegaRad_;
+  FrequenceEch_ = s. FrequenceEch_;
+  TempsFinal_ = s.TempsFinal_;
+  TempsInitial_ = s.TempsInitial_;
+  PhiZero_ = s. PhiZero_;
+  sPix_=new r_8[ NmaxPts_];
+  for (int_4 k=0; k<NmaxPts_; k++)  sPix_[k]=s.sPix_[k];
+  for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k];
+}
+//++
+// Titre	Destructeur
+//--
+//++
+Scan::~Scan()
+
+//     
+//--
+{
+  delete[] sPix_ ;
+}
+//++
+void Scan::WriteSelf(POutPersist& s) const
+
+//    créer un fichier persistant
+//--
+{
+  char strg[256];
+  if (mInfo_) {sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f HasInfo", (float)OmegaTeta_, (float)OmegaPhi_, (float)OmegaRad_);
+  }
+  else {sprintf(strg, "Scan: Theta=%9f Phi=%9f omega=%9f ", (float)OmegaTeta_, (float)OmegaPhi_, (float)OmegaRad_); 
+  }
+  s.PutLine(strg);
+  if (mInfo_)  mInfo_->Write(s);
+   s.PutI4( NmaxPts_);
+  s.PutI4(NmaxTrs_);
+  s.PutI4( NPts1Tr_);
+  s.PutR4(Ouverture_);
+  s.PutR4(OmegaTeta_);
+  s.PutR4(OmegaPhi_);
+  s.PutR4(OmegaRad_);
+  s.PutR4(FrequenceEch_);
+  s.PutR4(TempsFinal_);
+  s.PutR4( TempsInitial_);
+  s.PutR4( PhiZero_);
+  s.PutR8s(sPix_,NmaxPts_);
+  s.PutR8s(Rota_, 9);
+  return;
+}
+
+
+//++
+void Scan::ReadSelf(PInPersist& s) 
+
+//    relire un fichier persistant
+//--
+{
+  Clear();
+  char strg[256];
+  s.GetLine(strg,255);
+// Pour savoir s'il y avait un DVList Info associe
+  bool hadinfo = false;
+  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
+// Pour savoir s'il y avait un DVList Info associe
+  if (hadinfo) {    // Lecture eventuelle du DVList Info
+    if (mInfo_ == NULL)  mInfo_ = new DVList;
+  mInfo_->Read(s);
+  }
+  s.GetI4( NmaxPts_);
+  s.GetI4(NmaxTrs_);
+  s.GetI4( NPts1Tr_);
+  s.GetR4(Ouverture_);
+  s.GetR4(OmegaTeta_);
+  s.GetR4(OmegaPhi_);
+  s.GetR4(OmegaRad_);
+  s.GetR4(FrequenceEch_);
+  s.GetR4(TempsFinal_);
+  s.GetR4( TempsInitial_);
+  s.GetR4( PhiZero_);
+  sPix_=new r_8[NmaxPts_];
+  s.GetR8s(sPix_,NmaxPts_);
+  s.GetR8s(Rota_, 9);
+}
+
+
+//++
+// Titre	Méthodes
+//--
+
+//++
+int_4 Scan::NbPoints() const
+
+//    Retourne le nombre de Points du scan
+//--
+{
+  return(NmaxPts_);
+}
+
+/* --Methode-- */
+//++
+int_4 Scan::NbTours() const
+
+//    Retourne le nombre total de tours
+//--
+{
+  return(NmaxTrs_);
+}
+
+/* --Methode-- */
+//++
+int_4 Scan::NbPts1Tr() const
+
+//    Retourne le nombre de points pour 1 tour
+//--
+{
+  return(NPts1Tr_);
+}
+
+/* --Methode-- */
+//++
+int_4 Scan::ValueIndex(float t) const
+
+//    Retourne l'indice du pixel associé au temps t
+//--
+{
+  int_4 k;
+  float eps= 1.0E-06;
+
+  // verification si t est dans [TempsInitial_,TempsFinal_]
+  if( (t< TempsInitial_) || (t> TempsFinal_) ) {
+    printf("\n ValueIndex:: probleme sur le temps t= %f",t);
+    return(-1);
+  }
+
+  k= (int_4)((t-TempsInitial_)*FrequenceEch_+eps);
+  if ( (k< 0) || (k >= NmaxPts_) ) {
+    printf("\n ValueIndex:: probleme sur l'indice du point k= %d",k); 
+    return(-1);
+  }
+  return(k);
+}
+
+/* --Methode-- */
+//++
+void Scan::Direction(float t, float&  teta , float&  phi)
+
+//    Retourne les coordonnées (teta,phi) du pixel lié au temps t
+//--
+{
+  r_8 alfa,xs,ys,zs,x,y,z;
+
+  // coordonnees dans le systeme du satellite
+  alfa= OmegaRad_*(t-TempsInitial_)+PhiZero_;
+  xs = sin((double)Ouverture_)*cos(alfa);
+  ys = sin((double)Ouverture_)*sin(alfa);
+  zs = cos((double)Ouverture_);
+
+  // coordonnees dans le systeme fixe
+  x = Rota_[0]*xs+Rota_[1]*ys+Rota_[2]*zs;
+  y = Rota_[3]*xs+Rota_[4]*ys+Rota_[5]*zs;
+  z = Rota_[6]*xs+Rota_[7]*ys+Rota_[8]*zs;
+
+  // angles teta,phi
+  teta = acos(z);
+  phi  = atan2(y,x);
+  if( phi< 0. ) phi= DeuxPi+phi;
+}
+
+/* --Methode-- */
+//++
+void Scan::DirectionIndex(int_4 k, float&  teta, float&  phi)
+
+//    Retourne les coordonnées (teta,phi) du pixel d'indice k
+//--
+{
+  float t;
+
+  //recupere le temps associe a l'indice k du pixel
+  t= TempsInitial_+(float)k/FrequenceEch_;
+
+  // angles teta,phi 
+  Direction(t, teta, phi);
+}
+
+static r_8 dummy_pixel = 0;
+/* --Methode-- */
+//++
+r_8& Scan::PixelValue(int_4 k) const
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if ( (k<0) || (k >=  NmaxPts_) ) {
+    printf("\n ValueIndex::indice du pixel errone k= %d",k); 
+    return(dummy_pixel);
+  }
+  return(sPix_[k]);
+}
+
+void Scan::Clear() {
+  if (sPix_)  delete[] sPix_;
+}         
+void Scan::InitNull() {
+  sPix_=NULL;
+  mInfo_=NULL;
+}         
+
+
+//++
+// Titre	Opérateurs
+//--
+//++
+// Links	double&        operator()(int_4 k)
+//--
+//++
+//    Surcharge de la parenthese a un indice entier : remplit ou/et renvoie 
+//    la valeur du contenu du pixel d'indice k                               
+//--
+
+
+
+
+
Index: trunk/SophyaLib/Samba/scan.h
===================================================================
--- trunk/SophyaLib/Samba/scan.h	(revision 228)
+++ trunk/SophyaLib/Samba/scan.h	(revision 228)
@@ -0,0 +1,88 @@
+#ifndef SCAN_H_SEEN
+#define SCAN_H_SEEN
+
+
+// includes PEIDA
+#include "ppersist.h"
+#include <iostream.h>
+#include "dvlist.h"
+class Scan : public PPersist {
+
+  public  :
+    /* 
+arguments du constructeur : 
+  . ouverture (radians), 
+  . un tableau de 3 reels (deux angles-en radians- definissant la direction de 
+                      l'axe de rotation du satellite + vitesse de rotation 
+                      en radians/s),
+  . frequence d'echantillonnage (Hz), 
+  . instant final de prise de donnees (s),
+deux arguments optionnels :
+  . instant initial pour la prise de donnees
+  . offset d'antenne en angle phi
+  */
+Scan(float,float*,float,float,float,float);
+Scan(const Scan& );
+Scan() { 
+  InitNull();
+  cout << " constructeur de scan par defaut" << endl;
+}
+~Scan();
+
+// ------------ Persistence handling
+
+     enum {classId = 0xF003 };
+int_4                  ClassId() const        { return classId; }
+
+virtual void           WriteSelf(POutPersist&) const;
+virtual void           ReadSelf(PInPersist&);
+
+int_4     NbPoints() const;
+int_4     NbTours() const;
+int_4     NbPts1Tr() const; 
+int_4     ValueIndex(float) const;
+void    Direction(float, float& ,float& );
+void    DirectionIndex(int_4,float& ,float& );
+double& PixelValue(int_4 k) const;
+
+/* Surcharge de la parenthese a un indice entier : remplit ou/et renvoie */
+/* la valeur du contenu du pixel d'indice k                              */ 
+
+  inline  double&        operator()(int_4 k)
+                                      { return(PixelValue(k)) ; } ;
+//++
+DVList&  Info()
+//
+//      Renvoie une reference sur l'objet DVList Associe
+//--
+{
+if (mInfo_ == NULL)  mInfo_ = new DVList;
+return(*mInfo_);
+}
+
+
+  private :
+    
+void          Clear();
+void          InitNull();
+
+int_4 NmaxPts_;
+int_4 NmaxTrs_;
+int_4 NPts1Tr_;
+r_4   Ouverture_;
+r_4   OmegaTeta_;     /* direction de l'axe de rotation du satellite */
+r_4   OmegaPhi_;      /* (radians)                                   */
+r_4   OmegaRad_;      /* vitesse de rotation du satellite (radians/s) */
+r_4   FrequenceEch_;  /* frequence d'echantillonnage (Hz)*/
+r_4   TempsFinal_;
+r_4   TempsInitial_;
+r_4   PhiZero_;       /* Offset antenne (radians)  */
+r_8*  sPix_;          /* valeurs contenues dans chaque pixel */ 
+r_8   Rota_[9];       /* matrice de passage : [Xf]=[Rota][Xs] */
+                      /* Xs : coordonnees dans un systeme lie au satellite*/ 
+                      /* Xf : coordonnees dans un repere "fixe"           */ 
+DVList* mInfo_;        // Infos (variables) attachees 
+};
+
+
+#endif /* SCAN_H_SEEN */
Index: trunk/SophyaLib/Samba/spheregorski.cc
===================================================================
--- trunk/SophyaLib/Samba/spheregorski.cc	(revision 228)
+++ trunk/SophyaLib/Samba/spheregorski.cc	(revision 228)
@@ -0,0 +1,1360 @@
+//
+#include "spheregorski.h"
+#include "strutil.h"
+extern "C" {
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+}
+      
+extern "C" {
+  //void ang2pix_ring_(int&,double&,double&,int&);
+  //void pix2ang_ring_(int&, int&, double&, double&);
+//void ang2pix_nest_(int&,double&,double&,int&);
+  //void pix2ang_nest_(int&, int&, double&, double&);
+//void nest2ring_(int&, int&, int&);
+//void ring2nest_(int&, int&, int&);
+void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,
+	      float*,float*,float*);
+void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,
+	      double*, double*,double*,double*,double*,float*);
+void  input_map_(char*,float*,int&);
+void  ecrire_fits_(int&,int&,int&,float*,char*,char*,int&,float&,float&);
+
+}
+//*******************************************************************
+//++ 
+// Class	SphereGorski
+//
+// include      spheregorski.h strutil.h
+//
+//    Pixelisation Gorski  
+//
+//
+//|    -----------------------------------------------------------------------
+//|     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
+//|    -----------------------------------------------------------------------
+//
+//    the sphere is split in 12 diamond-faces containing nside**2 pixels each
+//
+//    the numbering of the pixels (in the nested scheme) is similar to
+//    quad-cube
+//    In each face the first pixel is in the lowest corner of the diamond
+//
+//    the faces are                    (x,y) coordinate on each face
+//|        .   .   .   .   <--- North Pole
+//|       / \ / \ / \ / \                          ^        ^     
+//|      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      /      
+//|       \ / \ / \ / \ /                        y   \    /  x   
+//|      4 . 5 . 6 . 7 . 4 <--- equator               \  /                     
+//|       / \ / \ / \ / \                              \/        
+//|      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner  
+//|       \ / \ / \ / \ /                                  
+//|        .   .   .   .   <--- South Pole
+//|
+//    phi:0               2Pi                                 
+//
+//    in the ring scheme pixels are numbered along the parallels
+//    the first parallel is the one closest to the north pole and so on
+//    on each parallel, pixels are numbered starting from the one closest
+//    to phi = 0
+//
+//    nside DOIT OBLIGATOIREMENT ETRE UNE PUISSANCE DE 2 (<= 8192)
+//--
+//++
+//
+// Links	Parents
+//
+//    SphericalMap 
+//-- 
+//++
+//
+// Links	Descendants
+//
+//     
+//-- 
+
+/* --Methode-- */
+//++
+// Titre	Constructeurs
+//--
+//++
+
+SphereGorski::SphereGorski()
+
+//--
+{
+  InitNul();
+}
+
+/* --Methode-- */
+//++
+SphereGorski::SphereGorski(char* flnm)
+
+//    Constructeur : charge une image à partir d'un fichier 
+//--
+{
+  InitNul();
+  PInPersist s(flnm);
+  Read(s);
+}
+
+//++
+SphereGorski::SphereGorski(int_4 m)
+
+//    Constructeur : m est la variable nside de l'algorithme de Gorski
+//    le nombre total de pixels sera Npix =  12*nside**2
+//    m DOIT OBLIGATOIREMENT ETRE UNE PUISSANCE DE 2 (<= 8192)
+//--
+{
+  //printf(" initialisation par defaut SphereGorski \n");
+  if (m<=0 || m> 8192) {
+    cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
+    exit(1);
+  }
+// verifier que m est une puissance de deux 
+  int_4 x=m;
+  while (x%2==0) x/=2;
+  if (x!=1) {  
+    cout << "SphereGorski : m doit etre une puissance de deux, m= " << m << endl;
+    exit(1);
+  }
+
+  InitNul();
+  Pixelize(m);
+  if (pix2x_==NULL) {
+    pix2x_=new int[1024];
+    pix2y_=new int[1024];
+    mk_pix2xy(pix2x_,pix2y_); 
+  }
+  if (x2pix_==NULL) {
+    x2pix_=new int[128];
+    y2pix_=new int[128];
+    mk_xy2pix(x2pix_,y2pix_); 
+  }
+}
+
+
+
+
+
+
+//++
+// Titre	Destructeur
+//--
+//++
+SphereGorski::~SphereGorski()
+
+//--
+{
+  Clear();
+}
+
+//++
+// Titre	Méthodes
+//--
+
+
+
+
+void  SphereGorski::Pixelize( int_4 m) 
+
+//    prépare la pixelisation Gorski (m a la même signification 
+//    que pour le constructeur)
+//
+//    
+//--
+{
+  
+  // On memorise les arguments d'appel
+  nSide_ = m;  
+  
+  
+  // Nombre total de pixels sur la sphere entiere
+  nPix_=12*nSide_*nSide_;
+
+// pour le moment les tableaux qui suivent seront ranges dans l'ordre 
+// de l'indexation GORSKY "RING"
+// on pourra ulterieurement changer de strategie et tirer profit
+// de la dualite d'indexation GORSKY (RING er NEST) : tout dependra 
+// de pourquoi c'est faire
+
+  // Creation et initialisation du vecteur des contenus des pixels 
+  mPix_ = new r_8[nPix_];
+  for(int i=0; i<nPix_; i++)  mPix_[i] = 0.;
+
+  // solid angle per pixel   
+  omeg_=4*Pi/nPix_;
+}
+
+void SphereGorski::InitNul()
+//
+//    initialise à zéro les variables de classe 
+{
+  nlmax_=0;
+  nmmax_=0;
+  iseed_=0;
+  fwhm_=0.;
+  quadrupole_=0.;
+  sym_cut_deg_=0.;
+  for (int k=0; k<128;k++) powFile_[k]=' ';
+  nSide_=0;
+  nPix_ =0;
+  mPix_ = NULL;
+  pix2x_=NULL;
+  pix2y_=NULL;
+  x2pix_=NULL;
+  y2pix_=NULL;
+  pix2xy_=NULL;
+}
+/* --Methode-- */
+void SphereGorski::Clear()
+
+{
+  if (mPix_)   delete[] mPix_;
+  if (pix2x_)  delete[] pix2x_;
+  if (pix2y_)  delete[] pix2y_;
+  if (x2pix_)  delete[] x2pix_;
+  if (y2pix_)  delete[] y2pix_;
+  if (pix2xy_) delete pix2xy_;
+}
+//++
+void SphereGorski::WriteSelf(POutPersist& s) const
+
+//    créer un fichier image
+//--
+{
+  char strg[256];
+  if (mInfo_) sprintf(strg, "SphereGorski: NSlices=%6d  NPix=%9d HasInfo",  
+                   (int_4)nSide_, (int_4)nPix_);
+  else
+    sprintf(strg, "SphereGorski: nSide=%6d  nPix=%9d ", 
+	  (int_4)nSide_, (int_4)nPix_);
+  s.PutLine(strg);
+  if (mInfo_)  mInfo_->Write(s);
+  s.PutI4(nlmax_);
+  s.PutI4(nmmax_);
+  s.PutI4(iseed_);
+  s.PutI4(nSide_);
+  s.PutI4(nPix_);
+  s.PutR4(fwhm_);
+  s.PutR4(quadrupole_);
+  s.PutR4(sym_cut_deg_);
+  s.PutR8(omeg_);
+  s.PutR8s(mPix_, nPix_);
+  s.PutLine(powFile_);
+  return;
+}
+
+/* --Methode-- */
+//++
+void SphereGorski::ReadSelf(PInPersist& s)
+
+//    relit un fichier d'image
+//--
+{
+  Clear();
+  char strg[256];
+  s.GetLine(strg, 255);
+// Pour savoir s'il y avait un DVList Info associe
+  bool hadinfo = false;
+  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
+  if (hadinfo) {    // Lecture eventuelle du DVList Info
+    if (mInfo_ == NULL)  mInfo_ = new DVList;
+    mInfo_->Read(s);
+  }
+  s.GetI4(nlmax_);
+  s.GetI4(nmmax_);
+  s.GetI4(iseed_);
+  s.GetI4(nSide_);
+  s.GetI4(nPix_);
+  s.GetR4(fwhm_);
+  s.GetR4(quadrupole_);
+  s.GetR4(sym_cut_deg_);
+  mPix_=new r_8[nPix_];
+  s.GetR8(omeg_);
+  s.GetR8s(mPix_, nPix_);
+  s.GetLine(powFile_, 127);
+
+  return;
+}
+
+
+/* --Methode-- */
+//++
+void SphereGorski::ReadFits(char flnm[])
+
+//    remplit la sphere a partir d'un fichier FITS
+//--
+{
+  strip(flnm,'B',' ');
+  if (access(flnm,F_OK) != 0) {perror(flnm); exit(1);}
+  if(!nPix_) {
+    cout << " ReadFits : SphereGorski non pixelisee " << endl;
+    exit(1);
+  }
+  int npixtot=nPix_;
+
+// quand map et mPix_ auront le meme type, map ne sera plus necessaire
+
+  float* map=new float[12*nSide_*nSide_];
+  input_map_(flnm,map,npixtot);
+  // Remplissage de la sphère 
+
+  for(int  j=0; j<nPix_; j++ ) mPix_[j]=(r_8)map[j];
+  delete [] map;
+  
+
+}
+
+/* --Methode-- */
+//++
+void SphereGorski::WriteFits(char flnm[])
+
+//    ecrit la sphere sur un fichier FITS
+
+//--
+{
+  strip(flnm,'B',' ');
+  if (access(flnm,F_OK) == 0) {
+    cout << " SphereGorski::WriteFits : le fichier existe deja" << endl; 
+    exit(1);
+  }
+  //else cout << "un fichier sera cree " << endl;
+  if(!nPix_) {
+    cout << " WriteFits : SphereGorski non pixelisee " << endl;
+    exit(1);
+  }
+
+
+  char infile[128];
+  for (int k=0; k< 128; k++) infile[k]=powFile_[k];
+
+// quand map et mPix_ auront le meme type, map ne sera plus necessaire
+
+  float* map=new float[12*nSide_*nSide_];
+
+  for(int  j=0; j<nPix_; j++ ) map[j]=(float)mPix_[j];
+  int nlmax=nlmax_;
+  int nsmax=nSide_;
+  int nmmax=nmmax_;
+  int iseed=iseed_;
+  float fwhm=fwhm_;
+  float quadrupole=quadrupole_;
+  ecrire_fits_(nlmax,nsmax,nmmax,map,infile,flnm,iseed,fwhm,quadrupole);
+  delete [] map;
+  
+}
+
+/* --Methode-- */
+//++
+int_4 SphereGorski::NbPixels() const
+
+//    Retourne le nombre de pixels du découpage
+//--
+{
+  return(nPix_);
+}
+
+//++
+int_4 SphereGorski::NbThetaSlices() const 
+
+//    Retourne le nombre de tranches en theta sur la sphere
+//--
+{return int_4(4*nSide_-1);}
+
+
+//++
+void  SphereGorski::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const 
+
+//    Retourne, pour la tranche en theta d'indice 'index' le theta 
+//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
+//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
+//    correspondantes
+//--
+{
+
+cout << "entree GetThetaSlice, couche no " << index << endl;
+  if (index<0 || index > NbThetaSlices()) {
+    // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
+    cout << " SphereGorski::GetThetaSlice : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  int_4 iring=0;
+  int lring=0;
+  if (index<nSide_-1) {
+    iring=2*index*(index+1);
+    lring=4*(index+1);
+  }else
+    if (index<3*nSide_) {
+      iring=2*nSide_*(2*index-nSide_+1);
+      lring=4*nSide_;
+    }else
+      {
+	int nc=4*nSide_-1-index;
+	iring=nPix_-2*nc*(nc+1);
+	lring=4*nc;
+      }
+  phi.Realloc(lring);
+  value.Realloc(lring);
+  float T=0.;
+  float F=0.;
+  for (int kk=0; kk<lring;kk++) {
+    PixThetaPhi(kk+iring,T,F);
+    phi(kk)=F;
+    value(kk)=PixVal(kk+iring);
+  }
+  theta=T;
+
+}
+
+
+
+/* --Methode-- */
+//++
+r_8& SphereGorski::PixVal(int_4 k)
+
+//    Retourne la valeur du contenu du pixel d'indice "RING" k 
+//--
+{
+  if ( (k<0) || (k >= nPix_) ) {
+    // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
+    cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  return(mPix_[k]);
+}
+
+/* --Methode-- */
+//++
+r_8 const& SphereGorski::PixVal(int_4 k) const
+
+//    Retourne la valeur du contenu du pixel d'indice "RING" k 
+//--
+{
+  if ( (k<0) || (k >= nPix_) ) {
+    //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range")); 
+    cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  return(mPix_[k]);
+}
+
+
+//++
+r_8& SphereGorski::PixValNest(int_4 k)
+
+//    Retourne la valeur du contenu du pixel d'indice "NESTED" k 
+//--
+{
+  if ( (k<0) || (k >= nPix_) ) { 
+   //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range")); 
+    cout << " SphereGorski::PIxValNest : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  return mPix_[nest2ring(nSide_,k)];
+}
+//++
+
+r_8 const& SphereGorski::PixValNest(int_4 k) const
+
+//    Retourne la valeur du contenu du pixel d'indice "NESTED" k 
+//--
+{
+  if ( (k<0) || (k >= nPix_) ) { 
+   //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range")); 
+    cout << " SphereGorski::PIxValNest : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  int_4 pix=nest2ring(nSide_,k);
+  return mPix_[pix];
+}
+
+
+
+/* --Methode-- */
+//++
+int_4 SphereGorski::PixIndexSph(r_4 theta, r_4 phi) const
+
+//    Retourne l'indice "RING" du pixel vers lequel pointe une direction 
+//    définie par ses coordonnées sphériques
+//--
+{
+  return ang2pix_ring(nSide_, double(theta), double(phi));
+}
+
+//++
+int_4 SphereGorski::PixIndexSphNest(r_4 theta, r_4 phi) const
+
+//    Retourne l'indice NESTED" du pixel vers lequel pointe une direction 
+//    définie par ses coordonnées sphériques
+//--
+{
+  return ang2pix_nest(nSide_, double(theta), double(phi));
+}
+
+
+/* --Methode-- */
+//++
+void SphereGorski::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const
+
+//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice "RING" k
+//--
+{
+  double t;
+  double p;
+  pix2ang_ring(nSide_,k, t,  p);
+  teta=(r_4)t;
+  phi=(r_4)p;
+}
+
+//++
+r_8       SphereGorski::PixSolAngle(int_4 dummy) const
+//    Pixel Solid angle  (steradians)
+//    All the pixels have the same solid angle. The dummy argument is
+//    for compatibility with eventual pixelizations which would not 
+//    fulfil this requirement.
+//--
+{
+  return omeg_;
+}
+
+
+//++
+void SphereGorski::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const
+
+//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice 
+//    NESTED k
+//--
+{   
+  double t;
+  double p;
+  pix2ang_nest(nSide_, k, t, p);
+  teta=(r_4)t;
+  phi=(r_4)p;
+}
+
+//++
+int_4 SphereGorski::NestToRing(int_4 k) 
+
+//    conversion d'index NESTD en un index RING
+//
+//--
+{
+  return  nest2ring(nSide_,k);
+}
+//++
+int_4 SphereGorski::RingToNest(int_4 k) 
+//
+//    conversion d'index RING en un index NESTED
+//
+//--
+{
+  return  ring2nest(nSide_,k);
+}
+//++
+void SphereGorski::anharm(int nlmax, float sym_c,float* powspec) 
+//
+//    analyse en harmoniques spheriques des valeurs des pixels de la 
+//    sphere : appel du module anafast (Gorski-Hivon)
+//
+//    "nlmax"   : multipole maximum, nlmax <= 2*nsmax (cf. Nyquist)
+//
+//    "sym c"   : coupure symetrique autour de l'equateur (degres)
+//
+//    "powspec" : tableau resultat (a reserver avant l'appel) de C(l) 
+//    (spectre de puissance)
+//
+//--
+//
+// Pb a resoudre : dans cette classe les valeurs de pixel sont "double"
+// dans anafast le tableau correspondant est "float"
+// pour l'instant on duplique les tableaux, il faudra decider quelque chose
+//
+{ 
+  if (nlmax > 2*nSide_) {
+    cout << " anharm : nlmax= " << nlmax << 
+            " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
+    exit(1);
+  }
+  else {
+    nlmax_=nlmax;
+    nmmax_=nlmax_;
+  }
+  sym_cut_deg_=sym_c;
+  float* map=new float[nPix_];
+  for (int k=0; k<nPix_; k++) map[k]=(float)mPix_[k];
+  int nsmax=nSide_;
+  int nmmax=nmmax_;
+  double sc=(double)sym_cut_deg_;
+  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
+  if (powspec==NULL) {
+    
+    cout << 
+      " anharm : un tableau de C_l doit etre alloue avant appel " << endl;
+    exit(1);
+  }
+  float* phas_n=new float[2*(nmmax+1)];
+  float* phas_s=new float[2*(nmmax+1)];
+  float*  dataw =new float[16*nsmax];
+  float*  work =new float[16*nsmax];
+
+  anafast_(nsmax,nlmax,nmmax,sc,map,alm_T, powspec,phas_n,phas_s,dataw,work);
+  quadrupole_=powspec[2];
+  delete [] map;
+  delete [] alm_T;
+  delete [] phas_n;
+  delete [] phas_s;
+  delete [] dataw;
+  delete [] work;
+  
+}
+
+
+//++
+void SphereGorski::synharm(int nlmax, int iseed,float fwhm, float* powspec) 
+//
+//    synthese  des valeurs des pixels de la sphere par l'intermediaire 
+//    des coefficients en harmoniques spheriques reconstitues apartir d'un 
+//    spectre en puissance : appel du module synfast (Gorski-Hivon)
+//
+//    powspec est un tableau (a fournir) de C(l) (spectre de puissance)
+//    Ce tableau doit contenir les valeur de C(l) par ordre 
+//    SEQUENTIEL de l (de l=0 a l=nlmax). IL SERA MODIFIE PAR L'ALGORITHME
+//
+//    nlmax : multipole maximum (nlmax <= 2*nsmax (cf. Nyquist)
+//    iseed : initialisation generation aleatoire (negatif, suggere : -1)
+//    fwhm  : largeur totale a mi-hauteur (minutes d'arc, >=0, ex: 5)
+//--
+// Pb a resoudre : dans cette classe les valeurs de pixel sont "double"
+// dans anafast le tableau correspondant est "float"
+// pour l'instant on duplique les tableaux, il faudra decider quelque chose
+
+{
+  if (nlmax > 2*nSide_) {
+    cout << " sphereGorski::synharm: nlmax= " << nlmax << 
+            " doit etre <= 2*nsmax (cf. Nyquist), soit : " << 2*nSide_ << endl;
+    exit(1);
+  }
+  else {
+    nlmax_=nlmax;
+    nmmax_=nlmax_;
+    quadrupole_=powspec[2];
+  }
+  if (powspec==NULL) {
+    
+    cout << 
+      "sphereGorski::synharm : un tableau de C_l doit etre alloue avant appel"
+	 << endl;
+    exit(1);
+  }
+  iseed_=iseed;
+  fwhm_ =fwhm; 
+  float* map=new float[nPix_];
+  int nsmax=nSide_;
+  int nmmax=nmmax_;
+  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
+  
+
+//    tableaux de travail 
+  double* b_north=new double[2*(2*nmmax+1)];
+  double* b_south=new double[2*(2*nmmax+1)];
+  double* bw=new double[2*4*nsmax];
+  double* data=new double[2*4*nsmax];
+  double* work=new double[2*4*nsmax];
+  float* lread=new float[nlmax+1];
+  synfast_(nsmax,nlmax,nmmax,iseed,fwhm, map,alm_T, powspec,
+                  b_north,b_south,bw,data,work,lread);
+  for (int k=0; k<nPix_; k++) mPix_[k]=(double)map[k];
+  delete [] map;
+  delete [] alm_T;
+  delete [] b_north;
+  delete [] b_south;
+  delete [] bw; 
+  delete [] data;
+  delete [] work;
+  delete [] lread;
+
+}
+
+int SphereGorski::nest2ring(int nside, int ipnest) const {
+  /*
+    c=======================================================================
+    subroutine nest2ring(nside, ipnest, ipring)
+    c=======================================================================
+    c     conversion from NESTED to RING pixel number
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+      int npix, npface, face_num, ncap, n_before;
+      int ipf, ip_low, ip_trunc, ip_med, ip_hi;
+      int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
+      int ns_max=8192;
+      int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
+      int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
+
+      if(  nside<1 ||  nside>ns_max ) {
+	cout << "nside out of range" << endl;
+	exit(0);
+      }
+      npix = 12 *  nside* nside;
+      if( ipnest<0 || ipnest>npix-1 ) {
+	cout << "ipnest out of range" << endl;
+	exit(0);
+      }
+
+
+      ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
+      nl4   = 4* nside;
+
+      //c     finds the face, and the number in the face
+      npface =  nside* nside;
+      //cccccc      ip = ipnest - 1         ! in {0,npix-1}
+
+      face_num = ipnest/npface;//  ! face number in {0,11}
+      ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
+      //c     finds the x,y on the face (starting from the lowest corner)
+      //c     from the pixel number
+      ip_low=ipf%1024;                //   ! content of the last 10 bits
+      ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
+      ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
+      ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
+
+      ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
+      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
+
+      //c     transforms this in (horizontal, vertical) coordinates
+      jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
+      jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
+
+      //c     computes the z coordinate on the sphere
+      //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
+      jr =  jrll[face_num]*nside - jrt - 1;
+      nr = nside;//                  ! equatorial region (the most frequent)
+      n_before = ncap + nl4 * (jr - nside);
+      kshift=(jr - nside)%2;
+      if( jr<nside ) {//then     ! north pole region
+         nr = jr;
+         n_before = 2 * nr * (nr - 1);
+         kshift = 0;
+      }
+      else if( jr>3*nside ) {//then ! south pole region
+         nr = nl4 - jr;
+         n_before = npix - 2 * (nr + 1) * nr;
+         kshift = 0;
+      }
+
+      //c     computes the phi coordinate on the sphere, in [0,2Pi]
+      jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
+
+      if( jp>nl4 ) jp = jp - nl4;
+      if( jp<1 )   jp = jp + nl4;
+
+      int aux=n_before + jp - 1;
+      return (n_before + jp - 1);// ! in {0, npix-1}
+
+}
+
+void SphereGorski::mk_pix2xy(int *pix2x,int *pix2y) {
+  /*
+    c=======================================================================
+    subroutine mk_pix2xy
+    c=======================================================================
+    c     constructs the array giving x and y in the face from pixel number
+    c     for the nested (quad-cube like) ordering of pixels
+    c
+    c     the bits corresponding to x and y are interleaved in the pixel number
+    c     one breaks up the pixel number by even and odd bits
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+      int kpix, jpix, IX, IY, IP, ID;
+      for (int i=0;i<1023;i++ ) pix2x[i]=0;
+
+      for( kpix=0;kpix<1024;kpix++ ) {
+	jpix = kpix;
+	IX = 0;
+	IY = 0;
+	IP = 1 ;//              ! bit position (in x and y)
+	while( jpix!=0 ){// ! go through all the bits
+	  ID=jpix%2;//  ! bit value (in kpix), goes in ix
+	  jpix = jpix/2;
+	  IX = ID*IP+IX;
+	  
+	  ID=jpix%2;//  ! bit value (in kpix), goes in iy
+	  jpix = jpix/2;
+	  IY = ID*IP+IY;
+	  
+	  IP = 2*IP;//         ! next bit (in x and y)
+	}
+	
+	pix2x[kpix] = IX;//     ! in 0,31
+	pix2y[kpix] = IY;//     ! in 0,31
+      }
+}
+
+int SphereGorski::ring2nest(int nside, int ipring) const {
+  /*
+    c=======================================================================
+    subroutine ring2nest(nside, ipring, ipnest)
+    c=======================================================================
+    c     conversion from RING to NESTED pixel number
+    c=======================================================================
+    */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  double fihip, hip;
+  int npix, nl2, nl4, ncap, ip, iphi, ipt, ipring1;
+  int     kshift, face_num, nr;
+  int irn, ire, irm, irs, irt, ifm , ifp;
+  int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf;
+  int ns_max(8192);
+  
+  // coordinate of the lowest corner of each face
+  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};// ! in unit of nside
+  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};//! in unit of nside/2
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12 * nside*nside;
+  if( ipring<0 || ipring>npix-1 ) {
+    cout << "ipring out of range" << endl;
+    exit(0);
+  }
+  
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  npix = 12*nside*nside;//      ! total number of points
+  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
+  ipring1 = ipring + 1;
+  
+  //c     finds the ring number, the position of the ring and the face number
+  if( ipring1<=ncap ) {//then
+    
+    hip   = ipring1/2.;
+    fihip = (int)floor ( hip );
+    irn   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
+    iphi  = ipring1 - 2*irn*(irn - 1);
+    
+    kshift = 0;
+    nr = irn   ;//               ! 1/4 of the number of points on the current ring
+    face_num = (iphi-1) / irn;// ! in {0,3}
+  }
+  else if( ipring1<=nl2*(5*nside+1) ) {//then
+    
+    ip    = ipring1 - ncap - 1;
+    irn   = (int)floor( ip / nl4 ) + nside;//               ! counted from North pole
+    iphi  = (int)fmod(ip,nl4) + 1;
+    
+    kshift  = (int)fmod(irn+nside,2);//  ! 1 if irn+nside is odd, 0 otherwise
+    nr = nside;
+    ire =  irn - nside + 1;// ! in {1, 2*nside +1}
+    irm =  nl2 + 2 - ire;
+    ifm = (iphi - ire/2 + nside -1) / nside;// ! face boundary
+    ifp = (iphi - irm/2 + nside -1) / nside;
+    if( ifp==ifm ) {//then          ! faces 4 to 7
+      face_num = (int)fmod(ifp,4) + 4;
+    }
+    else if( ifp + 1==ifm ) {//then ! (half-)faces 0 to 3
+      face_num = ifp;
+    }
+    else if( ifp - 1==ifm ) {//then ! (half-)faces 8 to 11
+      face_num = ifp + 7;
+    }
+  }
+  else {
+    
+    ip    = npix - ipring1 + 1;
+    hip   = ip/2.;
+    fihip = floor ( hip );
+    irs   = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//  ! counted from South pole
+    iphi  = 4*irs + 1 - (ip - 2*irs*(irs-1));
+    
+    kshift = 0;
+    nr = irs;
+    irn   = nl4 - irs;
+    face_num = (iphi-1) / irs + 8;// ! in {8,11}
+  }
+  
+  //c     finds the (x,y) on the face
+  irt =   irn  - jrll[face_num]*nside + 1;//       ! in {-nside+1,0}
+  ipt = 2*iphi - jpll[face_num]*nr - kshift - 1;// ! in {-nside+1,nside-1}
+
+
+  if( ipt>=nl2 ) ipt = ipt - 8*nside;// ! for the face #4
+  
+  ix =  (ipt - irt ) / 2;
+  iy = -(ipt + irt ) / 2;
+  
+  ix_low = (int)fmod(ix,128);
+  ix_hi  = ix/128;
+  iy_low = (int)fmod(iy,128);
+  iy_hi  = iy/128;
+  ipf =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)
+    + (x2pix_[ix_low]+y2pix_[iy_low]);
+
+  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
+}
+
+void SphereGorski::mk_xy2pix(int *x2pix, int *y2pix) {
+  /*
+    c=======================================================================
+    subroutine mk_xy2pix
+    c=======================================================================
+    c     sets the array giving the number of the pixel lying in (x,y)
+    c     x and y are in {1,128}
+    c     the pixel number is in {0,128**2-1}
+    c
+    c     if  i-1 = sum_p=0  b_p * 2^p
+    c     then ix = sum_p=0  b_p * 4^p
+    c          iy = 2*ix
+    c     ix + iy in {0, 128**2 -1}
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  int K,IP,I,J,ID;
+
+  for( int i(0);i<127;i++ ) x2pix[i]=0;
+  for( I=1;I<=128;I++ ) {
+    J  = I-1;//            !pixel numbers
+    K  = 0;//
+    IP = 1;//
+    truc : if( J==0 ) {
+      x2pix[I-1] = K;
+      y2pix[I-1] = 2*K;
+    }
+    else {
+      ID = (int)fmod(J,2);
+      J  = J/2;
+      K  = IP*ID+K;
+      IP = IP*4;
+      goto truc;
+    }
+  }     
+  //c      endif
+  
+}
+
+int SphereGorski::ang2pix_ring(int nside, double theta, double phi) const {
+    /*
+    c=======================================================================
+    c     gives the pixel number ipix (RING) 
+    c     corresponding to angles theta and phi
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+
+  int nl2, nl4, ncap, npix, jp, jm, ipix1;
+  double  z, za, tt, tp, tmp;
+  int ir, ip, kshift;
+
+  double piover2(Pi/2.);
+  double twopi(2.*Pi);
+  double z0(2./3.);
+  int ns_max(8192);
+
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+
+  if( theta<0. || theta>Pi) {
+    cout << "theta out of range" << endl;
+    exit(0);
+  }
+
+  z = cos(theta);
+  za = fabs(z);
+  if( phi >= twopi)  phi = phi - twopi;
+  if (phi < 0.)     phi = phi + twopi;
+  tt = phi / piover2;//  ! in [0,4)
+
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  ncap  = nl2*(nside-1);// ! number of pixels in the north polar cap
+  npix  = 12*nside*nside;
+
+  if( za <= z0 ) {
+    
+    jp = (int)floor(nside*(0.5 + tt - z*0.75));// ! index of  ascending edge line 
+    jm = (int)floor(nside*(0.5 + tt + z*0.75));// ! index of descending edge line
+
+    ir = nside + 1 + jp - jm;// ! in {1,2n+1} (ring number counted from z=2/3)
+    kshift = 0;
+    if (fmod(ir,2)==0.) kshift = 1;// ! kshift=1 if ir even, 0 otherwise
+    
+    ip = (int)floor( ( jp+jm - nside + kshift + 1 ) / 2 ) + 1;// ! in {1,4n}
+    if( ip>nl4 ) ip = ip - nl4;
+    
+    ipix1 = ncap + nl4*(ir-1) + ip ;
+  }
+  else {
+
+    tp = tt - floor(tt);//      !MOD(tt,1.d0)
+    tmp = sqrt( 3.*(1. - za) );
+
+    jp = (int)floor( nside * tp * tmp );// ! increasing edge line index
+    jm = (int)floor( nside * (1. - tp) * tmp );// ! decreasing edge line index
+
+    ir = jp + jm + 1;//        ! ring number counted from the closest pole
+    ip = (int)floor( tt * ir ) + 1;// ! in {1,4*ir}
+    if( ip>4*ir ) ip = ip - 4*ir;
+      
+    ipix1 = 2*ir*(ir-1) + ip;
+      if( z<=0. ) {
+	ipix1 = npix - 2*ir*(ir+1) + ip;
+      }
+  }
+    return (ipix1 - 1);// ! in {0, npix-1}
+  
+}
+
+int SphereGorski::ang2pix_nest(int nside, double theta, double phi) const {
+  /*
+    c=======================================================================
+    subroutine ang2pix_nest(nside, theta, phi, ipix)
+    c=======================================================================
+    c     gives the pixel number ipix (NESTED) 
+    c     corresponding to angles theta and phi
+    c
+    c     the computation is made to the highest resolution available (nside=8192)
+    c     and then degraded to that required (by integer division)
+    c     this doesn't cost more, and it makes sure 
+    c     that the treatement of round-off will be consistent 
+    c     for every resolution
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+    
+      double    z, za, z0, tt, tp, tmp;
+      int face_num,jp,jm;
+      int ifp, ifm;
+      int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
+      double piover2(Pi/2.), twopi(2.*Pi);
+      int ns_max(8192);
+
+      if( nside<1 || nside>ns_max ) {
+	cout << "nside out of range" << endl;
+	exit(0);
+      }
+      if( theta<0 || theta>Pi ) {
+	cout << "theta out of range" << endl;
+	exit(0);
+      }
+      z  = cos(theta);
+      za = fabs(z);
+      z0 = 2./3.;
+      if( phi>=twopi ) phi = phi - twopi;
+      if( phi<0. )    phi = phi + twopi;
+      tt = phi / piover2;// ! in [0,4[
+      if( za<=z0 ) { // then ! equatorial region
+	
+	//(the index of edge lines increase when the longitude=phi goes up)
+	jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
+	jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
+
+	//c        finds the face
+	ifp = jp / ns_max;//  ! in {0,4}
+	ifm = jm / ns_max;
+        if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then          ! faces 4 to 7
+	else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
+	else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
+
+	ix = (int)fmod(jm, ns_max);
+	iy = ns_max - (int)fmod(jp, ns_max) - 1;
+      }
+      else { //! polar region, za > 2/3
+	     
+	ntt = (int)floor(tt);
+	if( ntt>=4 ) ntt = 3;
+	tp = tt - ntt;
+	tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
+
+	//(the index of edge lines increase when distance from the closest pole goes up)
+	jp = (int)floor( ns_max * tp          * tmp );// ! line going toward the pole as phi increases
+	jm = (int)floor( ns_max * (1. - tp) * tmp );// ! that one goes away of the closest pole
+	jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
+	jm = (int)min(ns_max-1, jm);
+
+	 // finds the face and pixel's (x,y)
+         if( z>=0 ) {
+	   face_num = ntt;//  ! in {0,3}
+            ix = ns_max - jm - 1;
+            iy = ns_max - jp - 1;
+	 }
+         else {
+	   face_num = ntt + 8;// ! in {8,11}
+            ix =  jp;
+            iy =  jm;
+	 }
+      }
+
+      ix_low = (int)fmod(ix,128);
+      ix_hi  =     ix/128;
+      iy_low = (int)fmod(iy,128);
+      iy_hi  =     iy/128;
+      ipf =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)+ (x2pix_[ix_low]+y2pix_[iy_low]);
+      ipf = ipf / pow(ns_max/nside,2);//  ! in {0, nside**2 - 1}
+      return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
+}
+
+
+void SphereGorski::pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const {
+  /*
+    c=======================================================================
+    c     gives theta and phi corresponding to pixel ipix (RING) 
+    c     for a parameter nside
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+  
+  int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
+  double  fact1, fact2, fodd, hip, fihip;
+  
+  int ns_max(8192);
+  
+  if( nside<1 || nside>ns_max ) {
+    cout << "nside out of range" << endl;
+    exit(0);
+  }
+  npix = 12*nside*nside;      // ! total number of points
+  if( ipix<0 || ipix>npix-1 ) {
+    cout << "ipix out of range" << endl;
+    exit(0);
+  }
+  
+  ipix1 = ipix + 1; // in {1, npix}
+  nl2 = 2*nside;
+  nl4 = 4*nside;
+  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
+  fact1 = 1.5*nside;
+  fact2 = 3.0*nside*nside;
+  
+  if( ipix1 <= ncap ) {  //! North Polar cap -------------
+    
+    hip   = ipix1/2.;
+    fihip = floor(hip);
+    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
+    iphi  = ipix1 - 2*iring*(iring - 1);
+    
+    theta = acos( 1. - iring*iring / fact2 );
+    phi   = (1.*iphi - 0.5) * Pi/(2.*iring);
+    //    cout << theta << " " << phi << endl;
+  }
+  else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
+    
+    ip    = ipix1 - ncap - 1;
+    iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
+    iphi  = (int)fmod(ip,nl4) + 1;
+    
+    fodd  = 0.5 * (1 + fmod((double)(iring+nside),2));//  ! 1 if iring+nside is odd, 1/2 otherwise
+    theta = acos( (nl2 - iring) / fact1 );
+    phi   = (1.*iphi - fodd) * Pi /(2.*nside);
+  }
+  else {//! South Polar cap -----------------------------------
+    
+    ip    = npix - ipix1 + 1;
+    hip   = ip/2.;
+    fihip = 1.*hip;
+    iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;//     ! counted from South pole
+    iphi  = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1)));
+    
+    theta = acos( -1. + iring*iring / fact2 );
+    phi   = (1.*iphi - 0.5) * Pi/(2.*iring);
+    //    cout << theta << " " << phi << endl;
+  }
+}
+
+void SphereGorski::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
+  /*
+    c=======================================================================
+    subroutine pix2ang_nest(nside, ipix, theta, phi)
+    c=======================================================================
+    c     gives theta and phi corresponding to pixel ipix (NESTED) 
+    c     for a parameter nside
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+    
+      int npix, npface, face_num;
+      int  ipf, ip_low, ip_trunc, ip_med, ip_hi;
+      int     ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
+      double z, fn, fact1, fact2;
+      double piover2(Pi/2.);
+      int ns_max(8192);
+      
+      
+      // ! coordinate of the lowest corner of each face
+      int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
+      int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
+      
+      if( nside<1 || nside>ns_max ) {
+	cout << "nside out of range" << endl;
+	exit(0);
+      }
+      npix = 12 * nside*nside;
+      if( ipix<0 || ipix>npix-1 ) {
+	cout << "ipix out of range" << endl;
+	exit(0);
+      }
+
+      fn = 1.*nside;
+      fact1 = 1./(3.*fn*fn);
+      fact2 = 2./(3.*fn);
+      nl4   = 4*nside;
+
+      //c     finds the face, and the number in the face
+      npface = nside*nside;
+
+      face_num = ipix/npface;//  ! face number in {0,11}
+      ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
+
+      //c     finds the x,y on the face (starting from the lowest corner)
+      //c     from the pixel number
+      ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
+      ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
+      ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
+      ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
+
+      ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
+      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
+
+      //c     transforms this in (horizontal, vertical) coordinates
+      jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
+      jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
+
+      //c     computes the z coordinate on the sphere
+      //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
+      jr =  jrll[face_num]*nside - jrt - 1;
+      nr = nside;//                  ! equatorial region (the most frequent)
+      z  = (2*nside-jr)*fact2;
+      kshift = (int)fmod(jr - nside, 2);
+      if( jr<nside ) { //then     ! north pole region
+         nr = jr;
+         z = 1. - nr*nr*fact1;
+         kshift = 0;
+      }
+      else {
+	if( jr>3*nside ) {// then ! south pole region
+         nr = nl4 - jr;
+         z = - 1. + nr*nr*fact1;
+         kshift = 0;
+	}
+      }
+      theta = acos(z);
+      
+      //c     computes the phi coordinate on the sphere, in [0,2Pi]
+      //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
+      jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
+      if( jp>nl4 ) jp = jp - nl4;
+      if( jp<1 )   jp = jp + nl4;
+
+      phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
+
+}
+
+
+void SphereGorski::Pix2XY::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
+  /*
+    c=======================================================================
+    subroutine pix2ang_nest(nside, ipix, theta, phi)
+    c=======================================================================
+    c     gives theta and phi corresponding to pixel ipix (NESTED) 
+    c     for a parameter nside
+    c=======================================================================
+  */
+  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
+  //  (16/12/98)
+    
+      int npix, npface, face_num;
+      int  ipf, ip_low, ip_trunc, ip_med, ip_hi;
+      int     ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
+      double z, fn, fact1, fact2;
+      double piover2(Pi/2.);
+      int ns_max(8192);
+      
+      
+      // ! coordinate of the lowest corner of each face
+      int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
+      int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
+      
+      if( nside<1 || nside>ns_max ) {
+	cout << "nside out of range" << endl;
+	exit(0);
+      }
+      npix = 12 * nside*nside;
+      if( ipix<0 || ipix>npix-1 ) {
+	cout << "ipix out of range" << endl;
+	exit(0);
+      }
+
+      fn = 1.*nside;
+      fact1 = 1./(3.*fn*fn);
+      fact2 = 2./(3.*fn);
+      nl4   = 4*nside;
+
+      //c     finds the face, and the number in the face
+      npface = nside*nside;
+
+      face_num = ipix/npface;//  ! face number in {0,11}
+      ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
+
+      //c     finds the x,y on the face (starting from the lowest corner)
+      //c     from the pixel number
+      ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
+      ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
+      ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
+      ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
+
+      ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
+      iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
+
+      //c     transforms this in (horizontal, vertical) coordinates
+      jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
+      jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
+
+      //c     computes the z coordinate on the sphere
+      //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
+      jr =  jrll[face_num]*nside - jrt - 1;
+      nr = nside;//                  ! equatorial region (the most frequent)
+      z  = (2*nside-jr)*fact2;
+      kshift = (int)fmod(jr - nside, 2);
+      if( jr<nside ) { //then     ! north pole region
+         nr = jr;
+         z = 1. - nr*nr*fact1;
+         kshift = 0;
+      }
+      else {
+	if( jr>3*nside ) {// then ! south pole region
+         nr = nl4 - jr;
+         z = - 1. + nr*nr*fact1;
+         kshift = 0;
+	}
+      }
+      theta = acos(z);
+      
+      //c     computes the phi coordinate on the sphere, in [0,2Pi]
+      //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
+      jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
+      if( jp>nl4 ) jp = jp - nl4;
+      if( jp<1 )   jp = jp + nl4;
+
+      phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
+
+}
+
+
+
Index: trunk/SophyaLib/Samba/spheregorski.h
===================================================================
--- trunk/SophyaLib/Samba/spheregorski.h	(revision 228)
+++ trunk/SophyaLib/Samba/spheregorski.h	(revision 228)
@@ -0,0 +1,154 @@
+#ifndef SPHEREGORSKI_SEEN
+#define SPHEREGORSKI_SEEN
+
+#include "sphericalmap.h"
+#include "cvector.h"
+
+// les attributs de classe pix2x_ et pix2y_ sont relatifs a la traduction des 
+// indices  RING en indices NESTED (ou l'inverse). Ce sont des tableaux 
+// d'entiers initialises et remplis par le constructeur. Dans la version 
+// FORTRAN de healpix, il s'agissait de tableaux mis en COMMON. Ils etaient 
+// initialises au premier appel d'une conversion d'indice. Je ne peux pas 
+// garder cette procedure car, par exemple, la fonction PixValNest() const 
+// n'autorisera pas la constitution de ces tableaux (a cause du const). Ainsi, 
+// des la creation d'un objet SphereGorski ces tableaux sont calcules, ce qui 
+// est, certes, inutile si on ne se sert pas des indices NESTED. Ca ne me 
+// rejouit pas, mais c'est le seul moyen que j'ai trouve pour temir compte de 
+// toutes les demandes et des options prises.
+// 
+//                                                G. Le Meur 
+
+class SphereGorski : public SphericalMap {
+public :
+
+
+SphereGorski();
+SphereGorski(int_4 m);
+SphereGorski(char* flnm);
+virtual ~SphereGorski();
+
+// ------------ Persistence handling
+     enum {classId = 0xF002 };
+
+int_4                  ClassId() const        { return classId; }
+
+virtual void           WriteSelf(POutPersist&) const;
+virtual void           ReadSelf(PInPersist&);
+
+// -----------   FITS IO (1D FITS ARRAY)
+void ReadFits(char flnm[] );
+void WriteFits(char flnm[]);
+
+// ------------------ Definition of PixelMap abstract methods
+
+/* Nombre de pixels du decoupage                                         */
+virtual  int_4           NbPixels() const;
+
+/* Valeur du contenu du pixel d'indice "RING" k                                 */
+virtual  r_8&       PixVal(int_4 k);
+virtual  r_8  const&     PixVal(int_4 k) const;
+
+int_4 NbThetaSlices() const;
+void  GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const;
+
+/* Indice "RING" du pixel vers lequel pointe une direction definie par 
+ses  coordonnees spheriques */                                    
+virtual  int_4           PixIndexSph(float theta, float phi) const;
+
+/* Coordonnees spheriques du milieu du pixel d'indice "RING" k       */
+virtual  void          PixThetaPhi(int_4 k, float& teta, float& phi) const;
+
+// Pixel Solid angle  (steradians)
+virtual r_8         PixSolAngle(int_4 dummy) const; 
+
+
+// --------------- Specific methods 
+
+    // NEST indexing
+/* Valeur du contenu du pixel d'indice "NEST" k                                 */
+virtual  r_8&       PixValNest(int_4 k);
+virtual  r_8 const&       PixValNest(int_4 k) const;
+
+/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
+ses  coordonnees spheriques */                                    
+virtual  int_4      PixIndexSphNest(float theta, float phi) const;
+
+/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
+virtual  void       PixThetaPhiNest(int_4 k, float& teta, float& phi) const;
+
+   // algorithme de pixelisation
+void Pixelize(int_4); 
+
+
+/* convertit index nested en ring  */
+int_4 NestToRing(int_4 );
+
+
+/* convertit index ring en nested" */
+int_4 RingToNest(int_4 );
+
+
+/*    analyse en harmoniques spheriques des valeurs des pixels de la 
+   sphere : appel du module anafast (Gorski-Hivon) */
+void anharm(int, float, float*);
+ 
+/*synthese  des valeurs des pixels de la sphere par l'intermediaire 
+  des coefficients en harmoniques spheriques reconstitues apartir d'un 
+  spectre en puissance : appel du module synfast (Gorski-Hivon) */
+void synharm(int , int, float, float* ); 
+
+
+// ------------- méthodes internes ----------------------
+          
+private :
+
+void	      InitNul();
+void          Clear();
+int           nest2ring(int nside, int ipnest) const;
+void          mk_pix2xy(int *pix2x,int *pix2y);
+int           ring2nest(int nside, int ipring) const;
+void          mk_xy2pix(int *x2pix, int *y2pix);
+int           ang2pix_ring(int nside, double theta, double phi) const;
+int           ang2pix_nest(int nside, double theta, double phi) const;
+void          pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const;
+void          pix2ang_nest(int nside, int ipix, double& theta, double& phi) const;
+// ------------- variables internes -----------------------
+		   private :
+class Pix2XY {
+  public :
+Pix2XY() { 
+    pix2x_=new int[1024];
+    pix2y_=new int[1024];
+    mk_pix2xy(pix2x_,pix2y_); 
+  }
+~Pix2XY() {
+  if (pix2x_)  delete[] pix2x_;
+  if (pix2y_)  delete[] pix2y_;
+  }
+void          pix2ang_nest(int nside, int ipix, double& theta, double& phi) const;
+
+private : 
+int *pix2x_;
+int *pix2y_;
+void          mk_pix2xy(int *pix2x,int *pix2y);
+
+};
+
+int nlmax_;
+int nmmax_;
+int iseed_;
+int  nSide_;
+int_4  nPix_;
+int *pix2x_;
+int *pix2y_;
+int *x2pix_;
+int *y2pix_;
+float fwhm_;
+float quadrupole_;
+float sym_cut_deg_;
+r_8  omeg_;
+r_8*    mPix_;
+char powFile_[128];
+Pix2XY *pix2xy_;
+};
+#endif
Index: trunk/SophyaLib/Samba/spherethetaphi.cc
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 228)
+++ trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 228)
@@ -0,0 +1,635 @@
+//
+// 
+#include "spherethetaphi.h"
+#include "nbmath.h"
+#include <iostream.h>
+//***************************************************************
+//++ 
+// Class	SphereThetaPhi
+// 
+//
+// include      spherethetaphi.h nbmath.h 
+//
+//    Découpage de la sphère selon theta et phi, chaque 
+//    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du 
+//    beurre), chacune des m bandes de theta ainsi définies étant découpée par 
+//    des méridiens équirepartis, ce découpage étant fait de sorte que tous 
+//    les pixels aient la même surface et soient le plus carré possible.
+//    On commence par découper l'hémisphère de z positif en partant du pôle et 
+//    en allant vers l'équateur. Le premier pixel est la calotte polaire, 
+//    il est circulaire et centré sur theta=0.
+//--
+//++
+//
+// Links	Parents
+//
+//    SphericalMap 
+//-- 
+//++
+//
+// Links	Descendants
+//
+//     
+//-- 
+
+/* --Methode-- */
+//++
+// Titre	Constructeurs
+//--
+//++
+
+
+SphereThetaPhi::SphereThetaPhi()
+
+//--
+{
+  InitNul();
+}
+
+/* --Methode-- */
+//++
+SphereThetaPhi::SphereThetaPhi(char* flnm)
+
+//    Constructeur : charge une image à partir d'un fichier 
+//--
+{
+  InitNul();
+  PInPersist s(flnm);
+  Read(s);
+}
+
+/* --Methode-- */
+
+//++
+SphereThetaPhi::SphereThetaPhi(int_4 m, int_4 pet)
+
+//    Constructeur : m est le nombre de bandes en theta sur un hémisphère 
+//    (la calotte constituant la premiere bande).
+//    pet est le nombre de pixels (pétales) de la bande en contact avec la 
+//    calotte polaire. Pour l'instant pet est inopérant! 
+//--
+{
+  InitNul();
+  Pixelize(m,pet);
+}
+
+
+//++
+// Titre	Destructeur
+//--
+//++
+SphereThetaPhi::~SphereThetaPhi()
+
+//--
+{
+  Clear();
+}
+
+//++
+// Titre	Méthodes
+//--
+void SphereThetaPhi::InitNul()
+//
+//    initialise à zéro les variables de classe pointeurs
+{
+  mTheta = NULL;
+  mNphi = NULL;
+  mTNphi = NULL;
+  // mPix = NULL;
+  _pixel=NULL;
+  mNTheta = mNPet = 0;
+  mNPix = 0;
+}
+
+/* --Methode-- */
+void SphereThetaPhi::Clear()
+
+{
+  if (mTheta)  delete[] mTheta;
+  if (mNphi)  delete[] mNphi;
+  if (mTNphi) delete[] mTNphi;
+  //if (mPix)   delete[] mPix;
+  if (_pixel) _pixel->Detach();
+  mTheta = NULL;
+  mNphi = NULL;
+  mTNphi = NULL;
+  //mPix = NULL;
+  _pixel=NULL;
+  mNTheta = mNPet = 0;
+  mNPix = 0;
+}
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::WriteSelf(POutPersist& s) const
+
+//    créer un fichier image
+//--
+{
+  char strg[256];
+  if (mInfo_) {sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo", (int_4)mNTheta, (int_4)mNPix);
+  }
+  else { sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d ", 
+               (int_4)mNTheta, (int_4)mNPix);
+  }
+  s.PutLine(strg);
+  if (mInfo_)  mInfo_->Write(s);
+   s.PutI4(mNTheta);
+  s.PutI4(mNPet);
+  s.PutI4(mNPix);
+  s.PutR8(mOmeg);
+  s.PutR4s(mTheta, mNTheta+1);
+  s.PutI4s(mNphi, mNTheta);
+  s.PutI4s(mTNphi, mNTheta+1);
+  s.PutR8s(_pixel->Data(), mNPix);
+  //s.PutR8s(mPix, mNPix);
+  return;
+}
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::ReadSelf(PInPersist& s)
+
+//    relit un fichier d'image
+//--
+{
+  Clear();
+  char strg[256];
+  s.GetLine(strg, 255);
+// Pour savoir s'il y avait un DVList Info associe
+  bool hadinfo = false;
+  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
+  if (hadinfo) {    // Lecture eventuelle du DVList Info
+    if (mInfo_ == NULL)  mInfo_ = new DVList;
+    mInfo_->Read(s);
+  }
+  s.GetI4(mNTheta);
+  s.GetI4(mNPet);
+  s.GetI4(mNPix);
+  s.GetR8(mOmeg);
+  mTheta = new r_4[mNTheta+1];
+  mNphi = new int_4[mNTheta];
+  mTNphi = new int_4[mNTheta+1];
+  //mPix = new r_8[mNPix+1];
+  _pixel=new PDataArray<r_8>(mNPix,true);
+  s.GetR4s(mTheta, mNTheta+1);
+  s.GetI4s(mNphi, mNTheta);
+  s.GetI4s(mTNphi, mNTheta+1);
+  s.GetR8s(_pixel->Data(), mNPix);
+  //s.GetR8s(mPix, mNPix);
+  return;
+}
+
+
+/* --Methode-- */
+//++
+int_4 SphereThetaPhi::NbPixels() const
+
+//    Retourne le nombre de pixels du découpage
+//--
+{
+  return(mNPix);
+}
+
+
+static r_8 dummy_pixel = 0;
+
+/* --Methode-- */
+//++
+r_8& SphereThetaPhi::PixVal(int_4 k)
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if ( (k<0) || (k >= mNPix) ) {
+    //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+
+    //throw "SphereThetaPhi::PIxVal Pixel index out of range";
+  }
+  return (*(_pixel->Data()+k));
+  //return(mPix[k]);
+}
+//++
+r_8 const& SphereThetaPhi::PixVal(int_4 k) const
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if ( (k<0) || (k >= mNPix) ) {
+    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
+  //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+
+    throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
+  } 
+  return *(_pixel->Data()+k);
+  //return(mPix[k]);
+}
+
+
+/* --Methode-- */
+//++
+int_4 SphereThetaPhi::PixIndexSph(r_4 theta, r_4 phi) const
+
+//    Retourne l'indice du pixel vers lequel pointe une direction définie par 
+//    ses coordonnées sphériques
+//--
+{
+  r_4 dphi;
+  int_4 i,j,k;
+  bool fgzn = false;
+
+  if( (theta > (r_4)Pi) || (theta < 0. ) ) return(-1);
+  if( (phi < 0.) || (phi > DeuxPi) )  return(-1);
+  if( theta > (r_4)Pi*0.5)  { fgzn = true;  theta = Pi-theta; }
+  
+  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
+  // mTheta[kt] et mTheta[kt+1]
+  for( i=1; i< mNTheta; i++ ) 
+    if( theta < mTheta[i] ) break;
+  
+  dphi= (r_4)DeuxPi/(r_4)mNphi[i-1];
+  
+  if (fgzn)  k= mNPix-mTNphi[i]+(int_4)(phi/dphi);
+  else k= mTNphi[i-1]+(int_4)(phi/dphi); 
+  
+  return(k);
+}
+
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const
+
+//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
+//--
+{
+  int_4 i; 
+  bool fgzn = false;
+  
+  if( (k < 0) || (k >= mNPix) ) {theta = -99999.; phi = -99999.;  return; }
+  if( k >= mNPix/2)  {fgzn = true;  k = mNPix-1-k; }
+
+  // recupère l'indice i de la tranche qui contient le pixel k
+  for( i=0; i< mNTheta; i++ ) 
+    if( k < mTNphi[i+1] ) break;
+
+  // angle theta
+  theta= 0.5*(mTheta[i]+mTheta[i+1]);
+  if (fgzn) theta= Pi-theta;
+  
+  // angle phi
+  k -= mTNphi[i];
+  phi= (r_4)DeuxPi/(r_4)mNphi[i]*(r_4)(k+.5);
+  if (fgzn) phi= DeuxPi-phi;
+
+  return;
+}
+
+
+//++
+r_8         SphereThetaPhi::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 mOmeg;
+}
+
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::Limits(int_4 k, r_4& tetMin, r_4& tetMax, 
+			   r_4& phiMin, r_4& phiMax) 
+
+//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
+//--
+  {
+  int_4 j;
+  r_4 dphi;
+  bool fgzn= false;
+  
+  if( (k< 0) || (k >= mNPix) ) {
+    tetMin= -99999.; 
+    phiMin= -99999.;  
+    tetMax= -99999.; 
+    phiMax= -99999.;  
+    return;
+  }
+  
+  // si on se trouve dans l'hémisphère Sud
+  if( k >= mNPix/2 ) { 
+    fgzn= true; 
+    k= mNPix-1-k; 
+  }
+  
+  // on recupere l'indice i de la tranche qui contient le pixel k
+  int i;
+  for( i=0; i< mNTheta; i++ ) 
+    if( k< mTNphi[i+1] ) break;
+  
+  // valeurs limites de theta dans l'hemisphere Nord
+  tetMin= mTheta[i];
+  tetMax= mTheta[i+1];
+  // valeurs limites de theta dans l'hemisphere Sud
+  if (fgzn) {
+    tetMin= Pi-mTheta[i+1];
+    tetMax= Pi-mTheta[i];
+  }
+  
+  // pixel correspondant dans l'hemisphere Nord
+  if (fgzn) k= mTNphi[i+1]-k+mTNphi[i]-1;
+  
+  // indice j de discretisation ( phi= j*dphi )
+  j= k-mTNphi[i];
+  dphi= (r_4)DeuxPi/(r_4)mNphi[i];
+  
+  // valeurs limites de phi 
+  phiMin= dphi*(r_4)(j);
+  phiMax= dphi*(r_4)(j+1);
+  return;
+}
+
+/* --Methode-- */
+//++
+int_4 SphereThetaPhi::NbThetaSlices() const
+
+//    Retourne le nombre de tranches en theta sur la sphere
+//--
+{
+  int_4 nbslices;
+  nbslices= 2*mNTheta;
+  return(nbslices);
+}
+
+/* --Methode-- */
+//++
+int_4 SphereThetaPhi::NPhi(int_4 kt) const
+
+//    Retourne le nombre de pixels en phi de la tranche kt
+//-- 
+{
+  
+  int_4 nbpix;
+  
+  // verification
+  if( (kt< 0) || (kt>= 2*mNTheta) ) return(-1);
+  
+  // si on se trouve dans l'hemisphere Sud
+  if( kt >= mNTheta ) {
+    kt= 2*mNTheta-1-kt; 
+  }
+  
+  // nombre de pixels
+  nbpix= mNphi[kt];
+  return(nbpix);
+}
+
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::Theta(int_4 kt,r_4& tetMin,r_4& tetMax) 
+
+//    Retourne les valeurs de theta limitant la tranche kt
+//--
+{
+  
+  bool fgzn= false;
+  
+  // verification
+  if( (kt< 0) || (kt>= 2*mNTheta) ) {
+    tetMin= -99999.; 
+    tetMax= -99999.; 
+    return;
+  }
+
+  // si on se trouve dans l'hemisphere Sud
+  if( kt >= mNTheta ) {
+    fgzn= true;
+    kt= 2*mNTheta-1-kt; 
+  }
+  
+  // valeurs limites de theta dans l'hemisphere Nord
+  tetMin= mTheta[kt];
+  tetMax= mTheta[kt+1];
+  // valeurs limites de theta dans l'hemisphere Sud
+  if (fgzn) {
+    tetMin= Pi-mTheta[kt+1];
+    tetMax= Pi-mTheta[kt];
+  }  
+}
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax) 
+
+//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
+//--
+{
+  
+  r_4 dphi;
+  
+  // verification
+  if( (kt< 0) || (kt>= 2*mNTheta) ) {
+    phiMin= -99999.; 
+    phiMax= -99999.; 
+    return;
+  }
+  
+  // si on se trouve dans l'hemisphere Sud
+  if( kt >= mNTheta ) kt= 2*mNTheta-1-kt;
+  
+  // verifie si la tranche kt contient au moins jp pixels
+  if( (jp< 0) || (jp >= mNphi[kt]) ) {
+    phiMin= -88888.; 
+    phiMax= -88888.; 
+    return;
+  }
+  
+  dphi= (r_4)DeuxPi/(r_4)mNphi[kt];
+  phiMin= dphi*(r_4)(jp);
+  phiMax= dphi*(r_4)(jp+1);
+  return;
+}
+
+/* --Methode-- */
+//++
+int_4 SphereThetaPhi::Index(int_4 kt,int_4 jp) const
+
+//    Retourne l'indice du pixel d'indice jp dans la tranche kt 
+//--
+{
+  
+  int_4 k;
+  bool fgzn= false;
+  
+  // si on se trouve dans l'hemisphere Sud
+  if( kt >= mNTheta ) {
+    fgzn= true; 
+    kt= 2*mNTheta-1-kt;
+  }
+  
+  // si la tranche kt contient au moins i pixels
+  if( (jp>=0) && (jp<mNphi[kt]) ) {
+    // dans l'hemisphere Sud
+    if (fgzn) k= mNPix-mTNphi[kt+1]+jp;
+    // dans l'hemisphere Nord
+    else k= mTNphi[kt]+jp;
+  }
+  else{
+    k= 9999;
+    printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
+  }
+  return(k);
+}
+
+/* --Methode-- */
+//++
+void SphereThetaPhi::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 
+
+//    Retourne les indices kt et jp du pixel d'indice k
+//--
+{
+  
+  bool fgzn= false;
+  
+  // si on se trouve dans l'hemisphere Sud
+  if( k >= mNPix/2 ) { 
+    fgzn= true; 
+    k= mNPix-1-k; 
+  }
+  
+  // on recupere l'indice kt de la tranche qui contient le pixel k
+  int i;
+  for( i=0; i< mNTheta; i++ ) 
+    if( k< mTNphi[i+1] ) break;
+  
+  // indice  kt de tranche
+  if (fgzn) kt= 2*mNTheta-1-i;
+  else kt= i;
+  
+  // indice jp de pixel
+  if (fgzn) jp= mTNphi[i+1]-k-1;
+  else jp= k-mTNphi[i];
+  
+}
+//++
+void  SphereThetaPhi::Pixelize( int_4 m, int_4 pet) 
+
+//    effectue le découpage en pixels (m et pet ont la même signification 
+//    que pour le constructeur)
+//
+//    Chaque bande de theta sera découpée en partant de phi=0 ...
+//    L'autre hémisphère est parcourue dans le même sens en phi et de
+//    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
+//    proche de l'équateur a z>0 est celui de plus petit phi de la bande la 
+//    plus proche de l'equateur a z<0).
+//--
+{
+  int_4 ntotpix,i,j;
+  //r_8 x, omeg, omeg2pi, teti, tetm, tet, nphi, costeti;
+  r_8 omeg2pi;
+  // Decodage et controle des arguments d'appel :
+  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
+  if (m < 2) m = 2;
+  if (m > 16384) m = 16384;
+  
+  // Au moins 4 et au plus 256 pixels dans la premiere bande decoupee en phi
+  if (pet < 4) pet = 4;
+  if (pet > 256) pet = 256;
+  
+  // On memorise les arguments d'appel
+  mNTheta = m;  mNPet = pet;
+  
+  // On commence par decouper l'hemisphere z>0.
+  // Creation des vecteurs contenant :
+  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
+  mTheta = new r_4[m+1];
+  
+  // Le nombre de pixels en phi de chacune des bandes en theta
+  mNphi = new int_4[m];
+  
+  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
+  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
+  mTNphi = new int_4[m+1];
+  
+  // Calcul du nombre total de pixels dans chaque bande pour optimiser
+  // le rapport largeur/hauteur des pixels 
+  
+  //calotte polaire
+  mTNphi[0]=0;
+  mNphi[0] = 1;
+  
+  //bandes jusqu'a l'equateur
+  for(j=1; j < m; j++) {
+    mTNphi[j] = mTNphi[j-1]+mNphi[j-1];
+    mNphi[j]  = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
+  }
+  
+  // Nombre total de pixels sur l'hemisphere
+  ntotpix = mTNphi[m-1]+mNphi[m-1];
+  mTNphi[m] = ntotpix;
+  // et sur la sphere entiere
+  mNPix=2*ntotpix;
+  
+  // Creation et initialisation du vecteur des contenus des pixels 
+  //mPix = new r_8[mNPix];
+  //for(i=0; i<mNPix; i++)  mPix[i] = 0.;
+  _pixel=new PDataArray<r_8>(mNPix,true);
+  for(i=0; i<mNPix; i++)  *(_pixel->Data()+i)=0.;
+  // Determination des limites des bandes en theta :
+  // omeg est l'angle solide couvert par chaque pixel,
+  // une bande donnee kt couvre un angle solide mNphi[kt]*omeg
+  // egal a 2* Pi*(cos mTheta[kt]-cos mTheta[kt+1]). De meme, l'angle solide 
+  //de la
+  // calotte allant du pole a la limite haute de la bande kt vaut
+  // 2* Pi*(1.-cos mTheta[kt+1])= mTNphi[kt]*omeg...
+  
+  omeg2pi = 1./(r_8)ntotpix;
+  mOmeg = omeg2pi*DeuxPi;
+  
+  for(j=0; j <= m; j++) {
+    mTheta[j] = acos(1.-(double)mTNphi[j]*omeg2pi);
+  }
+}
+
+
+//++
+void  SphereThetaPhi::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const 
+
+//    Retourne, pour la tranche en theta d'indice 'index' le theta 
+//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
+//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
+//    correspondantes
+//--
+
+{
+  cout << "entree GetThetaSlice, couche no " << index << endl;
+  if (index<0 || index > NbThetaSlices()) {
+    // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
+    cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
+  THROW(rangeCheckErr);
+  }
+  int_4 iring=Index(index,0);
+  int_4 bid=this->NPhi(index);
+  int lring=bid;
+  cout << " iring= " << iring << " lring= " << lring << endl;
+  phi.Realloc(lring);
+  value.Realloc(lring);
+  float T=0.;
+  float F=0.;
+  for (int kk=0; kk<lring;kk++) {
+    PixThetaPhi(kk+iring,T,F);
+    phi(kk)=F;
+    value(kk)=PixVal(kk+iring);
+  }
+  theta=T;
+
+}
+
+
+/* --Methode-- */
+
Index: trunk/SophyaLib/Samba/spherethetaphi.h
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.h	(revision 228)
+++ trunk/SophyaLib/Samba/spherethetaphi.h	(revision 228)
@@ -0,0 +1,100 @@
+#ifndef SPHERETHETAPHI_SEEN
+#define SPHERETHETAPHI_SEEN
+
+#include "sphericalmap.h"
+#include "pdataarray.h"
+#include "cvector.h"
+
+// ***************** CLASSE SphereThetaPhi *****************************
+
+
+class SphereThetaPhi : public SphericalMap {
+public :
+
+
+SphereThetaPhi(int_4 m, int_4 pet);
+SphereThetaPhi();
+SphereThetaPhi(char* flnm);
+virtual ~SphereThetaPhi();
+
+// ------------ Persistence handling
+
+     enum {classId = 0xF001 };
+int_4                  ClassId() const        { return classId; }
+
+virtual void           WriteSelf(POutPersist&) const;
+virtual void           ReadSelf(PInPersist&);
+
+// ------------------ Definition of PixelMap abstract methods
+
+// number of pixels
+virtual  int_4           NbPixels() const;
+
+// Value of pixel number k
+virtual  r_8&       PixVal(int_4 k);
+virtual  r_8 const&       PixVal(int_4 k) const;
+
+// Index of pixel at (theta,phi)
+virtual  int_4           PixIndexSph(float theta, float phi) const;
+
+// Spherical coordinates of center of pixel number k
+virtual  void          PixThetaPhi(int_4 k, float& theta, float& phi) const;
+
+// Pixel Solid angle  (steradians)
+virtual r_8         PixSolAngle(int_4 dummy) const; 
+// --------------- Specific methods 
+
+
+/* Valeurs de theta des paralleles et phi des meridiens limitant          */
+/* le pixel d'indice k                                                   */
+virtual  void          Limits(int_4 k, float& tet1, float& tet2, 
+				float& phi1, float& phi2 );
+
+/* Nombre de tranches en theta                                           */
+int_4           NbThetaSlices() const;
+
+/* Nombre de pixels en phi de la tranche d'indice kt                     */
+int_4           NPhi(int_4 kt) const;
+
+/* Renvoie dans t1,t2 les valeurs respectives de theta min et theta max  */
+/* de la tranche d'indice kt                                             */
+void          Theta(int_4 kt, float& t1, float& t2);
+
+/* Renvoie dans p1,p2 les valeurs phimin et phimax du pixel d'indice jp  */
+/* dans la tranche d'indice kt                                           */
+void          Phi(int_4 kt, int_4 jp, float& p1, float& p2);
+
+/* Renvoie l'indice k du pixel d'indice jp dans la tranche d'indice kt   */
+int_4           Index(int_4 kt, int_4 jp) const;
+
+
+/* Indice kt de la tranche et indice jp du pixel d'indice k            */
+void          ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp);
+
+void Pixelize(int_4,int_4); 
+void GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const; 
+
+
+private :
+// ------------- méthodes internes ----------------------
+          
+
+void	      InitNul();
+void          Clear();
+
+// ------------- variables internes -----------------------
+int_4 mNTheta,mNPet;
+int_4 mNPix;
+r_4* mTheta;
+r_8  mOmeg;
+int_4* mNphi;
+int_4* mTNphi;
+//r_8* mPix;
+PDataArray<r_8>* _pixel;
+
+
+};
+
+
+
+#endif
Index: trunk/SophyaLib/Samba/sphericalmap.h
===================================================================
--- trunk/SophyaLib/Samba/sphericalmap.h	(revision 228)
+++ trunk/SophyaLib/Samba/sphericalmap.h	(revision 228)
@@ -0,0 +1,45 @@
+#ifndef SPHERICALMAP_SEEN
+#define SPHERICALMAP_SEEN
+
+// valeurs de Pi, 2*Pi, etc                                              
+#include "nbmath.h"
+#include <math.h>
+#include "pixelmap.h"
+#include "cvector.h" 
+// Map of pixels on a whole sphere.
+// Class hierarchy :
+//   PixelMap
+//      SphericalMap
+//         SphereThetaPhi
+//         SphereGorski
+//         SphereIco
+//      LocalMap
+
+class SphericalMap : public PixelMap {
+
+public :
+SphericalMap() {};
+virtual	 ~SphericalMap() {};
+
+   // Overloading of () to access pixel number k.
+   inline  r_8&        operator()(int_4 k)
+                          {return(PixVal(k));}
+   inline  r_8 const&  operator()(int_4 k) const
+                          {return(PixVal(k));}
+inline  r_8&        operator()(float theta, float phi)
+                { return(PixValSph(theta, phi)) ; };
+inline  r_8 const&        operator()(float theta, float phi) const
+                { return(PixValSph(theta, phi)) ; };
+
+
+virtual int_4 NbThetaSlices() const=0;
+virtual void GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const=0; 
+};
+
+#endif
+
+
+
+
+
+
Index: trunk/SophyaLib/Samba/unitvec.h
===================================================================
--- trunk/SophyaLib/Samba/unitvec.h	(revision 228)
+++ trunk/SophyaLib/Samba/unitvec.h	(revision 228)
@@ -0,0 +1,57 @@
+#ifndef UNITVEC_H_SEEN
+#define UNIITVEC_H_SEEN
+
+#include <math.h>
+#include <iostream.h>
+
+
+
+//  Classe pour faire des vecteurs normes x,y,z a partir de Teta,Phi
+//  et l inverse.  On peut calculer le produit scalaire et les vecteurs 
+//  perpendiculaires.
+//  Exemple:
+//  UnitVector u1(0.5, 0.7);   // Vers Teta=0.5 Phi = 0.7 Rad.
+//  UnitVector u2(1.,2.,1.);   // Vecteur norme parallele a (x=1., y=2., z=1.)
+//  u2.Teta();                 // Direction teta de u2
+//  u1.psc(u2);                // u1 prod. scalaire u2
+//  UnitVector vp1 = u2.vperp();   // Vecteur perpendiculaire a u1 , meme phi (ou a 180 deg)
+//  UnitVector vp2 = u2.vperp2();  // Vecteur perpendiculaire a u1 , a phi+90 deg
+
+//  
+class UnitVector {
+public:
+  UnitVector(double teta=0., double phi=0.) 
+    {  set(teta, phi); }
+  UnitVector(double vx, double vy, double vz) 
+      { setxyz(vx, vy, vz); }
+
+  inline void set(double teta, double phi) 
+    { t=teta; p=phi; z = cos(teta); x = sin(teta)*cos(phi);  y =  sin(teta)*sin(phi); }
+  inline void setxyz(double vx, double vy, double vz) 
+    { double r = sqrt(vx*vx+vy*vy+vz*vz); x = vx/r; y = vy/r; z = vz/r;
+    t = acos(z); p = atan2(y,x); }
+
+  inline double psc(UnitVector& u) { return(x*u.x+y*u.y+z*u.z); }
+
+  inline double Teta() { return(t); }
+  inline double Phi() { return(p); }
+  inline double X() { return(x); }
+  inline double Y() { return(y); }
+  inline double Z() { return(z); }
+
+  inline UnitVector vperp() {double tp, pp; pp = p;  
+     tp = t+0.5*M_PI; if (tp > M_PI) { tp = tp-M_PI; pp += M_PI*0.5; }
+     return(UnitVector(tp, pp)); }
+  inline UnitVector vperp2() {double tp, pp; tp = M_PI*0.5;  pp = p + 0.5*M_PI; 
+     return(UnitVector(tp, pp)); }
+
+  inline friend ostream& operator << (ostream& s, UnitVector const& u)
+      { s << "UnitVector: X,Y,Z= (" << u.x << "," << u.y << "," << u.z 
+          << ") Teta,Phi= " << u.t << "," << u.p << endl ;  return(s); }
+
+  double x,y,z;
+  double t,p;
+};
+
+
+#endif
