Index: trunk/SophyaLib/Samba/alm.cc
===================================================================
--- trunk/SophyaLib/Samba/alm.cc	(revision 1682)
+++ trunk/SophyaLib/Samba/alm.cc	(revision 1683)
@@ -78,5 +78,5 @@
   return powsp;
 }
-
+ 
 #ifdef __CXX_PRAGMA_TEMPLATES__
 #pragma define_template Alm<r_8>
Index: trunk/SophyaLib/Samba/alm.h
===================================================================
--- trunk/SophyaLib/Samba/alm.h	(revision 1682)
+++ trunk/SophyaLib/Samba/alm.h	(revision 1683)
@@ -30,4 +30,13 @@
 TVector<T> powerSpectrum() const;
 
+inline  Alm<T>& SetComplexT(complex<T> a)
+   {
+     this->SetT(a);
+     return *this;
+   }
+ inline Alm<T>& operator = (complex<T> a)
+ {
+   return SetComplexT(a);
+ } 
 
 };
Index: trunk/SophyaLib/Samba/lambdaBuilder.cc
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 1682)
+++ trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 1683)
@@ -71,4 +71,10 @@
       cth_=cos(theta);
       sth_=sin(theta);
+      array_init(lmax, mmax);
+    }
+LambdaLMBuilder::LambdaLMBuilder(r_8 costet, r_8 sintet,int_4 lmax, int_4 mmax)
+    {
+      cth_=costet;
+      sth_=sintet;
       array_init(lmax, mmax);
     }
Index: trunk/SophyaLib/Samba/lambdaBuilder.h
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 1682)
+++ trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 1683)
@@ -41,8 +41,16 @@
 LambdaLMBuilder() {}
 LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax);
+LambdaLMBuilder(r_8 ct, r_8 st,int_4 lmax, int_4 mmax);
 virtual ~LambdaLMBuilder() {};
 
  /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
 inline double lamlm(int l, int m) const { return lambda_(l,m); }
+
+inline double lamlm(int index) const { return lambda_(index); }
+
+
+  //Return pointer to first element address of the alm's
+//inline r_8* Data()  {return lambda_.Data();}
+
 
  private:
Index: trunk/SophyaLib/Samba/sphericaltransformserver.cc
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 1682)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 1683)
@@ -7,4 +7,7 @@
 #include "nbrandom.h"
 #include "nbmath.h"
+#include "timing.h"
+//#include "spherehealpix.h"
+
 
 /*! \class SOPHYA::SphericalTransformServer
@@ -184,12 +187,13 @@
 	throw IOExc(" unknown type of sphere: " + (string)sphere_type );
       }
-  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
-  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
-  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
+  //  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
+  //  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
+  //  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
   if (nlmax>3*nsmax-1) 
     {
-      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
+      //     cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
       if (strncmp(sphere_type,"TETAFI",6) == 0)
 	{
+	  cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
 	  cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
 	}
@@ -221,21 +225,8 @@
       for (m = 0; m <= nmmax; m++)
 	{
-	  //  somme sur l de m a l'infini
 	  b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
-	  //    if (ith==0 && m==0)
-	  //	{
-	  //	  cout << " guy: lmm= " <<  lb.lamlm(m,m) << " alm " << alm(m,m) << "b00= " <<  b_m_theta(m) << endl;
-	  //	} 
 	  for (int l = m+1; l<= nlmax; l++)
 	    {
 	       b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
-
-
-	       //     if (ith==0 && m==0)
-	       //	 {
-	       //	   cout << " guy:l=" << l << " m= " << m << " lmm= " <<  lb.lamlm(l,m) << " alm " << alm(l,m) << "b00= " <<  b_m_theta(m) << endl;
-
-	       //	 } 
-
 	    }
 	}
@@ -245,5 +236,4 @@
       for (m=1;m<=nmmax;m++)
 	{
-	  //compiler doesn't have conj()
 	  b_m_theta(-m) = conj(b_m_theta(m));
 	}
@@ -251,10 +241,8 @@
       //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
       // ---------------------------------------------------------------*/
-      //    TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0); 
-	  TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
+      TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
       for (int i=0;i< nph;i++)
 	{
-	  // 	  map(pixNumber(i))=Temp(i).real();
-	  	  	  map(pixNumber(i))=Temp(i);
+	  map(pixNumber(i))=Temp(i);
 	}
     }
@@ -421,5 +409,4 @@
     }
 
-  //sortie.ReSize(nph);
   TVector<T> sortie;
 
@@ -432,12 +419,61 @@
  /*! \fn  Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
 
-\return the Alm coefficients from analysis of a temperature map.
+\return the Alm coefficients from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
 
     \param<nlmax> : maximum value of the l index
 
      \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
+
+\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
   */ 
 template<class T>
- Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
+void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
+{
+  int_4  nmmax = nlmax;
+  //  PrtTim("appel  carteVersAlm");
+  carteVersAlm(map, nlmax, cos_theta_cut, alm);
+  //  PrtTim("retour  carteVersAlm");
+  if (iterationOrder > 0)
+    {
+      TVector<int_4> fact(iterationOrder+2);
+      fact(0) = 1;
+      for (int k=1; k <= iterationOrder+1; k++)
+	{
+	  fact(k) = fact(k-1)*k;
+	}
+      Alm<T> alm2(alm);
+      T Tzero = (T)0.;
+      complex<T> complexZero = complex<T>(Tzero, Tzero);
+      alm = complexZero;
+      int signe = 1;
+      int nbIteration = iterationOrder+1;
+      for (int k=1; k <= nbIteration; k++)
+	{
+	  T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
+	  for (int m = 0; m <= nmmax; m++)
+	    {
+	      for (int l = m; l<= nlmax; l++)
+		{
+		  alm(l,m) += facMult*alm2(l,m); 
+		}
+	    }
+	  if (k == nbIteration) break;
+	  signe = -signe;
+	  for (int k=0; k< map.NbPixels(); k++) map(k) = (T)0.;
+	  //        synthetize a map from the estimated alm
+	  //	  PrtTim("appel  GenerateFromAlm");
+	  GenerateFromAlm( map, map.SizeIndex(), alm2);
+	  //	  PrtTim("retour  GenerateFromAlm");
+	  alm2 = complexZero;
+	  //        analyse the new map
+	  //	  PrtTim("appel  carteVersAlm");
+	  carteVersAlm(map, nlmax, cos_theta_cut, alm2);
+	  //	  PrtTim("retour  carteVersAlm");
+	}
+    }
+}
+
+template<class T>
+ void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, Alm<T>& alm) const
 {
   
@@ -450,5 +486,5 @@
   int_4  nmmax = nlmax;
   TVector< complex<T> > phase(nmmax+1);
-  Alm<T> alm;
+  
   alm.ReSizeToLmax(nlmax);
   for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
@@ -456,43 +492,56 @@
       r_8 phi0;
       r_8 theta;
+      //  PrtTim("debut 1ere tranche ");
       map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
-      for (int i=0;i< nmmax+1;i++)
-	{
-	  phase(i)=0; 
-	}
+      phase = complex<T>((T)0.,(T)0.);
       double cth = cos(theta);
       
       //part of the sky out of the symetric cut
       bool keep_it = (fabs(cth) >= cos_theta_cut); 
+
+      //    PrtTim("fin 1ere tranche ");
+  
+      if (keep_it)
+	{
+	  //	  phase = CFromFourierAnalysis(nmmax,data,phi0);
+	  // 	  PrtTim("avant Fourier ");
+	  CFromFourierAnalysis(nmmax,data,phase, phi0);
+	  // 	  PrtTim("apres Fourier ");
+
+	}
       
-      if (keep_it)
-	{
-	  // tableau datain a supprimer
-	  //  TVector<complex<T> > datain(pixNumber.NElts());
-	  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);
-
-	  //  phase = CFromFourierAnalysis(nmmax,datain,phi0);
-	  phase = CFromFourierAnalysis(nmmax,data,phi0);
-
+//      ---------------------------------------------------------------------
+//	computes the a_lm by integrating over theta
+//	lambda_lm(theta) * phas_m(theta)
+//	for each m and l
+//	-----------------------------------------------------------------------
+      //        PrtTim("avant instanciation LM ");
+      LambdaLMBuilder lb(theta,nlmax,nmmax);
+      //        PrtTim("apres instanciation LM ");
+      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
+
+      //   PrtTim("avant mise a jour Alm ");
+      complex<T> fi;
+      T facteur;
+      int index;
+      for (int m = 0; m <= nmmax; m++)
+	{
+	  fi = phase(m);
+	  for (int l = m; l<= nlmax; l++)
+	    {
+	         index = alm.indexOfElement(l,m);
+		 //  facteur = (T)(lb.lamlm(l,m) * domega);
+		    facteur = (T)(lb.lamlm(index) * domega);
+		  // alm(l,m) += facteur * fi ; 
+		      alm(index) += facteur * fi ; 
+	    }
 	}
       
-      /*-----------------------------------------------------------------------
-	computes the a_lm by integrating over theta
-	lambda_lm(theta) * phas_m(theta)
-	for each m and l
-	-----------------------------------------------------------------------*/
-      //      LambdaBuilder lb(theta,nlmax,nmmax);
-      LambdaLMBuilder lb(theta,nlmax,nmmax);
-      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
-      for (int m = 0; m <= nmmax; m++)
-	{
-	  alm(m,m) += (T)lb.lamlm(m,m) * phase(m) * (T)domega; //m,m even
-	  for (int l = m+1; l<= nlmax; l++)
-	    {
-	      alm(l,m) += (T)lb.lamlm(l,m) * phase(m)*(T)domega; 
-	    }
-	}
-    }
-  return alm;
+
+      
+      //
+      //
+      //       PrtTim("apres mise a jour Alm ");
+    }
 }
   /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
@@ -518,16 +567,13 @@
   fftIntfPtr_-> FFTForward(datain, transformedData);
 
-  //dataout.ReSize(nmmax+1);
   TVector< complex<T> > dataout(nmmax+1);
 
   int im_max=min(nph,nmmax+1);
   int i;
-  for (i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
+  dataout = complex<T>((T)0.,(T)0.);
+  //  for (i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
   for (i=0;i<im_max;i++) dataout(i)=transformedData(i);
 
 
-  //  for (int i = 0;i <im_max;i++){
-  //    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
-  //  }
   for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
   for (i = 0;i <dataout.NElts();i++){
@@ -542,5 +588,5 @@
 same as previous one, but with a "datain" which is real (not complex) */
 template<class T>
-TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const
+void SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, TVector< complex<T> >& dataout, r_8 phi0) const
 {
   //=======================================================================
@@ -559,44 +605,47 @@
       throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
     }
+  if (nph%2 != 0 )
+    {
+      throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
+    }
   TVector<complex<T> > transformedData;
-    // a remodifier
-  //FFTPackServer ffts;
-  //ffts.setNormalize(false);
-  //ffts.FFTForward(datain, transformedData);
-
+
+  // la taille du vecteur complexe retourne est nph/2+1 (si la taille 
+  // du vecteur reel entre est nph)
   fftIntfPtr_-> FFTForward(datain, transformedData);
-    //
-
-  //dataout.ReSize(nmmax+1);
-  TVector< complex<T> > dataout(nmmax+1);
+
+  //  TVector< complex<T> > dataout(nmmax+1);
+  dataout.ReSize(nmmax+1);
 
   // on transfere le resultat de la fft dans dataout.
-  // on s'assure que ca ne depasse pas la taille de dataout
-  int sizeOfTransformToGet = min(transformedData.NElts(),nmmax+1);
-  //  int im_max=min(transformedData.NElts()-1,nmmax);
+
+  int maxFreqAccessiblesParFFT = min(nph/2,nmmax);
   int i;
-  for (i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);
-
-
-  // si dataout n'est pas plein, on complete jusqu'a  nph valeurs (a moins 
+  for (i=0;i<=maxFreqAccessiblesParFFT;i++) dataout(i)=transformedData(i);
+
+
+  // si dataout n'est pas plein, on complete jusqu'a  nph+1 valeurs (a moins 
   // que dataout ne soit plein avant d'atteindre nph)
-  if (sizeOfTransformToGet == (transformedData.NElts()))
-    {
-      for (i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
-	{
-
-	  //	  dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
+  if (maxFreqAccessiblesParFFT != nmmax )
+    {
+      int maxMfft = min(nph,nmmax);
+      for (i=maxFreqAccessiblesParFFT+1; i<=maxMfft; i++)
+	{
 	  dataout(i) = conj(dataout(nph-i) );
 	}
       // on conplete, si necessaire, par periodicite
-      for (int kk=nph; kk<dataout.NElts(); kk++) 
-	{
-	  dataout(kk)=dataout(kk%nph);
-	}
-    }
-  for (i = 0;i <dataout.NElts();i++){
-    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
-  }
-  return dataout;
+      if ( maxMfft != nmmax )
+	{
+	  for (int kk=nph+1; kk <= nmmax; kk++) 
+	    {
+	      dataout(kk)=dataout(kk%nph);
+	    }
+	}
+    }
+  for (i = 0;i <dataout.NElts();i++)
+    {
+      dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
+    }
+  //  return dataout;
 }
 
@@ -701,5 +750,60 @@
  */ 
 template<class T>
-void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
+void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& mapq,
+					      SphericalMap<T>& mapu,
+					      Alm<T>& alme,
+					      Alm<T>& almb,
+					      int_4 nlmax,
+					      r_8 cos_theta_cut, 
+					      int iterationOrder) const
+{
+  int_4  nmmax = nlmax;
+  carteVersAlm(mapq, mapu, alme, almb, nlmax, cos_theta_cut);
+  if (iterationOrder > 0)
+    {
+      TVector<int_4> fact(iterationOrder+2);
+      fact(0) = 1;
+      for (int k=1; k <= iterationOrder+1; k++)
+	{
+	  fact(k) = fact(k-1)*k;
+	}
+      Alm<T> alme2(alme);
+      Alm<T> almb2(almb);
+      T Tzero = (T)0.;
+      complex<T> complexZero = complex<T>(Tzero, Tzero);
+      alme = complexZero;
+      almb = complexZero;
+      int signe = 1;
+      int nbIteration = iterationOrder+1;
+      for (int k=1; k <= nbIteration; k++)
+	{
+	  T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
+	  for (int m = 0; m <= nmmax; m++)
+	    {
+	      for (int l = m; l<= nlmax; l++)
+		{
+		  alme(l,m) += facMult*alme2(l,m); 
+		  almb(l,m) += facMult*almb2(l,m); 
+		}
+	    }
+	  if (k == nbIteration) break;
+	  signe = -signe;
+	  for (int k=0; k< mapq.NbPixels(); k++)
+	    {
+	      mapq(k) = (T)0.;
+	      mapu(k) = (T)0.;
+	    }
+	  //        synthetize a map from the estimated alm
+	  GenerateFromAlm(mapq,mapu,mapq.SizeIndex(),alme2,almb2); 
+	  alme2 = complexZero;
+	  almb2 = complexZero;
+	  //        analyse the new map
+	  carteVersAlm(mapq, mapu, alme2, almb2, nlmax, cos_theta_cut);
+	}
+    }
+}
+
+template<class T>
+void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& mapq,
 					      const SphericalMap<T>& mapu,
 					      Alm<T>& alme,
@@ -797,11 +901,9 @@
   //  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
 
-  // phaseq = CFromFourierAnalysis(nmmax,datain,phi0); 
-  phaseq = CFromFourierAnalysis(nmmax,dataq,phi0); 
-
-  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(datau(kk),0.);
-
-  // phaseu=  CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
-  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 
+  //  phaseq = CFromFourierAnalysis(nmmax,dataq,phi0); 
+  CFromFourierAnalysis(nmmax,dataq,phaseq, phi0); 
+
+  //  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 
+  CFromFourierAnalysis(nmmax,datau,phaseu, phi0); 
 
   LambdaWXBuilder lwxb(theta,nlmax,nmmax);
@@ -1141,14 +1243,20 @@
 /*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
 
-\return power spectrum from analysis of a temperature map.
+\return power spectrum from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
 
      \param<nlmax> : maximum value of the l index
 
      \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
+
+\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
+
   */ 
 template <class T>
-TVector<T>  SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
-{
-    Alm<T> alm=DecomposeToAlm( sph, nlmax, cos_theta_cut);
+TVector<T>  SphericalTransformServer<T>::DecomposeToCl(SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
+{
+  Alm<T> alm;
+  DecomposeToAlm( sph, alm, nlmax, cos_theta_cut, iterationOrder);
+    //    cout << " SphericalTransformServer : impression alm " << endl;
+    //     alm.Print();
   // power spectrum
      return  alm.powerSpectrum();
Index: trunk/SophyaLib/Samba/sphericaltransformserver.h
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 1682)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 1683)
@@ -46,26 +46,43 @@
 
 
-  Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const;
+void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
 
-  void DecomposeToAlm(const SphericalMap<T>& mapq, const SphericalMap<T>& mapu,
-		     Alm<T>& a2lme, Alm<T>& a2lmb,
-		     int_4 nlmax, r_8 cos_theta_cut) const;
+void DecomposeToAlm(SphericalMap<T>& mapq,
+		    SphericalMap<T>& mapu,
+		    Alm<T>& alme,
+		    Alm<T>& almb,
+		    int_4 nlmax,
+		    r_8 cos_theta_cut, 
+		    int iterationOrder = 0) const;
 
-  TVector<T>  DecomposeToCl(const SphericalMap<T>& sph,  
-			   int_4 nlmax, r_8 cos_theta_cut) const;
+TVector<T>  DecomposeToCl(SphericalMap<T>& sph,  
+			   int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
 
 
  private:
-  TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+
+void carteVersAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, Alm<T>& alm) const;
+
+void carteVersAlm(const SphericalMap<T>& mapq,
+		  const SphericalMap<T>& mapu,
+		  Alm<T>& alme,
+		  Alm<T>& almb,
+		  int_4 nlmax,
+		  r_8 cos_theta_cut) const;
+
+
+TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
 				  int_4 nph, r_8 phi0) const;
-  TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
 				  int_4 nph, r_8 phi0) const;
-
-  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
+TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
 				  const TVector<complex<T> > datain,
 				  r_8 phi0) const;
-  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
-			 const TVector<T> datain,  
-				  r_8 phi0) const;
+  //  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
+  //			 const TVector<T> datain,  
+  //				  r_8 phi0) const;
+  void CFromFourierAnalysis(int_4 mmax, const TVector<T> datain,
+			    TVector< complex<T> >& dataout,
+			    r_8 phi0) const;
 
   void almFromWX(int_4 nlmax, int_4 nmmax, r_8 phi0, 
@@ -90,4 +107,5 @@
 
  FFTServerInterface* fftIntfPtr_;
+
 };
 } // Fin du namespace
Index: trunk/SophyaLib/TArray/triangmtx.h
===================================================================
--- trunk/SophyaLib/TArray/triangmtx.h	(revision 1682)
+++ trunk/SophyaLib/TArray/triangmtx.h	(revision 1683)
@@ -15,5 +15,5 @@
 
 namespace SOPHYA {
-
+  
 //! Class for inferior triangular matrix (base class for the class Alm)
 template <class T>
@@ -35,38 +35,66 @@
   }
 
-//! Equal operator
-inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
-{
-  elem_=a.elem_;
-  long_diag_ = a.long_diag_;
-  return *this;
-}
+TriangularMatrix<T>& SetT(T a)
+  {
+    if (long_diag_ < 1)
+    throw RangeCheckError("TriangularMatrix<T>::SetT(T )  - TriangularMatrix not dimensionned ! ");
+    elem_ = a;
+    return (*this);
+  }
 
 //! () operator : access to elements row \b l and column \b m
 inline T& operator()(int l, int m) 
   {
-      return  elem_(adr_ij(l,m));
+      return  elem_(indexOfElement(l,m));
   }
+
+inline T& operator()(int index) 
+  {
+      return  elem_(index);
+  }
+
+
 //! () operator : access to elements row \b l and column \b m
 inline T const& operator()(int l, int m) const 
   {
-      return *(elem_.Begin()+ adr_ij(l,m));
+      return *(elem_.Begin()+ indexOfElement(l,m));
   }
+
+inline T const& operator()(int index) const 
+  {
+      return *(elem_.Begin()+ index);
+  }
+
 
 //! Return number of rows
 inline  int_4  rowNumber() const {return (int_4)long_diag_;}
 
+void Print(int nbLignes=0)
+  {
+    if (nbLignes == 0 ) nbLignes = long_diag_;
+    cout << " ***** matrice triangulaire : ********* " << endl;
+    for (int k=0; k < nbLignes; k++)
+      {
+	for (int kc = 0; kc <= k ; kc++)
+	  {
+	    cout << " " << elem_(indexOfElement(k,kc));
+	  }
+	cout << endl;
+      }
+    cout << "---------------- fin matrice ------------" << endl;
+  }
+
+  //Return pointer to first element address
+  //inline T* Data()  {return elem_.Begin();}
+
+//! compute the address of an element in the single array representing the matrix
+inline uint_4 indexOfElement(int i,int j) const 
+{
+  //  return(i*(i+1)/2+j);
+  // the (inferior triangular )matrix is stored column by column
+  return(i+ long_diag_*j-j*(j+1)/2);
+}
+
 private: 
-//! compute the address of an element in the single array representing the matrix
-inline uint_4 adr_ij(int i,int j) const 
-{
-  //  int adr= i*(i+1)/2+j;
-  //  if ( adr >= elem_.Size() || adr <0 )
-  //{
-  // cout << " attention depassement dans triangularMatrix " << endl;
-  // cout << " l= " << i << " m= " << j << " tableau reserve longueur " << elem_.Size() << endl; 
-  //}
-  return(i*(i+1)/2+j);
-}
 
 uint_4 long_diag_;    //!< size of the square matrix
@@ -74,5 +102,5 @@
 
   };
-
+  
 }   // namespace SOPHYA
 
Index: trunk/SophyaProg/PrgMap/map2cl.cc
===================================================================
--- trunk/SophyaProg/PrgMap/map2cl.cc	(revision 1682)
+++ trunk/SophyaProg/PrgMap/map2cl.cc	(revision 1683)
@@ -50,5 +50,5 @@
   else {
     cout << " map2cl : Spherical harmonics analysis - HEALPix map -> Power spectrum C_l \n" 
-	 << " Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg] \n" 
+	 << " Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg] [-iter_order lval]\n" 
 	 << "        [-fitsin] [-fitsout] InFileName OutFileName \n" 
 	 << "   -float (-r4): single precision C_l and map (default = double) \n"
@@ -56,4 +56,5 @@
 	 << "   -thetacut dtdeg : Symmetric delta-theta cut (in degree) along equator \n" 
 	 << "                    (default=0 -> no cut)\n"
+         << "   -iter_order lval : 1,2,3,4... order of an iterative analysis , 3rd order is usually optimal (default=0 -> standard analysis)\n"
 	 << "   -fitsout: Select the FITS format for the output map (default PPF format) \n"
 	 << "   -fitsin : Select the FITS format for the input vector (default PPF format) \n"
@@ -69,9 +70,9 @@
 public :
 static void ComputeCl(string & infile, string & outfile, int lmax, double tcut, 
-		      bool fgfitsin, bool fgfitsout)
+		      int iterationOrder, bool fgfitsin, bool fgfitsout)
 {
   double deg2rad =  M_PI/180.;
   double minute2rad =  M_PI/(180.*60.);
-
+ 
   SphereHEALPix<T> sph;
   if (fgfitsin) {
@@ -84,5 +85,5 @@
     PInPersist ppi(infile);
     ppi >> sph;
-  }
+  } 
 
   cout << " Input map : NbPixels= " <<  sph.NbPixels() << " NSide= " 
@@ -92,8 +93,8 @@
   double ctcut = (tcut < 1.e-19) ? 0. : cos((90.-tcut)*deg2rad);
   cout << "--- Calling  DecomposeToCl() (lmax= " << lmax 
-       << " cos_theta_cut= " << ctcut << ") theta_cut=" << tcut << " deg" << endl;
+       << " cos_theta_cut= " << ctcut << " iter_order= " << iterationOrder << ") theta_cut=" << tcut << " deg"  << endl;
   // Decomposition de la carte en C_l
   SphericalTransformServer<T> sphtr;
-  TVector<T> clvec = sphtr.DecomposeToCl(sph, lmax, ctcut);
+  TVector<T> clvec = sphtr.DecomposeToCl(sph, lmax, ctcut, iterationOrder);
  
   T min, max;
@@ -125,4 +126,5 @@
   
   int lmax = 255;
+  int iterationOrder = 0;
   double tcut = 0.;
   string infile;
@@ -142,4 +144,8 @@
       if (k == narg-1) Usage(true);  // -thetacut est suivi d'un argument 
       tcut = atof(arg[k+1]);  k++;       // k++ pour sauter au suivant
+    }
+    else if (strcmp(arg[k], "-iter_order") == 0)  {
+      if (k == narg-1) Usage(true);  // -iter_order est suivi d'un argument 
+      iterationOrder = atof(arg[k+1]);  k++;   // k++ pour sauter au suivant
     }
     else if (strcmp(arg[k], "-fitsin") == 0) {
@@ -166,9 +172,9 @@
     if (fgr4) {
       cout << " SphereHEALPix<r_4>  --> Power spectrum C_l<r_4> (float)" << endl;
-      _Map2Cl<r_4>::ComputeCl(infile, outfile, lmax, tcut, fgfitsin, fgfitsout);
+      _Map2Cl<r_4>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder, fgfitsin, fgfitsout);
     }
     else {
       cout << " SphereHEALPix<r_8>  --> Power spectrum C_l<r_8> (double)" << endl;
-      _Map2Cl<r_8>::ComputeCl(infile, outfile, lmax, tcut, fgfitsin, fgfitsout);
+      _Map2Cl<r_8>::ComputeCl(infile, outfile, lmax, tcut, iterationOrder,fgfitsin, fgfitsout);
     }
   }
