#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <math.h>
#include <complex>
#include "sphericaltransformserver.h"
#include "tvector.h"
#include "nbrandom.h"
#include "nbmath.h"
#include "timing.h"
//#include "spherehealpix.h"


/*! 
  \ingroup Samba
  \class SOPHYA::SphericalTransformServer
  
  \brief Analysis/synthesis in spherical harmonics server.

 Class for performing analysis and synthesis of sky maps using spin-0 or spin-2 spherical harmonics.

Maps must be SOPHYA SphericalMaps (SphereHEALPix or SphereThetaPhi or SphereECP).

Temperature and polarization (Stokes parameters) can be developped on spherical harmonics : 
\f[
\frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n})
\f]
\f[
Q(\hat{n})=\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EW_{lm}(\hat{n})+a_{lm}^BX_{lm}(\hat{n})\right)
\f]
\f[
U(\hat{n})=-\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EX_{lm}(\hat{n})-a_{lm}^BW_{lm}(\hat{n})\right)
\f]
\f[
\left(Q \pm iU\right)(\hat{n})=\sum_{lm}a_{\pm 2lm}\, _{\pm 2}Y_l^m(\hat{n})
\f]

\f[
Y_l^m(\hat{n})=\lambda_l^m(\theta)e^{im\phi}
\f]
\f[
_{\pm}Y_l^m(\hat{n})=_{\pm}\lambda_l^m(\theta)e^{im\phi}
\f]
\f[
W_{lm}(\hat{n})=\frac{1}{N_l}\,_{w}\lambda_l^m(\theta)e^{im\phi}
\f]
\f[
X_{lm}(\hat{n})=\frac{-i}{N_l}\,_{x}\lambda_l^m(\theta)e^{im\phi}
\f]

(see LambdaLMBuilder, LambdaPMBuilder, LambdaWXBuilder classes)

power spectra : 

\f[
C_l^T=\frac{1}{2l+1}\sum_{m=0}^{+ \infty }\left|a_{lm}^T\right|^2=\langle\left|a_{lm}^T\right|^2\rangle 
\f]
\f[
C_l^E=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^E\right|^2=\langle\left|a_{lm}^E\right|^2\rangle 
\f]
\f[
C_l^B=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^B\right|^2=\langle\left|a_{lm}^B\right|^2\rangle 
\f]

\arg
\b Synthesis : Get temperature and polarization maps  from \f$a_{lm}\f$ coefficients or from power spectra, (methods GenerateFrom...).

\b Temperature:
\f[
\frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n}) = \sum_{-\infty}^{+\infty}b_m(\theta)e^{im\phi}
\f]

with 
\f[
b_m(\theta)=\sum_{l=\left|m\right|}^{+\infty}a_{lm}^T\lambda_l^m(\theta)
\f]

\b Polarisation
\f[
Q \pm iU = \sum_{-\infty}^{+\infty}b_m^{\pm}(\theta)e^{im\phi}
\f]

where :
\f[
b_m^{\pm}(\theta) = \sum_{l=\left|m\right|}^{+\infty}a_{\pm 2lm}\,_{\pm}\lambda_l^m(\theta)
\f]

or :
\f[
Q  = \sum_{-\infty}^{+\infty}b_m^{Q}(\theta)e^{im\phi}
\f]
\f[
U  = \sum_{-\infty}^{+\infty}b_m^{U}(\theta)e^{im\phi}
\f]

where: 
\f[
b_m^{Q}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(a_{lm}^E\,_{w}\lambda_l^m(\theta)-ia_{lm}^B\,_{x}\lambda_l^m(\theta)\right)
\f]
\f[
b_m^{U}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(ia_{lm}^E\,_{x}\lambda_l^m(\theta)+a_{lm}^B\,_{w}\lambda_l^m(\theta)\right)
\f]

Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed  on \f$2\pi\f$  \f$\frac{\Delta T}{T}\f$, \f$Q\f$,\f$U\f$  can be computed by FFT.


\arg
\b Analysis :  Get \f$a_{lm}\f$ coefficients or  power spectra from temperature and polarization maps   (methods DecomposeTo...). 

\b Temperature:
\f[
a_{lm}^T=\int\frac{\Delta T}{T}(\hat{n})Y_l^{m*}(\hat{n})d\hat{n}
\f]

approximated as : 
\f[
a_{lm}^T=\sum_{\theta_k}\omega_kC_m(\theta_k)\lambda_l^m(\theta_k)
\f]
where :
\f[
C_m (\theta _k)=\sum_{\phi _{k\prime}}\frac{\Delta T}{T}(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
\f]
Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed  on \f$2\pi\f$ (\f$\omega_k\f$ is the solid angle of each pixel of the slice \f$\theta_k\f$) \f$C_m\f$ can be computed by FFT.

\b polarisation:

\f[
a_{\pm 2lm}=\sum_{\theta_k}\omega_kC_m^{\pm}(\theta_k)\,_{\pm}\lambda_l^m(\theta_k)
\f]
where :
\f[
C_m^{\pm} (\theta _k)=\sum_{\phi _{k\prime}}\left(Q \pm iU\right)(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
\f]
or :

\f[
a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(C_m^{Q}(\theta_k)\,_{w}\lambda_l^m(\theta_k)-iC_m^{U}(\theta_k)\,_{x}\lambda_l^m(\theta_k)\right)
\f]
\f[
a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(iC_m^{Q}(\theta_k)\,_{x}\lambda_l^m(\theta_k)+C_m^{U}(\theta_k)\,_{w}\lambda_l^m(\theta_k)\right)
\f]

where : 
\f[
C_m^{Q} (\theta _k)=\sum_{\phi _{k\prime}}Q(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
\f]
\f[
C_m^{U} (\theta _k)=\sum_{\phi _{k\prime}}U(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
\f]

 */

 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const

 synthesis of a temperature  map from  Alm coefficients 
*/
template<class T>
void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
{
  /*=======================================================================
    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)}
    
    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;
  // le Resize est suppose mettre a zero
  map.Resize(pixelSizeIndex);
  string sphere_type=map.TypeOfMap();
  int premiereTranche = 0;
  int derniereTranche = map.NbThetaSlices()-1;

  Bm<complex<T> > b_m_theta(nmmax);

  // pour chaque tranche en theta
  for (int_4 ith = premiereTranche; ith <= derniereTranche;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();
    if (nph < 2) continue;  // On laisse tomber les tranches avec un point
    //       -----------------------------------------------------
    //	      for each theta, and each m, computes
    //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
    //	      ------------------------------------------------------
    // ===> Optimisation Reza, Mai 2006 
    /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
      qui calcule la somme au vol 
      LambdaLMBuilder lb(theta,nlmax,nmmax);
      //  somme sur m de 0 a l'infini
      for (int_4 m = 0; m <= nmmax; m++) {
      b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
      for (int l = m+1; l<= nlmax; l++)
        b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
      }
      ------- Fin version PRE-Mai2006 */
    LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
    //Fin Optimisation Reza, Mai 2006 <==== 

    //        obtains the negative m of b(m,theta) (= complex conjugate)
    for (int_4 m=1;m<=nmmax;m++)
      b_m_theta(-m) = conj(b_m_theta(m));
    // ---------------------------------------------------------------
    //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
    // ---------------------------------------------------------------*/

    /* ----- Reza, Juin 2006 : 
       En verifiant la difference entre deux cartes 
       cl -> map -> alm -> map2 et mapdiff = map-map2 
       je me suis apercu qu'il y avait des differences importantes - dans les
       deux zones 'polar cap' de HEALPix - qui utilisait RfourierSynthesisFromB
       TF complex -> reel . Le probleme venant de l'ambiguite de taille, lie 
       a la partie imaginaire de la composante a f_nyquist , j'ai corrige et 
       tout mis en TF complexe -> reel 
    */
    TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
    // Si on peut acceder directement les pixels d'un tranche, on le fait
    T* pix = map.GetThetaSliceDataPtr(ith);
    if (pix != NULL) 
      for (int_4 i=0;i< nph;i++) pix[i] = Temp(i);
    else 
      for (int_4 i=0;i< nph;i++) map(pixNumber(i))=Temp(i); 
  }
}



  /*! \fn TVector< complex<T> >  SOPHYA::SphericalTransformServer::fourierSynthesisFromB(const Bm<complex<T> >& b_m,  int_4 nph, r_8 phi0) const

\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. 
  */
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_m' (exp(i*m'*phi(j)) somme_k 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.;
  int m;
  for (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 (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);

  fftIntfPtr_-> FFTBackward(data, sortie);
  
  return sortie;
}

//********************************************
/*! \fn TVector<T>  SOPHYA::SphericalTransformServer::RfourierSynthesisFromB(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) */
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_m' (exp(i*m'*phi(j)) somme_k 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_k 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> > data(nph/2+1);

  for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
  int m;
  for (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 (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++)  
    data(mprime)=bw(mprime)*complex<T>((T)cos(mprime*phi0),(T)sin(mprime*phi0));    
  TVector<T> sortie(nph);
// On met la partie imaginaire du dernier element du data a zero pour nph pair
  if (nph%2 == 0) data(nph/2) = complex<T>(data(nph/2).real(), (T)0.);
// et on impose l'utilisation de la taille en sortie pour FFTBack (..., ..., true)
  fftIntfPtr_-> FFTBackward(data, sortie, true); 
  return sortie;
}
//*******************************************

 /*! \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. 

    \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>
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
{
  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;
      int k;
      for (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 (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
{
  
  /*-----------------------------------------------------------------------
    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.ReSizeToLmax(nlmax);
  for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
    {
      r_8 phi0;
      r_8 theta;
      //  PrtTim("debut 1ere tranche ");
      map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
      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 ");

	}
      
//      ---------------------------------------------------------------------
//	computes the a_lm by integrating over theta
//	lambda_lm(theta) * phas_m(theta)
//	for each m and l
//	-----------------------------------------------------------------------

      // ===> Optimisation Reza, Mai 2006 
      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
        qui calcule la somme au vol 
      //        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 ; 
	    }
	}
      ------- Fin version PRE-Mai2006 */
      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
      phase *= complex<T>((T)domega, 0.);
      LambdaLMBuilder::ComputeAlmFrPhase(theta,nlmax,nmmax, phase, alm);
      //Fin Optimisation Reza, Mai 2006 <==== 
      

      
      //
      //
      //       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

\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.
   */
template<class T>
TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(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);

  TVector< complex<T> > dataout(nmmax+1);

  int im_max=min(nph,nmmax+1);
  int i;
  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 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;
}

//&&&&&&&&& nouvelle version
/* \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const

same as previous one, but with a "datain" which is real (not complex) */
template<class T>
void SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, TVector< complex<T> >& dataout, 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)");
    }
  // 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);

  // on transfere le resultat de la fft dans dataout.

  int maxFreqAccessiblesParFFT = min(nph/2,nmmax);
  int i;
  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 (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
      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;
}

 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm(SphericalMap<T>& mapq,
					       SphericalMap<T>& mapu, 
					       int_4 pixelSizeIndex,
					       const Alm<T>& alme,
					       const Alm<T>& almb) const

synthesis of a polarization map from  Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */
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);
  string sphere_type=mapq.TypeOfMap();
  if (sphere_type != mapu.TypeOfMap())
    {
      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");
      
    }
  bool healpix = true;
  if (sphere_type.substr(0,4) == "RING")
    {
      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.
    healpix = false;
    if (sphere_type.substr(0,6) == "TETAFI")
      {
	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 (sphere_type.substr(0,6) == "TETAFI")
	{
	  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, healpix);
    // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
}
 /*! \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) 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.


 */ 
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>
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;
      int k;
      for (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 (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,
					      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;

  string sphere_type=mapq.TypeOfMap();
  if (sphere_type != mapu.TypeOfMap())
    {
      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_4 ith = 0; ith < mapq.NbThetaSlices(); ith++)
    {
      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");
	}
      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 = (fabs(cth) >= cos_theta_cut); 
      if (keep_it)
	{
	  //  almFromPM(pixNumber.NElts(), nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
	  almFromWX(nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
	}
    }
}


 /*! \fn void SOPHYA::SphericalTransformServer::almFromWX(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}{\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.

  */
template<class T>
void SphericalTransformServer<T>::almFromWX(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(nmmax,dataq,phi0); 
  CFromFourierAnalysis(nmmax,dataq,phaseq, phi0); 

  //  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 
  CFromFourierAnalysis(nmmax,datau,phaseu, 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);
	}
    }
}


 /*! \fn void SOPHYA::SphericalTransformServer::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

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.

  */
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; 
    }
  int kk;
  for(kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));

  phasep = CFromFourierAnalysis(nmmax,datain,phi0); 

  for(kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),-datau(kk));
  phasem = CFromFourierAnalysis(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);
	}
    }
}


/*! \fn void SOPHYA::SphericalTransformServer::mapFromWX(int_4 nlmax, int_4 nmmax,
					 SphericalMap<T>& mapq,
					 SphericalMap<T>& mapu, 
					 const Alm<T>& alme,
					 const Alm<T>& almb, bool healpix) 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]
 */
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, bool healpix) const
{
  int i;

  Bm<complex<T> > b_m_theta_q(nmmax);
  Bm<complex<T> > b_m_theta_u(nmmax);

  for (int_4 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;
      int m;
      for (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 (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));
	}
      if (healpix) 
	{
	  TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);  
	  TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0); 
	  for (i=0;i< nph;i++)
	    {
	      mapq(pixNumber(i))=Tempq(i);
	      mapu(pixNumber(i))=Tempu(i);
	    }
	}
      else
	// pour des pixelisations quelconques (autres que HEALPix
	//  nph n'est pas toujours pair
	// ca fait des problemes pour les transformees de Fourier
	// car le server de TF ajuste la longueur du vecteur reel 
	// en sortie de TF, bref, la securite veut qu'on prenne une 
	// TF complexe
	{
          TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);  
	  TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0); 
	  for (i=0;i< nph;i++)
	    {
	      mapq(pixNumber(i))=Tempq(i).real();
	      mapu(pixNumber(i))=Tempu(i).real();
	    }
	}
    } 
}
/*! \fn void SOPHYA::SphericalTransformServer::mapFromPM(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]
 */
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_4 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);
      int m;
      for (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 (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();
	}
    }
}


  /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sphq, 
					      SphericalMap<T>& sphu, 
					      int_4 pixelSizeIndex,
					      const TVector<T>& Cle, 
					      const TVector<T>& Clb, 
					      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) 
*/
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); 
}
 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sph,
						 int_4 pixelSizeIndex, 
					      const TVector<T>& Cl, 
						 const r_8 fwhm)  const

synthesis of a temperature  map from  power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */
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(Cl, fwhm);
  GenerateFromAlm(sph,pixelSizeIndex, alm ); 
}



/*! \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) 

     \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. If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !

  */ 
template <class T>
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);
  // 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();
}

#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 SOPHYA::SphericalTransformServer<r_8>;
template class SOPHYA::SphericalTransformServer<r_4>;
#endif