Index: trunk/SophyaLib/Samba/alm.cc
===================================================================
--- trunk/SophyaLib/Samba/alm.cc	(revision 729)
+++ trunk/SophyaLib/Samba/alm.cc	(revision 729)
@@ -0,0 +1,88 @@
+#include "alm.h"
+template <class T>
+Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm) 
+
+{
+
+
+  /*=======================================================================
+     creates the a_lm from the power spectrum, 
+     assuming they are gaussian  and complex 
+     with a variance given by C(l)
+
+     the input file should contain : l and C(l) with *consecutive* l's
+     (missing C(l) are put to 0.)
+
+     because the map is real we have : a_l-m = (-)^m conjug(a_lm)
+     so we actually compute them only for m >= 0
+
+=======================================================================*/
+  int_4 nlmax= clin.NElts()-1;
+
+  //alm.ReSize(nlmax);
+  ReSizeRow(nlmax+1);
+
+  r_8 sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;
+  int_4 n_l = nlmax+1;
+
+
+  //    --- smoothes the initial power spectrum ---
+  TVector<T> cl=clin;
+  for (int l=0;l<n_l;l++)
+    {
+      r_8 gauss=exp(-l*(l+1.)*sig_smooth*sig_smooth);
+      cl(l)*=(T)gauss;
+    }
+  
+
+  //     --- generates randomly the alm according to their power spectrum ---
+  r_8 hsqrt2 = 1.0 / Rac2;
+
+  for (int l=0;l<n_l;l++)
+    {
+      T rms=sqrt(cl(l));
+      //        ------ m = 0 ------
+      complex<T> zeta1(NorRand());
+      (*this)(l,0)   = zeta1 * rms;
+
+      //------ m > 0 ------
+      for (int m=1;m<=l;m++)
+	{
+	  complex<T> aux1(hsqrt2);
+	  complex<T> aux2(NorRand(),NorRand());
+	  zeta1=aux1*aux2;
+	  (*this)(l,m)=rms*zeta1;
+	}
+    }
+}
+
+
+template <class T>
+TVector<T> Alm<T>::powerSpectrum() const
+{
+  int_4 nlmax=Lmax();
+
+  TVector<T>  pow(nlmax+1);
+  pow.SetTemp(true);
+
+  for (int l=0; l<=nlmax;l++)
+    {
+      pow(l)=( (*this)(l,0) ).real()*( (*this)(l,0) ).real();
+      for (int m=1; m<=l; m++)
+	{
+	  pow(l)+=2.*( (*this)(l,m) ).real()*( (*this)(l,m) ).real()+
+	    2.*( (*this)(l,m) ).imag()*( (*this)(l,m) ).imag();
+	}
+      pow(l)/=(2.*l+1.);
+    }
+  return pow;
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template Alm<r_8>
+#pragma define_template Alm<r_4>
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class Alm<r_8>;
+template class Alm<r_4>;
+#endif
Index: trunk/SophyaLib/Samba/alm.h
===================================================================
--- trunk/SophyaLib/Samba/alm.h	(revision 729)
+++ trunk/SophyaLib/Samba/alm.h	(revision 729)
@@ -0,0 +1,127 @@
+#ifndef ALM_SEEN
+#define ALM_SEEN
+#include "nbrandom.h"
+#include "nbmath.h"
+#include "tvector.h"
+#include "pexceptions.h"
+/*! Class for inferior triangular matrix (base class for the class Alm) */
+template <class T>
+class TriangularMatrix 
+  {
+    
+  public :
+
+TriangularMatrix()   {};
+  /* instanciate a triangular matrix from the number of rows */
+TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) {elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); };
+TriangularMatrix(const TriangularMatrix<T>& a,  bool share=false)  : elem_(a.elem_, share),  long_diag_(a.long_diag_) {;}
+/*! resize the matrix with a new number of rows */
+inline void ReSizeRow(int rowSize) 
+  {
+    long_diag_=(uint_4)rowSize;
+    elem_.ReSize(long_diag_*(long_diag_+1)/2);
+  }
+ inline void SetTemp(bool temp=false) const {elem_.SetTemp(temp);}
+
+inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
+{
+  elem_=a.elem_;
+  long_diag_ = a.long_diag_;
+  return *this;
+}
+inline T& operator()(int l, int m) 
+  {
+      return  elem_(adr_ij(l,m));
+  }
+inline T const& operator()(int l, int m) const 
+  {
+      return *(elem_.Begin()+ adr_ij(l,m));
+  }
+ inline  int_4  rowNumber() const {return (int_4)long_diag_;}
+  private: 
+ /*! compute the address of an element in the single array representing the matrix */
+inline uint_4 adr_ij(int i,int j) const 
+{
+  int adr= i*(i+1)/2+j;
+  if ( adr >= elem_.Size() || adr <0 )
+{
+ cout << " attention depassement dans triangularMatrix " << endl;
+ cout << " l= " << i << " m= " << j << " tableau reserve longueur " << elem_.Size() << endl; 
+}
+  return(i*(i+1)/2+j);
+}
+
+ uint_4 long_diag_;
+ NDataBlock<T> elem_;
+
+  };
+
+
+/*! class for the coefficients \f$a_{lm}\f$ of the development of a function defined on a sphere, in spherical harmonics */
+template <class T>
+class Alm : public TriangularMatrix<complex<T> >
+  {
+    public :
+ 
+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) {;}
+Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {;}
+
+Alm(const TVector<T>& clin, const r_8 fwhm) ;
+/*! resize with a new lmax */
+inline void ReSizeToLmax(int_4 nlmax) {ReSizeRow(nlmax+1);}
+inline int_4 Lmax() const {return rowNumber()-1;}
+TVector<T> powerSpectrum() const;
+
+
+};
+/*! class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ (then the size of the vector will be actually \f$2m_{max}+1)\f$ */
+template <class T>
+class Bm
+  {
+    public :
+Bm(): nmmax_(0) {;};
+    /* instanciate from the maximum absolute value of the index */
+Bm(int mmax) : nmmax_((uint_4)mmax)   { bm_.ReSize( (uint_4)(2*mmax+1) );}
+ Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;}
+ /*! resize with a new value of mmax */
+inline void ReSizeToMmax(int_4 mmax) 
+  {
+    nmmax_=(uint_4)l;
+    bm_.ReSize(2* nmmax_+1);
+  }
+//inline int_4 Size() const {return 2*nmmax_+1;};
+inline T& operator()(int m) {return bm_( adr_i(m));};
+inline T const& operator()(int m) const {return *(bm_.Begin()+ adr_i(m));};
+/*! return the current value of the maximum absolute value of the index */
+inline int_4 Mmax() const 
+   {
+     return (int_4)nmmax_;
+   }
+inline Bm<T>& operator = (const Bm<T>& b)
+ {
+   if (this != &b)
+     {
+       nmmax_= b.nmmax_;
+       bm_=  b.bm_;
+     }
+   return *this;
+ }
+
+  private:
+/*! return the address in the array representing the vector */
+inline uint_4 adr_i(int i) const 
+{
+  return(uint_4)(i+nmmax_);
+}
+
+
+
+ uint_4 nmmax_;
+ NDataBlock<T> bm_;
+
+  };
+
+
+#endif
Index: trunk/SophyaLib/Samba/lambdaBuilder.cc
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 729)
+++ trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 729)
@@ -0,0 +1,230 @@
+#include "lambdaBuilder.h"
+#include "nbconst.h"
+
+
+Legendre::Legendre(r_8 x, int_4 lmax) 
+{
+  if (abs(x) >1 ) 
+    {
+      throw RangeCheckError("variable for Legendre polynomials must have modules inferior to 1" );
+    }
+  x_ = x;
+  array_init(lmax);
+}
+void Legendre::array_init(int_4 lmax)
+{
+  lmax_ = lmax;
+  Pl_.ReSize(lmax_+1);
+  Pl_(0)=1.;
+  Pl_(1)=x_;
+  for (int k=2; k<Pl_.NElts(); k++)
+    {
+      Pl_(k) = ( (2.*k-1)*x_*Pl_(k-1)-(k-1)*Pl_(k-2) )/k;
+    }
+}
+
+TriangularMatrix<r_8>* LambdaLMBuilder::a_recurrence_ = NULL;
+TriangularMatrix<r_8>* LambdaLMBuilder::lam_fact_     = NULL;
+TVector<r_8>* LambdaLMBuilder::normal_l_     = NULL;
+
+LambdaLMBuilder::LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax)
+    {
+      cth_=cos(theta);
+      sth_=sin(theta);
+      array_init(lmax, mmax);
+    }
+void LambdaLMBuilder::array_init(int lmax, int mmax)
+   {
+     if (a_recurrence_ == NULL)
+       {
+	 a_recurrence_ = new TriangularMatrix<r_8>;
+	 updateArrayRecurrence(lmax);
+       }
+     else
+       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;
+	 updateArrayRecurrence(lmax);
+	 }
+     lmax_=lmax;
+     mmax_=mmax;
+     r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
+
+     lambda_.ReSizeRow(lmax_+1);
+
+     r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
+	  
+     for (int m=0; m<=mmax_;m++)
+       {
+
+
+	 lambda_(m,m)= lam_mm / bignorm2; 
+
+	 r_8 lam_0=0.;
+	 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 b_rec = 0.;
+	 for (int l=m+1; l<=lmax_; l++)
+	 {
+	   r_8 lam_2 = (cth_*lam_1-b_rec*lam_0)*a_rec;
+	   lambda_(l,m) = lam_2*lam_mm;
+	   b_rec=1./a_rec;
+	   //	   a_rec= LWK->a_recurr(l,m);
+	   a_rec= (*a_recurrence_)(l,m);
+	   lam_0 = lam_1;
+	   lam_1 = lam_2;
+	 }
+
+	 lam_mm = -lam_mm*sth_* sqrt( (2.*m+3.)/ (2.*m+2.) );
+
+       }
+   }
+
+
+void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
+   {
+     (*a_recurrence_).ReSizeRow(lmax+1);
+     for (int m=0; m<=lmax;m++)
+     {
+       (*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) );
+	 }
+     }
+   }
+
+
+void  LambdaLMBuilder::updateArrayLamNorm()
+     {
+       (*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)  );
+	     }
+	 }
+       (*normal_l_).ReSize(lmax_+1);
+       (*normal_l_)(0)=0.;
+       (*normal_l_)(1)=0.;
+       for (int l=2; l< (*normal_l_).NElts(); l++)
+	 {
+	   (*normal_l_)(l) =(r_8)sqrt( 2./( (l+2)*(l+1)*l*(l-1) ) );	
+	 }
+     }
+
+
+
+
+LambdaWXBuilder::LambdaWXBuilder(r_8 theta, int_4 lmax, int_4 mmax) : LambdaLMBuilder(theta, lmax, mmax)
+  {
+    array_init();
+  }
+
+
+void LambdaWXBuilder::array_init()
+   {
+     if (lam_fact_ == NULL || 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 )
+	 {
+	   updateArrayLamNorm();
+	 }
+
+     r_8 one_on_s2  = 1. / (sth_*sth_) ;    // 1/sin^2
+     r_8 c_on_s2    = cth_*one_on_s2;
+     lamWlm_.ReSizeRow(lmax_+1);
+     lamXlm_.ReSizeRow(lmax_+1);
+
+     // calcul des lambda
+     for(int m = 0;m<= mmax_; m++)
+       {
+	 for (int l=m; l<=lmax_; l++)
+	   {
+	     lamWlm_(l,m) =  0.;
+	     lamXlm_(l,m) =  0.;	  	
+	   }
+       }
+     for(int l = 2;l<= lmax_; l++)
+       {
+	 r_8 normal_l =  (*normal_l_)(l);	
+	 for (int m=0; m<=l; m++)
+	{
+	  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  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
+	  r_8  b_w =  c_on_s2 * lam_fact_l_m;
+	  r_8  a_x =  2. * cth_ * (l-1.);
+	  lamWlm_(l,m) =   normal_l * ( a_w * lam_lm - b_w * lam_lm1m );
+	  lamXlm_(l,m) = - normal_l * m* one_on_s2* ( a_x * lam_lm - lam_fact_l_m * lam_lm1m );
+	}
+       }
+
+   }
+
+
+LambdaPMBuilder::LambdaPMBuilder(r_8 theta, int_4 lmax, int_4 mmax) : LambdaLMBuilder(theta, lmax, mmax)
+  {
+    array_init();
+  }
+
+
+void LambdaPMBuilder::array_init()
+   {
+     if (lam_fact_ == NULL || 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 )
+	 {
+	   updateArrayLamNorm();
+	 }
+
+     r_8 one_on_s2 = 1. / (sth_*sth_) ;
+     r_8 c_on_s2 = cth_*one_on_s2;
+     lamPlm_.ReSizeRow(lmax_+1);
+     lamMlm_.ReSizeRow(lmax_+1);
+
+     // calcul des lambda
+     for(int m = 0;m<= mmax_; m++)
+       {
+	 for (int l=m; l<=lmax_; l++)
+	   {
+	     lamPlm_(l,m) =  0.;
+	     lamMlm_(l,m) =  0.;	  	
+	   }
+       }
+
+     for(int l = 2;l<= lmax_; l++)
+       {
+	 r_8 normal_l =  (*normal_l_)(l);	
+	 for (int m=0; m<=l; m++)
+	   {
+	     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  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
+	     r_8 f_w =  lam_fact_l_m/(sth_*sth_);
+	     r_8  c_w =  2*m*(l-1.) * c_on_s2;
+
+	     lamPlm_(l,m)  = normal_l * ( -(a_w+c_w) * lam_lm + f_w*( cth_ + m) * lam_lm1m )/Rac2;
+	     lamMlm_(l,m) = normal_l * ( -(a_w-c_w) * lam_lm + f_w*( cth_ - m) * lam_lm1m )/Rac2;
+	   }
+       }
+
+   }
+
+
Index: trunk/SophyaLib/Samba/lambdaBuilder.h
===================================================================
--- trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 729)
+++ trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 729)
@@ -0,0 +1,196 @@
+#ifndef LAMBDABUILDER_SEEN
+#define LAMBDABUILDER_SEEN
+
+#include <math.h>
+#include "ndatablock.h"
+#include "tvector.h" 
+#include "alm.h"
+
+/*! 
+generate Legendre polynomials : use in two steps : 
+
+a) instanciate Legendre(\f$x\f$, \f$lmax\f$) ; \f$x\f$ is the value for wich Legendre polynomials will be required (usually equal to \f$\cos \theta\f$) and \f$lmax\f$ is the MAXIMUM value of the order of polynomials wich will be required in the following code (all polynomials, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula).
+
+b) get the value of Legendre polynomial for a particular value of \f$l\f$ by calling the method getPl.
+*/
+class Legendre {
+
+ public :
+Legendre();
+Legendre(r_8 x, int_4 lmax);
+ inline r_8 getPl(int_4 l) 
+{
+  if (l>lmax_) 
+    {
+      throw (" illegal call of Legendre::getPl with index greater than lmax, which Legendre Class was instanciated with : instanciate Legendre with a greater lmax... ");
+    }
+return Pl_(l);
+}
+
+ private : 
+   /*! compute all \f$P_l(x,l_{max})\f$ for \f$l=1,l_{max}\f$ */
+  void array_init(int_4 lmax);
+
+  r_8 x_;
+  int_4 lmax_;
+  TVector<r_8> Pl_;
+};
+
+
+
+/*! 
+This class generate the coefficients :
+\f[
+            \lambda_l^m=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}
+	    P_l^m(\cos{\theta})
+\f]
+where \f$P_l^m\f$ are the associated Legendre polynomials.  The above coefficients contain the theta-dependance of spheric harmonics : 
+\f[
+            Y_{lm}(\cos{\theta})=\lambda_l^m(\cos{\theta}) e^{im\phi}.
+\f]
+
+Each object has a fixed theta (radians), and maximum l and m to be calculated
+(lmax and mmax).
+ use the class in two steps : 
+a) instanciate  LambdaLMBuilder(\f$\theta\f$, \f$lmax\f$, \f$mmax\f$) ;  \f$lmax\f$ and \f$mmax\f$ are  MAXIMUM values for which \f$\lambda_l^m\f$ will be required in the following code (all coefficients, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula).
+b) get the values of coefficients for  particular values of \f$l\f$ and \f$m\f$ by calling the method lamlm.
+*/
+ class LambdaLMBuilder {
+
+ public:
+
+LambdaLMBuilder() {}
+LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax);
+ ~LambdaLMBuilder() {};
+
+ /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
+inline double lamlm(int l, int m) const { return lambda_(l,m); }
+
+ private:
+/*! compute a static array of coefficients independant from theta (common to all instances of the LambdaBuilder Class */
+ void updateArrayRecurrence(int_4 lmax);
+
+ void array_init(int lmax, int mmax);
+
+
+ static TriangularMatrix<r_8>* a_recurrence_;
+ TriangularMatrix<r_8> lambda_;
+
+ protected :
+
+/*! compute  static arrays of coefficients independant from theta (common to all instances of the derived  classes */
+ void updateArrayLamNorm();
+
+ static  TriangularMatrix<r_8>* lam_fact_;
+ static  TVector<r_8>*  normal_l_;
+ int_4 lmax_;
+ int_4 mmax_;
+ r_8 cth_;
+ r_8 sth_;
+ 
+ };
+
+
+
+/*! 
+
+This class generates the coefficients :
+\f[
+            _{w}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}} G^+_{lm}
+\f]
+\f[
+            _{x}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}G^-_{lm}
+\f]
+where
+\f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})
+\f]
+and
+\f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)
+\f]
+ \f$P_l^m\f$ are the associated Legendre polynomials.
+
+The coefficients express the theta-dependance of the \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ functions :
+\f[W_{lm}(\cos{\theta}) = \sqrt{\frac{(l+2)!}{2(l-2)!}}_w\lambda_l^m(\cos{\theta})e^{im\phi}
+\f]
+\f[X_{lm}(\cos{\theta}) = -i\sqrt{\frac{(l+2)!}{2(l-2)!}}_x\lambda_l^m(\cos{\theta})e^{im\phi}
+\f]
+ where \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ are defined as :
+
+\f[
+W_{lm}(\cos{\theta})=-\frac{1}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(
+_{+2}Y_l^m(\cos{\theta})+_{-2}Y_l^m(\cos{\theta})\right)
+\f]
+\f[X_{lm}(\cos{\theta})=-\frac{i}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(
+_{+2}Y_l^m(\cos{\theta})-_{-2}Y_l^m(\cos{\theta})\right)
+\f]
+
+*/
+class LambdaWXBuilder : public LambdaLMBuilder
+{
+ public:
+
+
+LambdaWXBuilder() {}
+
+LambdaWXBuilder(r_8 theta, int_4 lmax, int_4 mmax);
+ 
+ /*! return the value of the coefficients \f$ _{w}\lambda_l^m\f$ and \f$_{x}\lambda_l^m\f$ */
+inline void lam_wx(int l, int m, r_8& w, r_8& x) const
+   {
+     w=lamWlm_(l,m);
+     x=lamXlm_(l,m);
+   }
+
+ private:
+
+ void array_init();
+
+
+  TriangularMatrix<r_8> lamWlm_;
+  TriangularMatrix<r_8> lamXlm_;
+  
+
+};
+
+/*! 
+
+This class generates the coefficients
+\f[
+            _{\pm}\lambda_l^m=2\sqrt{\frac{(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\left( G^+_{lm} \mp G^-_{lm}\right)
+\f]
+where
+\f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})
+\f]
+and
+\f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)
+\f]
+and \f$P_l^m\f$ are the associated Legendre polynomials.
+The coefficients express the theta-dependance of the  spin-2 spherical harmonics :
+\f[_{\pm2}Y_l^m(\cos{\theta})=_\pm\lambda_l^m(\cos{\theta})e^{im\phi}
+\f]
+*/
+class LambdaPMBuilder : public LambdaLMBuilder
+{
+ public:
+
+LambdaPMBuilder() {}
+
+ LambdaPMBuilder(r_8 theta, int_4 lmax, int_4 mmax);
+ /*! return the value of the coefficients \f$ _{+}\lambda_l^m\f$ and \f$_{-}\lambda_l^m\f$ */
+inline void lam_pm(int l, int m, r_8& lambda_plus, r_8& lambda_moins) const
+   {
+      lambda_plus = lamPlm_(l,m);
+     lambda_moins = lamMlm_(l,m);
+   }
+
+ private:
+ void array_init();
+
+  TriangularMatrix<r_8> lamPlm_;
+  TriangularMatrix<r_8> lamMlm_;
+
+
+
+};
+
+#endif
Index: trunk/SophyaLib/Samba/sphericaltransformserver.cc
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 729)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 729)
@@ -0,0 +1,862 @@
+#include "machdefs.h"
+#include <iostream.h>
+#include <math.h>
+#include <complex>
+#include "sphericaltransformserver.h"
+#include "tvector.h"
+#include "nbrandom.h"
+#include "nbmath.h"
+
+
+template<class T>
+void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
+{
+  /*=======================================================================
+    computes a map form 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)}
+    
+    where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
+    
+    * the recurrence of Ylm is the standard one (cf Num Rec)
+    * the sum over m is done by FFT
+    
+    =======================================================================*/
+  int_4 nlmax=alm.Lmax();
+  int_4 nmmax=nlmax;
+  int_4 nsmax=0;
+  map.Resize(pixelSizeIndex);
+  char* sphere_type=map.TypeOfMap();
+  if (strncmp(sphere_type,"RING",4) == 0)
+    {
+      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
+    // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
+    // c'est approximatif ; a raffiner.
+    if (strncmp(sphere_type,"TETAFI",6) == 0)
+      {
+	nsmax=(int_4)sqrt(map.NbPixels()/12.);
+      }
+    else
+      {
+	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) 
+    {
+      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
+      if (strncmp(sphere_type,"TETAFI",6) == 0)
+	{
+	  cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
+	}
+    }
+  Bm<complex<T> > b_m_theta(nmmax);
+
+  //  map.Resize(nsmax);
+
+
+  // pour chaque tranche en theta
+  for (int ith = 0; ith < map.NbThetaSlices();ith++)
+    {
+      int_4 nph;
+      r_8 phi0;
+      r_8 theta;
+      TVector<int_4> pixNumber; 
+      TVector<T> datan;
+      
+      map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
+      nph = pixNumber.NElts();
+
+      //       -----------------------------------------------------
+      //	      for each theta, and each m, computes
+      //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
+      //	      ------------------------------------------------------
+      LambdaLMBuilder lb(theta,nlmax,nmmax);
+      //  somme sur m de 0 a l'infini
+      for (int 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;
+
+	       //	 } 
+
+	    }
+	}
+
+    //        obtains the negative m of b(m,theta) (= complex conjugate)
+
+      for (int m=1;m<=nmmax;m++)
+	{
+	  //compiler doesn't have conj()
+	  b_m_theta(-m) = conj(b_m_theta(m));
+	}
+      // ---------------------------------------------------------------
+      //    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);
+      for (int i=0;i< nph;i++)
+	{
+	  // 	  map(pixNumber(i))=Temp(i).real();
+	  	  	  map(pixNumber(i))=Temp(i);
+	}
+    }
+}
+
+
+
+template<class T>
+TVector< complex<T> >  SphericalTransformServer<T>::fourierSynthesisFromB(const Bm<complex<T> >& b_m,  int_4 nph, r_8 phi0) const
+{
+  /*=======================================================================
+     dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 
+     with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
+
+     as the set of frequencies {m} is larger than nph, 
+     we wrap frequencies within {0..nph-1}
+     ie  m = k*nph + m' with m' in {0..nph-1}
+     then
+     noting bw(m') = exp(i*m'*phi0) 
+                   * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
+        with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
+     dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
+                = Fourier Transform of bw
+        is real
+
+         NB nph is not necessarily a power of 2
+
+=======================================================================*/
+  //**********************************************************************
+  // pour une valeur de phi (indexee par j) la temperature est la transformee 
+  // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
+  // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
+  //      DT/T(j) = sum_m b(m) * exp(i*m*phi(j)) 
+  // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
+  // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
+  // DT/T(j) = somme_k somme_m'  b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
+  // somme_k : de -infini a +infini
+  // somme_m' : de 0 a nph-1
+  // On echange les sommations :
+  // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
+  // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle 
+  // vaut 1.
+  // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
+  // si phi0 n'est pas nul, il y a juste un decalage a faire.
+  //**********************************************************************
+
+  TVector< complex<T> > bw(nph);
+  TVector< complex<T> > dataout(nph);
+  TVector< complex<T> > data(nph);
+
+
+  for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
+  for (int m=-b_m.Mmax();m<=-1;m++)
+    {
+      int maux=m;
+      while (maux<0) maux+=nph;
+      int iw=maux%nph;
+      double aux=(m-iw)*phi0;
+      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
+    }
+  for (int m=0;m<=b_m.Mmax();m++)
+    {
+      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
+      int iw=m%nph;
+      double aux=(m-iw)*phi0;
+      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
+    }
+
+  //     applies the shift in position <-> phase factor in Fourier space
+  for (int mprime=0; mprime < nph; mprime++)
+    {
+      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
+      data(mprime)=bw(mprime)*
+	               (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
+    }
+
+  //sortie.ReSize(nph);
+  TVector< complex<T> > sortie(nph);
+  //  sortie.SetTemp(true);
+
+  fftIntfPtr_-> FFTBackward(data, sortie);
+  
+  return sortie;
+}
+
+//********************************************
+template<class T>
+TVector<T>  SphericalTransformServer<T>::RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  int_4 nph, r_8 phi0) const
+{
+  /*=======================================================================
+     dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 
+     with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
+
+     as the set of frequencies {m} is larger than nph, 
+     we wrap frequencies within {0..nph-1}
+     ie  m = k*nph + m' with m' in {0..nph-1}
+     then
+     noting bw(m') = exp(i*m'*phi0) 
+                   * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
+        with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
+     dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
+                = Fourier Transform of bw
+        is real
+
+         NB nph is not necessarily a power of 2
+
+=======================================================================*/
+  //**********************************************************************
+  // pour une valeur de phi (indexee par j) la temperature est la transformee 
+  // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
+  // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
+  //      DT/T(j) = sum_m b(m) * exp(i*m*phi(j)) 
+  // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
+  // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
+  // DT/T(j) = somme_k somme_m'  b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
+  // somme_k : de -infini a +infini
+  // somme_m' : de 0 a nph-1
+  // On echange les sommations :
+  // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
+  // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle 
+  // vaut 1.
+  // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
+  // si phi0 n'est pas nul, il y a juste un decalage a faire.
+  //**********************************************************************
+
+  TVector< complex<T> > bw(nph);
+  TVector< complex<T> > dataout(nph);
+  TVector< complex<T> > data(nph/2+1);
+
+
+  for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
+  for (int m=-b_m.Mmax();m<=-1;m++)
+    {
+      int maux=m;
+      while (maux<0) maux+=nph;
+      int iw=maux%nph;
+      double aux=(m-iw)*phi0;
+      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
+    }
+  for (int m=0;m<=b_m.Mmax();m++)
+    {
+      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
+      int iw=m%nph;
+      double aux=(m-iw)*phi0;
+      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
+    }
+
+  //     applies the shift in position <-> phase factor in Fourier space
+  for (int mprime=0; mprime <= nph/2; mprime++)
+    {
+      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
+      data(mprime)=bw(mprime)*
+	               (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
+    }
+
+  //sortie.ReSize(nph);
+  TVector<T> sortie;
+  //  sortie.SetTemp(true);
+
+  fftIntfPtr_-> FFTBackward(data, sortie);
+  
+  return sortie;
+}
+//*******************************************
+
+template<class T>
+ Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
+{
+  
+  /*-----------------------------------------------------------------------
+    computes the integral in phi : phas_m(theta)
+    for each parallele from north to south pole
+    -----------------------------------------------------------------------*/
+  TVector<T> data;
+  TVector<int_4> pixNumber;
+  int_4  nmmax = nlmax;
+  TVector< complex<T> > phase(nmmax+1);
+  Alm<T> alm;
+  //  alm.SetTemp(true);
+  alm.ReSizeToLmax(nlmax);
+  //  cout << " nombre de tranches " << map.NbThetaSlices() << endl;
+  for (int ith = 0; ith < map.NbThetaSlices(); ith++)
+    {
+      int_4 nph;
+      r_8 phi0;
+      r_8 theta;
+      map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
+      for (int i=0;i< nmmax+1;i++)
+	{
+	  phase(i)=0; 
+	}
+      nph = pixNumber.NElts(); 
+      double cth = cos(theta);
+      
+      //part of the sky out of the symetric cut
+      bool keep_it = (abs(cth) >= cos_theta_cut); 
+      
+      if (keep_it)
+	{
+	  // tableau datain a supprimer
+	  //  TVector<complex<T> > datain(nph);
+	  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);
+
+	  //  phase = CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
+	  phase = CFromFourierAnalysis(nlmax,nmmax,data,phi0);
+
+	}
+      
+      /*-----------------------------------------------------------------------
+	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;
+}
+template<class T>
+TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nlmax,int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
+{
+  /*=======================================================================
+    integrates (data * phi-dependence-of-Ylm) over phi
+    --> function of m can be computed by FFT
+    
+    datain est modifie
+    =======================================================================*/
+  int_4 nph=datain.NElts();
+  if (nph <= 0) 
+    {
+      throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
+    }
+  TVector<complex<T> > transformedData(nph);
+  fftIntfPtr_-> FFTForward(datain, transformedData);
+
+  //dataout.ReSize(nmmax+1);
+  TVector< complex<T> > dataout(nmmax+1);
+  //  dataout.SetTemp(true);
+
+  int im_max=min(nph,nmmax+1);
+  for (int i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
+  for (int 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 (int i = 0;i <dataout.NElts();i++){
+    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
+  }
+  return dataout;
+}
+
+//&&&&&&&&& nouvelle version
+template<class T>
+TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nlmax,int_4 nmmax, const TVector<T> datain, r_8 phi0) const
+{
+  //=======================================================================
+  //    integrates (data * phi-dependence-of-Ylm) over phi
+  //    --> function of m can be computed by FFT
+  //   !     with  0<= m <= npoints/2 (: Nyquist)
+  //   !     because the data is real the negative m are the conjugate of the 
+  //   !     positive ones
+    
+  //    datain est modifie
+  //    
+  //    =======================================================================
+  int_4 nph=datain.NElts();
+  if (nph <= 0) 
+    {
+      throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
+    }
+  TVector<complex<T> > transformedData;
+    // a remodifier
+  //FFTPackServer ffts;
+  //ffts.setNormalize(false);
+  //ffts.FFTForward(datain, transformedData);
+
+  fftIntfPtr_-> FFTForward(datain, transformedData);
+    //
+
+  //dataout.ReSize(nmmax+1);
+  TVector< complex<T> > dataout(nmmax+1);
+  //  dataout.SetTemp(true); 
+
+  // 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);
+  for (int i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);
+
+
+  // si dataout n'est pas plein, on complete jusqu'a  nph valeurs (a moins 
+  // que dataout ne soit plein avant d'atteindre nph)
+  if (sizeOfTransformToGet == (transformedData.NElts()))
+    {
+      for (int i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
+	{
+
+	  //	  dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
+	  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 (int i = 0;i <dataout.NElts();i++){
+    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
+  }
+  return dataout;
+}
+
+template<class T>
+void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq,
+					       SphericalMap<T>& mapu, 
+					       int_4 pixelSizeIndex,
+					       const Alm<T>& alme,
+					       const Alm<T>& almb) const
+{
+  /*=======================================================================
+    computes a map form 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)}
+    
+    where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
+    
+    * the recurrence of Ylm is the standard one (cf Num Rec)
+    * the sum over m is done by FFT
+    
+    =======================================================================*/
+  int_4 nlmax=alme.Lmax();
+  if (nlmax != almb.Lmax())
+    {
+      cout << " SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille" << endl;
+      throw SzMismatchError("SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille");
+    }
+  int_4 nmmax=nlmax;
+  int_4 nsmax=0;
+  mapq.Resize(pixelSizeIndex);
+  mapu.Resize(pixelSizeIndex);
+  char* sphere_type=mapq.TypeOfMap();
+  if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
+    {
+      cout <<  " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
+      cout << " type 1 " << sphere_type << endl;
+      cout << " type 2 " << mapu.TypeOfMap() << endl;
+      throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
+      
+    }
+  if (strncmp(sphere_type,"RING",4) == 0)
+    {
+      nsmax=mapq.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
+    // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
+    // c'est approximatif ; a raffiner.
+    if (strncmp(sphere_type,"TETAFI",6) == 0)
+      {
+	nsmax=(int_4)sqrt(mapq.NbPixels()/12.);
+      }
+    else
+      {
+	cout << " unknown type of sphere : " << sphere_type << endl;
+	throw IOExc(" unknown type of sphere ");
+      }
+  cout << "GenerateFromAlm: the spheres are of type : " << sphere_type << endl;
+  cout << "GenerateFromAlm: size indices (nside) of  spheres= " << nsmax << endl;
+  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
+  if (nlmax>3*nsmax-1) 
+    {
+      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
+      if (strncmp(sphere_type,"TETAFI",6) == 0)
+	{
+	  cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
+	}
+    }
+  if (alme.Lmax()!=almb.Lmax())
+    {
+      cout << "GenerateFromAlm: arrays Alme and Almb have not the same size ? " << endl; 
+      throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ?  ");
+    }
+    mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
+    // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
+}
+
+
+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
+{
+  int_4  nmmax = nlmax;
+  // resize et remise a zero
+  alme.ReSizeToLmax(nlmax);
+  almb.ReSizeToLmax(nlmax);
+
+  
+  TVector<T> dataq;
+  TVector<T> datau;
+  TVector<int_4> pixNumber;
+
+  char* sphere_type=mapq.TypeOfMap();
+  if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
+    {
+      cout <<  " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
+      cout << " type 1 " << sphere_type << endl;
+      cout << " type 2 " << mapu.TypeOfMap() << endl;
+      throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
+      
+    }
+  if (mapq.NbPixels()!=mapu.NbPixels())
+    {
+      cout << " DecomposeToAlm: map Q and map U have not same size ?" << endl;
+      throw SzMismatchError("SphericalTransformServer::DecomposeToAlm: map Q and map U have not same size ");
+    }
+  for (int ith = 0; ith < mapq.NbThetaSlices(); ith++)
+    {
+      int_4 nph;
+      r_8 phi0;
+      r_8 theta;
+      mapq.GetThetaSlice(ith,theta,phi0, pixNumber,dataq);
+      mapu.GetThetaSlice(ith,theta,phi0, pixNumber,datau);
+      if (dataq.NElts() != datau.NElts() ) 
+	{
+	  throw  SzMismatchError("the spheres have not the same pixelization");
+	}
+      nph = pixNumber.NElts();
+      r_8 domega=mapq.PixSolAngle(mapq.PixIndexSph(theta,phi0));
+      double cth = cos(theta);
+      //part of the sky out of the symetric cut
+      bool keep_it = (abs(cth) >= cos_theta_cut); 
+      if (keep_it)
+	{
+	  //  almFromPM(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
+	  almFromWX(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
+	}
+    }
+}
+
+
+template<class T>
+void SphericalTransformServer<T>::almFromWX(int_4 nph, int_4 nlmax, int_4 nmmax,
+					 r_8 phi0, r_8 domega, 
+					 r_8 theta, 
+					 const TVector<T>& dataq, 
+					 const TVector<T>& datau,
+					 Alm<T>& alme,
+					 Alm<T>& almb) const
+{
+  TVector< complex<T> > phaseq(nmmax+1);
+  TVector< complex<T> > phaseu(nmmax+1);
+  //  TVector<complex<T> > datain(nph);
+  for (int i=0;i< nmmax+1;i++)
+    {
+      phaseq(i)=0; 
+      phaseu(i)=0; 
+    }
+  //  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
+
+  // phaseq = CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
+  phaseq = CFromFourierAnalysis(nlmax,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(nlmax,nmmax,datau,phi0); 
+
+  LambdaWXBuilder lwxb(theta,nlmax,nmmax);
+
+  r_8 sqr2inv=1/Rac2;
+  for (int m = 0; m <= nmmax; m++)
+    {
+      r_8 lambda_w=0.;
+      r_8 lambda_x=0.;
+      lwxb.lam_wx(m, m, lambda_w, lambda_x);
+      complex<T>  zi_lam_x((T)0., (T)lambda_x);
+      alme(m,m) +=  ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
+      almb(m,m) +=  ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
+      
+      for (int l = m+1; l<= nlmax; l++)
+	{
+	  lwxb.lam_wx(l, m, lambda_w, lambda_x);
+	  zi_lam_x = complex<T>((T)0., (T)lambda_x);
+	  alme(l,m) +=  ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
+	  almb(l,m) +=  ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
+	}
+    }
+}
+
+
+template<class T>
+void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax,
+					 r_8 phi0, r_8 domega,  
+					 r_8 theta, 
+					 const TVector<T>& dataq, 
+					 const TVector<T>& datau,
+					 Alm<T>& alme,
+					 Alm<T>& almb) const
+{
+  TVector< complex<T> > phasep(nmmax+1);
+  TVector< complex<T> > phasem(nmmax+1);
+  TVector<complex<T> > datain(nph);
+  for (int i=0;i< nmmax+1;i++)
+    {
+      phasep(i)=0; 
+      phasem(i)=0; 
+    }
+
+  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));
+
+  phasep = CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
+
+  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),-datau(kk));
+  phasem = CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
+  LambdaPMBuilder lpmb(theta,nlmax,nmmax);
+	  
+  for (int m = 0; m <= nmmax; m++)
+    {
+      r_8 lambda_p=0.;
+      r_8 lambda_m=0.;
+      complex<T> im((T)0.,(T)1.);
+      lpmb.lam_pm(m, m, lambda_p, lambda_m);
+	      
+      alme(m,m) +=   -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m)  )*(T)(domega*0.5);
+      almb(m,m) +=  im*( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
+      for (int l = m+1; l<= nlmax; l++)
+	{
+	  lpmb.lam_pm(l, m, lambda_p, lambda_m);
+	  alme(l,m) +=  -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m)  )*(T)(domega*0.5);
+	  almb(l,m) += im* ( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
+	}
+    }
+}
+
+
+template<class T>
+void SphericalTransformServer<T>::mapFromWX(int_4 nlmax, int_4 nmmax,
+					 SphericalMap<T>& mapq,
+					 SphericalMap<T>& mapu, 
+					 const Alm<T>& alme,
+					 const Alm<T>& almb) const
+{
+  Bm<complex<T> > b_m_theta_q(nmmax);
+  Bm<complex<T> > b_m_theta_u(nmmax);
+
+  for (int ith = 0; ith < mapq.NbThetaSlices();ith++)
+    {
+      int_4 nph;
+      r_8 phi0;
+      r_8 theta;
+      TVector<int_4>  pixNumber; 
+      TVector<T> datan;
+      
+      mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
+      nph =  pixNumber.NElts();
+      //       -----------------------------------------------------
+      //	      for each theta, and each m, computes
+      //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
+      //	      ------------------------------------------------------
+      LambdaWXBuilder lwxb(theta,nlmax,nmmax);
+      //      LambdaPMBuilder lpmb(theta,nlmax,nmmax);
+      r_8 sqr2inv=1/Rac2;
+      for (int m = 0; m <= nmmax; m++)
+	{
+	  r_8 lambda_w=0.;
+	  r_8 lambda_x=0.;
+	  lwxb.lam_wx(m, m, lambda_w, lambda_x);
+	  complex<T>  zi_lam_x((T)0., (T)lambda_x);
+	  
+	  b_m_theta_q(m) =  ( (T)(lambda_w) * alme(m,m) - zi_lam_x * almb(m,m))*(T)sqr2inv ;
+	  b_m_theta_u(m) =  ( (T)(lambda_w) * almb(m,m) + zi_lam_x * alme(m,m))*(T)sqr2inv;
+	  
+	  
+	  for (int l = m+1; l<= nlmax; l++)
+	    {
+	      
+	      lwxb.lam_wx(l, m, lambda_w, lambda_x);
+	      zi_lam_x= complex<T>((T)0., (T)lambda_x);
+	      
+	      b_m_theta_q(m) += ((T)(lambda_w)*alme(l,m)-zi_lam_x *almb(l,m))*(T)sqr2inv;
+	      b_m_theta_u(m) += ((T)(lambda_w)*almb(l,m)+zi_lam_x *alme(l,m))*(T)sqr2inv;
+	      
+	    } 
+	}
+      //        obtains the negative m of b(m,theta) (= complex conjugate)
+      for (int m=1;m<=nmmax;m++)
+	{
+	  b_m_theta_q(-m) = conj(b_m_theta_q(m));
+	  b_m_theta_u(-m) = conj(b_m_theta_u(m));
+	}
+
+      //    TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);  
+      //     TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0); 
+      TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);  
+      TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0); 
+      for (int i=0;i< nph;i++)
+	{
+	  //	  mapq(pixNumber(i))=Tempq(i).real();
+	  //	  mapu(pixNumber(i))=Tempu(i).real();
+	  mapq(pixNumber(i))=Tempq(i);
+	  mapu(pixNumber(i))=Tempu(i);
+	  
+	}
+    } 
+}
+template<class T>
+void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax,
+					 SphericalMap<T>& mapq,
+					 SphericalMap<T>& mapu, 
+					 const Alm<T>& alme,
+					 const Alm<T>& almb) const
+{
+  Bm<complex<T> > b_m_theta_p(nmmax);
+  Bm<complex<T> > b_m_theta_m(nmmax);
+  for (int ith = 0; ith < mapq.NbThetaSlices();ith++)
+    {
+      int_4 nph;
+      r_8 phi0;
+      r_8 theta;
+      TVector<int_4> pixNumber; 
+      TVector<T> datan;
+      
+      mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
+      nph =  pixNumber.NElts();
+
+      //       -----------------------------------------------------
+      //	      for each theta, and each m, computes
+      //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
+      //------------------------------------------------------
+
+      LambdaPMBuilder lpmb(theta,nlmax,nmmax);
+      for (int m = 0; m <= nmmax; m++)
+	{
+	  r_8 lambda_p=0.;
+	  r_8 lambda_m=0.;
+	  lpmb.lam_pm(m, m, lambda_p, lambda_m);
+	  complex<T> im((T)0.,(T)1.);
+	  
+	  b_m_theta_p(m) =  (T)(lambda_p )* (-alme(m,m) - im * almb(m,m));
+	  b_m_theta_m(m) =  (T)(lambda_m) * (-alme(m,m) + im * almb(m,m));
+	  
+	  
+	  for (int l = m+1; l<= nlmax; l++)
+	    {
+	      lpmb.lam_pm(l, m, lambda_p, lambda_m);
+	      b_m_theta_p(m) +=  (T)(lambda_p)*(-alme(l,m)-im *almb(l,m));
+	      b_m_theta_m(m) +=  (T)(lambda_m)*(-alme(l,m)+im *almb(l,m));
+	    }
+	}
+      
+      //        obtains the negative m of b(m,theta) (= complex conjugate)
+      for (int m=1;m<=nmmax;m++)
+	{
+	  b_m_theta_p(-m) = conj(b_m_theta_m(m));
+	  b_m_theta_m(-m) = conj(b_m_theta_p(m));
+	}
+
+      TVector<complex<T> > Tempp = fourierSynthesisFromB(b_m_theta_p,nph,phi0);  
+      TVector<complex<T> > Tempm = fourierSynthesisFromB(b_m_theta_m,nph,phi0); 
+
+      for (int i=0;i< nph;i++)
+	{
+	  	  mapq(pixNumber(i))=0.5*(Tempp(i)+Tempm(i)).real();
+	  	  mapu(pixNumber(i))=0.5*(Tempp(i)-Tempm(i)).imag();
+	}
+    }
+}
+
+
+template<class T>
+void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq, 
+					      SphericalMap<T>& sphu, 
+					      int_4 pixelSizeIndex,
+					      const TVector<T>& Cle, 
+					      const TVector<T>& Clb, 
+					      const r_8 fwhm) const
+{
+  if (Cle.NElts() != Clb.NElts())
+    {
+      cout << " SphericalTransformServer: les deux tableaux Cl n'ont pas la meme taille" << endl;
+      throw SzMismatchError("SphericalTransformServer::GenerateFromCl :  two Cl arrays have not same size");
+    }
+
+  //  Alm<T> a2lme,a2lmb;
+  //  almFromCl(a2lme, Cle, fwhm); 
+  //  almFromCl(a2lmb, Clb, fwhm); 
+  //  Alm<T> a2lme = almFromCl(Cle, fwhm);
+  // Alm<T> a2lmb = almFromCl(Clb, fwhm);
+  Alm<T> a2lme(Cle, fwhm);
+  Alm<T> a2lmb(Clb, fwhm);
+
+  GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb); 
+}
+template<class T>
+void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph,
+						 int_4 pixelSizeIndex, 
+					      const TVector<T>& Cl, 
+						 const r_8 fwhm)  const
+{
+
+  //  Alm<T> alm;
+  //  almFromCl(alm, Cl, fwhm); 
+  //Alm<T> alm = almFromCl(Cl, fwhm);
+  Alm<T> alm(Cl, fwhm);
+  GenerateFromAlm(sph,pixelSizeIndex, alm ); 
+}
+
+
+
+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);
+    Alm<T> alm=DecomposeToAlm( sph, nlmax, cos_theta_cut);
+  // power spectrum
+     return  alm.powerSpectrum();
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template SphericalTransformServer<r_8>
+#pragma define_template SphericalTransformServer<r_4>
+#endif
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class SphericalTransformServer<r_8>;
+template class SphericalTransformServer<r_4>;
+#endif
Index: trunk/SophyaLib/Samba/sphericaltransformserver.h
===================================================================
--- trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 729)
+++ trunk/SophyaLib/Samba/sphericaltransformserver.h	(revision 729)
@@ -0,0 +1,204 @@
+#ifndef SPHERICALTRANFORMSERVER_SEEN
+#define SPHERICALTRANFORMSERVER_SEEN
+
+#include "sphericalmap.h"
+#include "spheregorski.h"
+#include "fftservintf.h"
+#include "fftpserver.h"
+#include "alm.h"
+#include "lambdaBuilder.h"
+
+
+
+//
+/*! Class for performing analysis and synthesis of sky maps using spin-0 or spin-2 spherical harmonics.
+ */
+template <class T>
+class SphericalTransformServer
+{
+
+ public:
+
+ SphericalTransformServer()  
+{
+    fftIntfPtr_=new FFTPackServer;
+    fftIntfPtr_->setNormalize(false);
+};
+ ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;};
+/*!
+ Set a fft server. The constructor sets a default fft server (fft-pack). So it is not necessary to call this method for a standard use.
+*/
+ void SetFFTServer(FFTServerInterface* srv=NULL) 
+{
+  if (fftIntfPtr_!=NULL) delete fftIntfPtr_;
+  fftIntfPtr_=srv;
+  fftIntfPtr_->setNormalize(false);
+}
+ /*! synthesis of a temperature  map from  Alm coefficients */
+ void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; 
+ /*! synthesis of a polarization map from  Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */
+ void GenerateFromAlm(SphericalMap<T>& mapq, SphericalMap<T>& mapu, int_4 pixelSizeIndex, const Alm<T>& alme,    const Alm<T>& almb) const;
+
+ /*! synthesis of a temperature  map from  power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */
+ void GenerateFromCl(SphericalMap<T>& sph, int_4 pixelSizeIndex,
+		     const TVector<T>& Cl, const r_8 fwhm) const; 
+  /*! synthesis of a polarization  map from  power spectra electric-Cl and magnetic-Cl (Alm's are generated randomly, following a gaussian distribution). 
+  \param fwhm FWHM in arcmin for random generation of Alm's (eg. 5) 
+
+*/
+ void GenerateFromCl(SphericalMap<T>& sphq, SphericalMap<T>& sphu, 
+		     int_4 pixelSizeIndex,
+		     const TVector<T>& Cle, const TVector<T>& Clb, 
+		    const r_8 fwhm) const; 
+ /*!return the Alm coefficients from analysis of a temperature map.
+
+    nlmax : maximum value of the l index
+
+     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.
+  */ 
+
+
+Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const;
+ /*analysis of a polarization map into Alm coefficients.
+
+ The spheres mapq and mapu contain respectively the Stokes parameters. 
+
+ a2lme and a2lmb will receive respectively electric and magnetic Alm's
+    nlmax : maximum value of the l index
+
+    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.
+ */ 
+
+ 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;
+
+ // /*return power spectrum from analysis of a temperature map.
+
+ //    nlmax : maximum value of the l index
+
+ //    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.
+ // */ 
+ TVector<T>  DecomposeToCl(const SphericalMap<T>& sph,  
+			   int_4 nlmax, r_8 cos_theta_cut) const;
+
+
+  private:
+  /*! return a vector with nph elements  which are  sums :\f$\sum_{m=-mmax}^{mmax}b_m(\theta)e^{im\varphi}\f$ for nph values of \f$\varphi\f$ regularly distributed in \f$[0,\pi]\f$ ( calculated by FFT)
+
+  The object b_m (\f$b_m\f$) of the class Bm is a special vector which index goes from -mmax to mmax. 
+  */
+ TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+				  int_4 nph, r_8 phi0) const;
+/*! same as fourierSynthesisFromB, but return a real vector, taking into account the fact that b(-m) is conjugate of b(m) */
+ TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+				  int_4 nph, r_8 phi0) const;
+
+  /*! return a vector with mmax elements  which are  sums :
+\f$\sum_{k=0}^{nphi}datain(\theta,\varphi_k)e^{im\varphi_k}\f$ for (mmax+1) values of \f$m\f$ from 0 to mmax.
+   */
+ TVector< complex<T> > CFromFourierAnalysis(int_4 nlmax,int_4 mmax, 
+				  const TVector<complex<T> > datain,
+				  r_8 phi0) const;
+/* same as previous one, but with a "datain" which is real (not complex) */
+ TVector< complex<T> > CFromFourierAnalysis(int_4 nlmax,int_4 mmax, 
+			 const TVector<T> datain,  
+				  r_8 phi0) const;
+
+ /*! 
+Compute polarized Alm's as : 
+\f[
+a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(_{w}\lambda_l^m\tilde{Q}-i_{x}\lambda_l^m\tilde{U}\right)}
+\f]
+\f[
+a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(i_{x}\lambda_l^m\tilde{Q}+_{w}\lambda_l^m\tilde{U}\right)}
+\f]
+
+where \f$\tilde{Q}\f$ and \f$\tilde{U}\f$ are C-coefficients computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters.
+
+\f$\omega_{pix}\f$ are solid angle of each pixel.
+
+dataq, datau : Stokes parameters.
+
+  */
+void almFromWX(int_4 nph, int_4 nlmax, int_4 nmmax, r_8 phi0, 
+	       r_8 domega, r_8 theta, 
+	       const TVector<T>& dataq, const TVector<T>& datau,
+	       Alm<T>& alme, Alm<T>& almb) const;
+ /*! 
+Compute polarized Alm's as : 
+\f[
+a_{lm}^E=-\frac{1}{2}\sum_{slices}{\omega_{pix}\left(_{+}\lambda_l^m\tilde{P^+}+_{-}\lambda_l^m\tilde{P^-}\right)}
+\f]
+\f[
+a_{lm}^B=\frac{i}{2}\sum_{slices}{\omega_{pix}\left(_{+}\lambda_l^m\tilde{P^+}-_{-}\lambda_l^m\tilde{P^-}\right)}
+\f]
+
+where \f$\tilde{P^{\pm}}=\tilde{Q}\pm\tilde{U}\f$  computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters,\f$Q\f$ and \f$U\f$ .
+
+\f$\omega_{pix}\f$ are solid angle of each pixel.
+
+dataq, datau : Stokes parameters.
+
+  */
+void almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax, 
+	       r_8 phi0, r_8 domega, r_8 theta, 
+	       const TVector<T>& dataq, const TVector<T>& datau,
+	       Alm<T>& alme, Alm<T>& almb) const;
+
+/*! synthesis of Stokes parameters following formulae : 
+
+\f[
+Q=\sum_{m=-mmax}^{mmax}b_m^qe^{im\varphi}
+\f]
+\f[
+U=\sum_{m=-mmax}^{mmax}b_m^ue^{im\varphi}
+\f]
+
+computed by FFT (method fourierSynthesisFromB called by the present one)
+
+with :
+
+\f[
+b_m^q=-\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(_{w}\lambda_l^ma_{lm}^E-i_{x}\lambda_l^ma_{lm}^B\right) }
+\f]
+\f[
+b_m^u=\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(i_{x}\lambda_l^ma_{lm}^E+_{w}\lambda_l^ma_{lm}^B\right) }
+\f]
+ */
+void mapFromWX(int_4 nlmax, int_4 nmmax,
+	       SphericalMap<T>& mapq, SphericalMap<T>& mapu, 
+	       const Alm<T>& alme, const Alm<T>& almb) const;
+
+/*! synthesis of polarizations following formulae : 
+
+\f[
+P^+ = \sum_{m=-mmax}^{mmax} {b_m^+e^{im\varphi} }
+\f]
+\f[
+P^- = \sum_{m=-mmax}^{mmax} {b_m^-e^{im\varphi} }
+\f]
+
+computed by FFT (method fourierSynthesisFromB called by the present one)
+
+with :
+
+\f[
+b_m^+=-\sum_{l=|m|}^{lmax}{_{+}\lambda_l^m \left( a_{lm}^E+ia_{lm}^B \right) }
+\f]
+\f[
+b_m^-=-\sum_{l=|m|}^{lmax}{_{+}\lambda_l^m \left( a_{lm}^E-ia_{lm}^B \right) }
+\f]
+ */
+
+void mapFromPM(int_4 nlmax, int_4 nmmax,
+	       SphericalMap<T>& mapq, SphericalMap<T>& mapu, 
+	       const Alm<T>& alme, const Alm<T>& almb) const;
+
+
+
+ FFTServerInterface* fftIntfPtr_;
+};
+
+
+#endif
