source: Sophya/trunk/SophyaLib/Samba/sphericaltransformserver.cc@ 1827

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


/*! \class SOPHYA::SphericalTransformServer

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

Maps must be SOPHYA SphericalMaps (SphereGorski or SphereThetaPhi).

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);
  char* sphere_type=map.TypeOfMap();
      int premiereTranche = 0;
      int derniereTranche = map.NbThetaSlices()-1;
  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.);
      premiereTranche++;
      derniereTranche--;
    }
    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 << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
  //      cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
  //    }
  //}
  }
  Bm<complex<T> > b_m_theta(nmmax);

  //  map.Resize(nsmax);


  // 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();

      //       -----------------------------------------------------
      //              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
      int m;
      for (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);
            }
        }
    //        obtains the negative m of b(m,theta) (= complex conjugate)

      for (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)
      // ---------------------------------------------------------------*/
      TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
      for (int 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_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.;
  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_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.;
  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++)
    {
      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
      data(mprime)=bw(mprime)*
                       (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
    }

  TVector<T> sortie;
  fftIntfPtr_-> FFTBackward(data, sortie);
  
  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
//      -----------------------------------------------------------------------
      //        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 ; 
            }
        }
      

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

  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_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) 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) const
{
  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));
        }

      //    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);
          
        }
    } 
}
/*! \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 SphericalTransformServer<r_8>;
template class SphericalTransformServer<r_4>;
#endif