Index: trunk/SophyaLib/Samba/alm.h
===================================================================
--- trunk/SophyaLib/Samba/alm.h	(revision 1715)
+++ trunk/SophyaLib/Samba/alm.h	(revision 1756)
@@ -19,5 +19,5 @@
     public :
  
-Alm() : TriangularMatrix<complex<T> >()  {};
+Alm() : TriangularMatrix<complex<T> >()  {;};
     /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
 Alm(int nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
Index: trunk/SophyaLib/Samba/lambdaBuilder.cc
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 1715)
+++ trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 1756)
@@ -40,6 +40,6 @@
 }
 
-TriangularMatrix<r_8>* LambdaLMBuilder::a_recurrence_ = NULL;
-TriangularMatrix<r_8>* LambdaLMBuilder::lam_fact_     = NULL;
+TriangularMatrix<r_8> LambdaLMBuilder::a_recurrence_ = TriangularMatrix<r_8>();
+TriangularMatrix<r_8> LambdaLMBuilder::lam_fact_     = TriangularMatrix<r_8>();
 TVector<r_8>* LambdaLMBuilder::normal_l_     = NULL;
 
@@ -81,14 +81,14 @@
 void LambdaLMBuilder::array_init(int lmax, int mmax)
    {
-     if (a_recurrence_ == NULL)
-       {
-	 a_recurrence_ = new TriangularMatrix<r_8>;
+     if (a_recurrence_.Size() == 0)
+       {
+	 //	 a_recurrence_ = new TriangularMatrix<r_8>;
 	 updateArrayRecurrence(lmax);
        }
      else
-       if ( lmax > (*a_recurrence_).rowNumber()-1     )
+       if ( lmax > (a_recurrence_).rowNumber()-1     )
 	 {
 	   cout << " WARNING : The classes LambdaXXBuilder will be more efficient if instanciated with parameter lmax = maximum value of l index which will be needed in the whole application (arrays not recomputed) " << endl;
-	   cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (*a_recurrence_).rowNumber() << endl;
+	   cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (a_recurrence_).rowNumber() << endl;
 	 updateArrayRecurrence(lmax);
 	 }
@@ -110,5 +110,5 @@
 	 r_8 lam_1=1. /bignorm2 ;
 	 //	 r_8 a_rec = LWK->a_recurr(m,m);
-	 r_8 a_rec = (*a_recurrence_)(m,m);
+	 r_8 a_rec = a_recurrence_(m,m);
 	 r_8 b_rec = 0.;
 	 for (int l=m+1; l<=lmax_; l++)
@@ -118,5 +118,5 @@
 	   b_rec=1./a_rec;
 	   //	   a_rec= LWK->a_recurr(l,m);
-	   a_rec= (*a_recurrence_)(l,m);
+	   a_rec= a_recurrence_(l,m);
 	   lam_0 = lam_1;
 	   lam_1 = lam_2;
@@ -135,12 +135,12 @@
 void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
    {
-     (*a_recurrence_).ReSizeRow(lmax+1);
+     a_recurrence_.ReSizeRow(lmax+1);
      for (int m=0; m<=lmax;m++)
      {
-       (*a_recurrence_)(m,m) = sqrt( 2.*m +3.);
+       a_recurrence_(m,m) = sqrt( 2.*m +3.);
        for (int l=m+1; l<=lmax; l++)
 	 {
 	   r_8 fl2 = (l+1.)*(l+1.);
-	   (*a_recurrence_)(l,m)=sqrt( (4.*fl2-1.)/(fl2-m*m) );
+	   a_recurrence_(l,m)=sqrt( (4.*fl2-1.)/(fl2-m*m) );
 	 }
      }
@@ -153,10 +153,10 @@
 void  LambdaLMBuilder::updateArrayLamNorm()
      {
-       (*lam_fact_).ReSizeRow(lmax_+1);
+       lam_fact_.ReSizeRow(lmax_+1);
        for(int m = 0;m<= lmax_; m++)
 	 {
 	   for (int l=m; l<=lmax_; l++)
 	     {
-	       (*lam_fact_)(l,m) =2.*(r_8)sqrt( (2.*l+1)*(l+m)*(l-m)/(2.*l-1)  );
+	       lam_fact_(l,m) =2.*(r_8)sqrt( (2.*l+1)*(l+m)*(l-m)/(2.*l-1)  );
 	     }
 	 }
@@ -214,12 +214,12 @@
 void LambdaWXBuilder::array_init()
    {
-     if (lam_fact_ == NULL || normal_l_ == NULL)
-       {
-	 lam_fact_ = new  TriangularMatrix<r_8>;
+     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
+       {
+	 //	 lam_fact_ = new  TriangularMatrix<r_8>;
 	 normal_l_ = new TVector<r_8>;
 	 updateArrayLamNorm();
        }
      else
-       if ( lmax_ > (*lam_fact_).rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
+       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
 	 {
 	   updateArrayLamNorm();
@@ -247,5 +247,5 @@
 	  r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
 	  r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
-	  r_8 lam_fact_l_m = (*lam_fact_)(l,m);
+	  r_8 lam_fact_l_m = lam_fact_(l,m);
 	  r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
 	  r_8  b_w =  c_on_s2 * lam_fact_l_m;
@@ -284,12 +284,12 @@
 void LambdaPMBuilder::array_init()
    {
-     if (lam_fact_ == NULL || normal_l_ == NULL)
-       {
-	 lam_fact_ = new  TriangularMatrix<r_8>;
+     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
+       {
+	 //	 lam_fact_ = new  TriangularMatrix<r_8>;
 	 normal_l_ = new TVector<r_8>;
 	 updateArrayLamNorm();
        }
      else
-       if ( lmax_ > (*lam_fact_).rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
+       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
 	 {
 	   updateArrayLamNorm();
@@ -318,5 +318,5 @@
 	     r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
 	     r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
-	     r_8 lam_fact_l_m = (*lam_fact_)(l,m);
+	     r_8 lam_fact_l_m = lam_fact_(l,m);
 	     r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
 	     r_8 f_w =  lam_fact_l_m/(sth_*sth_);
Index: trunk/SophyaLib/Samba/lambdaBuilder.h
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 1715)
+++ trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 1756)
@@ -60,5 +60,6 @@
 
 
- static TriangularMatrix<r_8>* a_recurrence_;
+ // static TriangularMatrix<r_8>* a_recurrence_;
+ static TriangularMatrix<r_8> a_recurrence_;
  TriangularMatrix<r_8> lambda_;
 
@@ -67,5 +68,6 @@
  void updateArrayLamNorm();
 
- static  TriangularMatrix<r_8>* lam_fact_;
+ // static  TriangularMatrix<r_8>* lam_fact_;
+ static  TriangularMatrix<r_8> lam_fact_;
  static  TVector<r_8>*  normal_l_;
  int_4 lmax_;
Index: trunk/SophyaLib/Samba/sphericaltransformserver.cc
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 1715)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 1756)
@@ -154,5 +154,5 @@
 {
   /*=======================================================================
-    computes a map form its alm for the HEALPIX pixelisation
+    computes a map from its alm for the HEALPIX pixelisation
     map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
     = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
@@ -167,11 +167,15 @@
   int_4 nmmax=nlmax;
   int_4 nsmax=0;
+  // le Resize est suppose mettre a zero
   map.Resize(pixelSizeIndex);
   char* sphere_type=map.TypeOfMap();
+      int premiereTranche = 0;
+      int derniereTranche = map.NbThetaSlices()-1;
   if (strncmp(sphere_type,"RING",4) == 0)
-    {
-      nsmax=map.SizeIndex();
-    }
+  {
+    nsmax=map.SizeIndex();
+  }
   else
+  {
     // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
     // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
@@ -179,24 +183,27 @@
     // c'est approximatif ; a raffiner.
     if (strncmp(sphere_type,"TETAFI",6) == 0)
-      {
-	nsmax=(int_4)sqrt(map.NbPixels()/12.);
-      }
+    {
+      nsmax=(int_4)sqrt(map.NbPixels()/12.);
+      premiereTranche++;
+      derniereTranche--;
+    }
     else
-      {
-	cout << " unknown type of sphere : " << sphere_type << endl;
-	throw IOExc(" unknown type of sphere: " + (string)sphere_type );
-      }
+    {
+      cout << " unknown type of sphere : " << sphere_type << endl;
+      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;
-  if (nlmax>3*nsmax-1) 
-    {
+  //  if (nlmax>3*nsmax-1) 
+  //  {
       //     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;
-	}
-    }
+  //   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;
+  //	}
+  //}
+  }
   Bm<complex<T> > b_m_theta(nmmax);
 
@@ -205,5 +212,5 @@
 
   // pour chaque tranche en theta
-  for (int_4 ith = 0; ith < map.NbThetaSlices();ith++)
+  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)
     {
       int_4 nph;
@@ -231,5 +238,4 @@
 	    }
 	}
-
     //        obtains the negative m of b(m,theta) (= complex conjugate)
 
@@ -377,5 +383,4 @@
   // si phi0 n'est pas nul, il y a juste un decalage a faire.
   //**********************************************************************
-
   TVector< complex<T> > bw(nph);
   TVector< complex<T> > dataout(nph);
@@ -410,5 +415,4 @@
 
   TVector<T> sortie;
-
   fftIntfPtr_-> FFTBackward(data, sortie);
   
@@ -419,5 +423,5 @@
  /*! \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. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
+\return the Alm coefficients from analysis of a temperature map. 
 
     \param<nlmax> : maximum value of the l index
@@ -425,6 +429,21 @@
      \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>
+void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut) const
+{
+  DecomposeToAlm(const_cast< SphericalMap<T>& >(map), alm, nlmax, cos_theta_cut, 0);
+}
+//*******************************************
+
+ /*! \fn  Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
+
+\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). If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !  */ 
 template<class T>
 void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
@@ -606,14 +625,15 @@
       throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
     }
-  if (nph%2 != 0 )
-    {
-      throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
-    }
+  // if (nph%2 != 0 )
+  //  {
+      //  throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
+  //  }
   TVector<complex<T> > transformedData;
 
   // la taille du vecteur complexe retourne est nph/2+1 (si la taille 
   // du vecteur reel entre est nph)
+  //   cout << " longueur de datain  = " << nph << endl;
   fftIntfPtr_-> FFTForward(datain, transformedData);
-
+  //  cout <<  " taille de la transformee " << transformedData.Size() << endl;
   //  TVector< complex<T> > dataout(nmmax+1);
   dataout.ReSize(nmmax+1);
@@ -732,6 +752,4 @@
     // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
 }
-
-
  /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq,
 					      const SphericalMap<T>& mapu,
@@ -749,4 +767,39 @@
 
  \c 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.
+
+
+ */ 
+template<class T>
+void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
+					      const SphericalMap<T>& mapu,
+					      Alm<T>& alme,
+					      Alm<T>& almb,
+					      int_4 nlmax,
+					      r_8 cos_theta_cut) const
+{
+  DecomposeToAlm(const_cast< SphericalMap<T>& >(mapq), const_cast< SphericalMap<T>& >(mapu), alme, almb, nlmax, cos_theta_cut);
+}
+
+ /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq,
+					      const SphericalMap<T>& mapu,
+					      Alm<T>& alme,
+					      Alm<T>& almb,
+					      int_4 nlmax,
+					      r_8 cos_theta_cut,
+					      int iterationOrder) const
+
+analysis of a polarization map into Alm coefficients.
+
+ The spheres \c mapq and \c mapu contain respectively the Stokes parameters. 
+
+ \c a2lme and \c a2lmb will receive respectively electric and magnetic Alm's
+    nlmax : maximum value of the l index
+
+ \c 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). If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !
+
+THE INPUT MAPS CAN BE MODIFIED (only if iterationOrder >0)
+
  */ 
 template<class T>
@@ -1243,5 +1296,5 @@
 
 
-/*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
+/*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
 
 \return power spectrum from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
@@ -1251,5 +1304,5 @@
      \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)
+\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !
 
   */ 
@@ -1259,6 +1312,26 @@
   Alm<T> alm;
   DecomposeToAlm( sph, alm, nlmax, cos_theta_cut, iterationOrder);
-    //    cout << " SphericalTransformServer : impression alm " << endl;
-    //     alm.Print();
+  // power spectrum
+     return  alm.powerSpectrum();
+}
+
+
+/*! \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. 
+
+     \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.
+
+
+  */ 
+
+
+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, alm, nlmax, cos_theta_cut);
   // power spectrum
      return  alm.powerSpectrum();
Index: trunk/SophyaLib/Samba/sphericaltransformserver.h
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 1715)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 1756)
@@ -46,5 +46,6 @@
 
 
-void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
+void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const;
+void DecomposeToAlm(const SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut) const;
 
 void DecomposeToAlm(SphericalMap<T>& mapq,
@@ -54,8 +55,17 @@
 		    int_4 nlmax,
 		    r_8 cos_theta_cut, 
-		    int iterationOrder = 0) const;
+		    int iterationOrder) const;
+void DecomposeToAlm(const SphericalMap<T>& mapq,
+		    const SphericalMap<T>& mapu,
+		    Alm<T>& alme,
+		    Alm<T>& almb,
+		    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;
+			   int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const;
+
+TVector<T>  DecomposeToCl(const SphericalMap<T>& sph,  
+			   int_4 nlmax, r_8 cos_theta_cut) const;
 
 
