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


template<class T>
void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
{
  /*=======================================================================
    computes a map form its alm for the HEALPIX pixelisation
    map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
    = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
    
    where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
    
    * the recurrence of Ylm is the standard one (cf Num Rec)
    * the sum over m is done by FFT
    
    =======================================================================*/
  int_4 nlmax=alm.Lmax();
  int_4 nmmax=nlmax;
  int_4 nsmax=0;
  map.Resize(pixelSizeIndex);
  char* sphere_type=map.TypeOfMap();
  if (strncmp(sphere_type,"RING",4) == 0)
    {
      nsmax=map.SizeIndex();
    }
  else
    // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
    // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
    // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
    // c'est approximatif ; a raffiner.
    if (strncmp(sphere_type,"TETAFI",6) == 0)
      {
	nsmax=(int_4)sqrt(map.NbPixels()/12.);
      }
    else
      {
	cout << " unknown type of sphere : " << sphere_type << endl;
	throw IOExc(" unknown type of sphere: " + (string)sphere_type );
      }
  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
  if (nlmax>3*nsmax-1) 
    {
      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
      if (strncmp(sphere_type,"TETAFI",6) == 0)
	{
	  cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
	}
    }
  Bm<complex<T> > b_m_theta(nmmax);

  //  map.Resize(nsmax);


  // pour chaque tranche en theta
  for (int_4 ith = 0; ith < map.NbThetaSlices();ith++)
    {
      int_4 nph;
      r_8 phi0;
      r_8 theta;
      TVector<int_4> pixNumber; 
      TVector<T> datan;
      
      map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
      nph = pixNumber.NElts();

      //       -----------------------------------------------------
      //	      for each theta, and each m, computes
      //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
      //	      ------------------------------------------------------
      LambdaLMBuilder lb(theta,nlmax,nmmax);
      //  somme sur m de 0 a l'infini
      for (int m = 0; m <= nmmax; m++)
	{
	  //  somme sur l de m a l'infini
	  b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
	  //    if (ith==0 && m==0)
	  //	{
	  //	  cout << " guy: lmm= " <<  lb.lamlm(m,m) << " alm " << alm(m,m) << "b00= " <<  b_m_theta(m) << endl;
	  //	} 
	  for (int l = m+1; l<= nlmax; l++)
	    {
	       b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);


	       //     if (ith==0 && m==0)
	       //	 {
	       //	   cout << " guy:l=" << l << " m= " << m << " lmm= " <<  lb.lamlm(l,m) << " alm " << alm(l,m) << "b00= " <<  b_m_theta(m) << endl;

	       //	 } 

	    }
	}

    //        obtains the negative m of b(m,theta) (= complex conjugate)

      for (int m=1;m<=nmmax;m++)
	{
	  //compiler doesn't have conj()
	  b_m_theta(-m) = conj(b_m_theta(m));
	}
      // ---------------------------------------------------------------
      //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
      // ---------------------------------------------------------------*/
      //    TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0); 
	  TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
      for (int i=0;i< nph;i++)
	{
	  // 	  map(pixNumber(i))=Temp(i).real();
	  	  	  map(pixNumber(i))=Temp(i);
	}
    }
}



template<class T>
TVector< complex<T> >  SphericalTransformServer<T>::fourierSynthesisFromB(const Bm<complex<T> >& b_m,  int_4 nph, r_8 phi0) const
{
  /*=======================================================================
     dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 
     with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1

     as the set of frequencies {m} is larger than nph, 
     we wrap frequencies within {0..nph-1}
     ie  m = k*nph + m' with m' in {0..nph-1}
     then
     noting bw(m') = exp(i*m'*phi0) 
                   * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
        with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
     dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
                = Fourier Transform of bw
        is real

         NB nph is not necessarily a power of 2

=======================================================================*/
  //**********************************************************************
  // pour une valeur de phi (indexee par j) la temperature est la transformee 
  // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
  // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
  //      DT/T(j) = sum_m b(m) * exp(i*m*phi(j)) 
  // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
  // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
  // DT/T(j) = somme_k somme_m'  b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
  // somme_k : de -infini a +infini
  // somme_m' : de 0 a nph-1
  // On echange les sommations :
  // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
  // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle 
  // vaut 1.
  // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
  // si phi0 n'est pas nul, il y a juste un decalage a faire.
  //**********************************************************************

  TVector< complex<T> > bw(nph);
  TVector< complex<T> > dataout(nph);
  TVector< complex<T> > data(nph);


  for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
  for (int m=-b_m.Mmax();m<=-1;m++)
    {
      int maux=m;
      while (maux<0) maux+=nph;
      int iw=maux%nph;
      double aux=(m-iw)*phi0;
      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
    }
  for (int m=0;m<=b_m.Mmax();m++)
    {
      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
      int iw=m%nph;
      double aux=(m-iw)*phi0;
      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
    }

  //     applies the shift in position <-> phase factor in Fourier space
  for (int mprime=0; mprime < nph; mprime++)
    {
      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
      data(mprime)=bw(mprime)*
	               (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
    }

  //sortie.ReSize(nph);
  TVector< complex<T> > sortie(nph);
  sortie.SetTemp(true);

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

//********************************************
template<class T>
TVector<T>  SphericalTransformServer<T>::RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  int_4 nph, r_8 phi0) const
{
  /*=======================================================================
     dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 
     with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1

     as the set of frequencies {m} is larger than nph, 
     we wrap frequencies within {0..nph-1}
     ie  m = k*nph + m' with m' in {0..nph-1}
     then
     noting bw(m') = exp(i*m'*phi0) 
                   * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
        with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
     dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
                = Fourier Transform of bw
        is real

         NB nph is not necessarily a power of 2

=======================================================================*/
  //**********************************************************************
  // pour une valeur de phi (indexee par j) la temperature est la transformee 
  // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
  // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
  //      DT/T(j) = sum_m b(m) * exp(i*m*phi(j)) 
  // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
  // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
  // DT/T(j) = somme_k somme_m'  b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
  // somme_k : de -infini a +infini
  // somme_m' : de 0 a nph-1
  // On echange les sommations :
  // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
  // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle 
  // vaut 1.
  // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
  // si phi0 n'est pas nul, il y a juste un decalage a faire.
  //**********************************************************************

  TVector< complex<T> > bw(nph);
  TVector< complex<T> > dataout(nph);
  TVector< complex<T> > data(nph/2+1);


  for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
  for (int m=-b_m.Mmax();m<=-1;m++)
    {
      int maux=m;
      while (maux<0) maux+=nph;
      int iw=maux%nph;
      double aux=(m-iw)*phi0;
      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
    }
  for (int m=0;m<=b_m.Mmax();m++)
    {
      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
      int iw=m%nph;
      double aux=(m-iw)*phi0;
      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
    }

  //     applies the shift in position <-> phase factor in Fourier space
  for (int mprime=0; mprime <= nph/2; mprime++)
    {
      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
      data(mprime)=bw(mprime)*
	               (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
    }

  //sortie.ReSize(nph);
  TVector<T> sortie;
  sortie.SetTemp(true);

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

template<class T>
 Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
{
  
  /*-----------------------------------------------------------------------
    computes the integral in phi : phas_m(theta)
    for each parallele from north to south pole
    -----------------------------------------------------------------------*/
  TVector<T> data;
  TVector<int_4> pixNumber;
  int_4  nmmax = nlmax;
  TVector< complex<T> > phase(nmmax+1);
  Alm<T> alm;
  alm.SetTemp(true);
  alm.ReSizeToLmax(nlmax);
  for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
    {
      int_4 nph;
      r_8 phi0;
      r_8 theta;
      map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
      for (int i=0;i< nmmax+1;i++)
	{
	  phase(i)=0; 
	}
      nph = pixNumber.NElts(); 
      double cth = cos(theta);
      
      //part of the sky out of the symetric cut
      bool keep_it = (abs(cth) >= cos_theta_cut); 
      
      if (keep_it)
	{
	  // tableau datain a supprimer
	  //  TVector<complex<T> > datain(nph);
	  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);

	  //  phase = CFromFourierAnalysis(nmmax,datain,phi0);
	  phase = CFromFourierAnalysis(nmmax,data,phi0);

	}
      
      /*-----------------------------------------------------------------------
	computes the a_lm by integrating over theta
	lambda_lm(theta) * phas_m(theta)
	for each m and l
	-----------------------------------------------------------------------*/
      //      LambdaBuilder lb(theta,nlmax,nmmax);
      LambdaLMBuilder lb(theta,nlmax,nmmax);
      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
      for (int m = 0; m <= nmmax; m++)
	{
	  alm(m,m) += (T)lb.lamlm(m,m) * phase(m) * (T)domega; //m,m even
	  for (int l = m+1; l<= nlmax; l++)
	    {
	      alm(l,m) += (T)lb.lamlm(l,m) * phase(m)*(T)domega; 
	    }
	}
    }
  return alm;
}
template<class T>
TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
{
  /*=======================================================================
    integrates (data * phi-dependence-of-Ylm) over phi
    --> function of m can be computed by FFT
    
    datain est modifie
    =======================================================================*/
  int_4 nph=datain.NElts();
  if (nph <= 0) 
    {
      throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
    }
  TVector<complex<T> > transformedData(nph);
  fftIntfPtr_-> FFTForward(datain, transformedData);

  //dataout.ReSize(nmmax+1);
  TVector< complex<T> > dataout(nmmax+1);
  //  dataout.SetTemp(true);

  int im_max=min(nph,nmmax+1);
  for (int i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
  for (int i=0;i<im_max;i++) dataout(i)=transformedData(i);


  //  for (int i = 0;i <im_max;i++){
  //    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
  //  }
  for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
  for (int i = 0;i <dataout.NElts();i++){
    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
  }
  return dataout;
}

//&&&&&&&&& nouvelle version
template<class T>
TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const
{
  //=======================================================================
  //    integrates (data * phi-dependence-of-Ylm) over phi
  //    --> function of m can be computed by FFT
  //   !     with  0<= m <= npoints/2 (: Nyquist)
  //   !     because the data is real the negative m are the conjugate of the 
  //   !     positive ones
    
  //    datain est modifie
  //    
  //    =======================================================================
  int_4 nph=datain.NElts();
  if (nph <= 0) 
    {
      throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
    }
  TVector<complex<T> > transformedData;
    // a remodifier
  //FFTPackServer ffts;
  //ffts.setNormalize(false);
  //ffts.FFTForward(datain, transformedData);

  fftIntfPtr_-> FFTForward(datain, transformedData);
    //

  //dataout.ReSize(nmmax+1);
  TVector< complex<T> > dataout(nmmax+1);
  //  dataout.SetTemp(true); 

  // on transfere le resultat de la fft dans dataout.
  // on s'assure que ca ne depasse pas la taille de dataout
  int sizeOfTransformToGet = min(transformedData.NElts(),nmmax+1);
  //  int im_max=min(transformedData.NElts()-1,nmmax);
  for (int i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);


  // si dataout n'est pas plein, on complete jusqu'a  nph valeurs (a moins 
  // que dataout ne soit plein avant d'atteindre nph)
  if (sizeOfTransformToGet == (transformedData.NElts()))
    {
      for (int i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
	{

	  //	  dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
	  dataout(i) = conj(dataout(nph-i) );
	}
      // on conplete, si necessaire, par periodicite
      for (int kk=nph; kk<dataout.NElts(); kk++) 
	{
	  dataout(kk)=dataout(kk%nph);
	}
    }
  for (int i = 0;i <dataout.NElts();i++){
    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
  }
  return dataout;
}

template<class T>
void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq,
					       SphericalMap<T>& mapu, 
					       int_4 pixelSizeIndex,
					       const Alm<T>& alme,
					       const Alm<T>& almb) const
{
  /*=======================================================================
    computes a map form its alm for the HEALPIX pixelisation
    map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
    = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
    
    where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
    
    * the recurrence of Ylm is the standard one (cf Num Rec)
    * the sum over m is done by FFT
    
    =======================================================================*/
  int_4 nlmax=alme.Lmax();
  if (nlmax != almb.Lmax())
    {
      cout << " SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille" << endl;
      throw SzMismatchError("SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille");
    }
  int_4 nmmax=nlmax;
  int_4 nsmax=0;
  mapq.Resize(pixelSizeIndex);
  mapu.Resize(pixelSizeIndex);
  char* sphere_type=mapq.TypeOfMap();
  if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
    {
      cout <<  " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
      cout << " type 1 " << sphere_type << endl;
      cout << " type 2 " << mapu.TypeOfMap() << endl;
      throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
      
    }
  if (strncmp(sphere_type,"RING",4) == 0)
    {
      nsmax=mapq.SizeIndex();
    }
  else
    // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
    // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
    // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
    // c'est approximatif ; a raffiner.
    if (strncmp(sphere_type,"TETAFI",6) == 0)
      {
	nsmax=(int_4)sqrt(mapq.NbPixels()/12.);
      }
    else
      {
	cout << " unknown type of sphere : " << sphere_type << endl;
	throw IOExc(" unknown type of sphere ");
      }
  cout << "GenerateFromAlm: the spheres are of type : " << sphere_type << endl;
  cout << "GenerateFromAlm: size indices (nside) of  spheres= " << nsmax << endl;
  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
  if (nlmax>3*nsmax-1) 
    {
      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
      if (strncmp(sphere_type,"TETAFI",6) == 0)
	{
	  cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
	}
    }
  if (alme.Lmax()!=almb.Lmax())
    {
      cout << "GenerateFromAlm: arrays Alme and Almb have not the same size ? " << endl; 
      throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ?  ");
    }
    mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
    // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
}


template<class T>
void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
					      const SphericalMap<T>& mapu,
					      Alm<T>& alme,
					      Alm<T>& almb,
					      int_4 nlmax,
					      r_8 cos_theta_cut) const
{
  int_4  nmmax = nlmax;
  // resize et remise a zero
  alme.ReSizeToLmax(nlmax);
  almb.ReSizeToLmax(nlmax);

  
  TVector<T> dataq;
  TVector<T> datau;
  TVector<int_4> pixNumber;

  char* sphere_type=mapq.TypeOfMap();
  if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
    {
      cout <<  " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
      cout << " type 1 " << sphere_type << endl;
      cout << " type 2 " << mapu.TypeOfMap() << endl;
      throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
      
    }
  if (mapq.NbPixels()!=mapu.NbPixels())
    {
      cout << " DecomposeToAlm: map Q and map U have not same size ?" << endl;
      throw SzMismatchError("SphericalTransformServer::DecomposeToAlm: map Q and map U have not same size ");
    }
  for (int_4 ith = 0; ith < mapq.NbThetaSlices(); ith++)
    {
      int_4 nph;
      r_8 phi0;
      r_8 theta;
      mapq.GetThetaSlice(ith,theta,phi0, pixNumber,dataq);
      mapu.GetThetaSlice(ith,theta,phi0, pixNumber,datau);
      if (dataq.NElts() != datau.NElts() ) 
	{
	  throw  SzMismatchError("the spheres have not the same pixelization");
	}
      nph = pixNumber.NElts();
      r_8 domega=mapq.PixSolAngle(mapq.PixIndexSph(theta,phi0));
      double cth = cos(theta);
      //part of the sky out of the symetric cut
      bool keep_it = (abs(cth) >= cos_theta_cut); 
      if (keep_it)
	{
	  //  almFromPM(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
	  almFromWX(nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb); 
	}
    }
}


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,datain,phi0); 
  phaseq = CFromFourierAnalysis(nmmax,dataq,phi0); 

  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(datau(kk),0.);

  // phaseu=  CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 

  LambdaWXBuilder lwxb(theta,nlmax,nmmax);

  r_8 sqr2inv=1/Rac2;
  for (int m = 0; m <= nmmax; m++)
    {
      r_8 lambda_w=0.;
      r_8 lambda_x=0.;
      lwxb.lam_wx(m, m, lambda_w, lambda_x);
      complex<T>  zi_lam_x((T)0., (T)lambda_x);
      alme(m,m) +=  ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
      almb(m,m) +=  ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
      
      for (int l = m+1; l<= nlmax; l++)
	{
	  lwxb.lam_wx(l, m, lambda_w, lambda_x);
	  zi_lam_x = complex<T>((T)0., (T)lambda_x);
	  alme(l,m) +=  ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
	  almb(l,m) +=  ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
	}
    }
}


template<class T>
void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax,
					 r_8 phi0, r_8 domega,  
					 r_8 theta, 
					 const TVector<T>& dataq, 
					 const TVector<T>& datau,
					 Alm<T>& alme,
					 Alm<T>& almb) const
{
  TVector< complex<T> > phasep(nmmax+1);
  TVector< complex<T> > phasem(nmmax+1);
  TVector<complex<T> > datain(nph);
  for (int i=0;i< nmmax+1;i++)
    {
      phasep(i)=0; 
      phasem(i)=0; 
    }

  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));

  phasep = CFromFourierAnalysis(nmmax,datain,phi0); 

  for(int 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);
	}
    }
}


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;
      for (int m = 0; m <= nmmax; m++)
	{
	  r_8 lambda_w=0.;
	  r_8 lambda_x=0.;
	  lwxb.lam_wx(m, m, lambda_w, lambda_x);
	  complex<T>  zi_lam_x((T)0., (T)lambda_x);
	  
	  b_m_theta_q(m) =  ( (T)(lambda_w) * alme(m,m) - zi_lam_x * almb(m,m))*(T)sqr2inv ;
	  b_m_theta_u(m) =  ( (T)(lambda_w) * almb(m,m) + zi_lam_x * alme(m,m))*(T)sqr2inv;
	  
	  
	  for (int l = m+1; l<= nlmax; l++)
	    {
	      
	      lwxb.lam_wx(l, m, lambda_w, lambda_x);
	      zi_lam_x= complex<T>((T)0., (T)lambda_x);
	      
	      b_m_theta_q(m) += ((T)(lambda_w)*alme(l,m)-zi_lam_x *almb(l,m))*(T)sqr2inv;
	      b_m_theta_u(m) += ((T)(lambda_w)*almb(l,m)+zi_lam_x *alme(l,m))*(T)sqr2inv;
	      
	    } 
	}
      //        obtains the negative m of b(m,theta) (= complex conjugate)
      for (int m=1;m<=nmmax;m++)
	{
	  b_m_theta_q(-m) = conj(b_m_theta_q(m));
	  b_m_theta_u(-m) = conj(b_m_theta_u(m));
	}

      //    TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);  
      //     TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0); 
      TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);  
      TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0); 
      for (int i=0;i< nph;i++)
	{
	  //	  mapq(pixNumber(i))=Tempq(i).real();
	  //	  mapu(pixNumber(i))=Tempu(i).real();
	  mapq(pixNumber(i))=Tempq(i);
	  mapu(pixNumber(i))=Tempu(i);
	  
	}
    } 
}
template<class T>
void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax,
					 SphericalMap<T>& mapq,
					 SphericalMap<T>& mapu, 
					 const Alm<T>& alme,
					 const Alm<T>& almb) const
{
  Bm<complex<T> > b_m_theta_p(nmmax);
  Bm<complex<T> > b_m_theta_m(nmmax);
  for (int_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);
      for (int m = 0; m <= nmmax; m++)
	{
	  r_8 lambda_p=0.;
	  r_8 lambda_m=0.;
	  lpmb.lam_pm(m, m, lambda_p, lambda_m);
	  complex<T> im((T)0.,(T)1.);
	  
	  b_m_theta_p(m) =  (T)(lambda_p )* (-alme(m,m) - im * almb(m,m));
	  b_m_theta_m(m) =  (T)(lambda_m) * (-alme(m,m) + im * almb(m,m));
	  
	  
	  for (int l = m+1; l<= nlmax; l++)
	    {
	      lpmb.lam_pm(l, m, lambda_p, lambda_m);
	      b_m_theta_p(m) +=  (T)(lambda_p)*(-alme(l,m)-im *almb(l,m));
	      b_m_theta_m(m) +=  (T)(lambda_m)*(-alme(l,m)+im *almb(l,m));
	    }
	}
      
      //        obtains the negative m of b(m,theta) (= complex conjugate)
      for (int m=1;m<=nmmax;m++)
	{
	  b_m_theta_p(-m) = conj(b_m_theta_m(m));
	  b_m_theta_m(-m) = conj(b_m_theta_p(m));
	}

      TVector<complex<T> > Tempp = fourierSynthesisFromB(b_m_theta_p,nph,phi0);  
      TVector<complex<T> > Tempm = fourierSynthesisFromB(b_m_theta_m,nph,phi0); 

      for (int i=0;i< nph;i++)
	{
	  	  mapq(pixNumber(i))=0.5*(Tempp(i)+Tempm(i)).real();
	  	  mapu(pixNumber(i))=0.5*(Tempp(i)-Tempm(i)).imag();
	}
    }
}


template<class T>
void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq, 
					      SphericalMap<T>& sphu, 
					      int_4 pixelSizeIndex,
					      const TVector<T>& Cle, 
					      const TVector<T>& Clb, 
					      const r_8 fwhm) const
{
  if (Cle.NElts() != Clb.NElts())
    {
      cout << " SphericalTransformServer: les deux tableaux Cl n'ont pas la meme taille" << endl;
      throw SzMismatchError("SphericalTransformServer::GenerateFromCl :  two Cl arrays have not same size");
    }

  //  Alm<T> a2lme,a2lmb;
  //  almFromCl(a2lme, Cle, fwhm); 
  //  almFromCl(a2lmb, Clb, fwhm); 
  //  Alm<T> a2lme = almFromCl(Cle, fwhm);
  // Alm<T> a2lmb = almFromCl(Clb, fwhm);
  Alm<T> a2lme(Cle, fwhm);
  Alm<T> a2lmb(Clb, fwhm);

  GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb); 
}
template<class T>
void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph,
						 int_4 pixelSizeIndex, 
					      const TVector<T>& Cl, 
						 const r_8 fwhm)  const
{

  Alm<T> alm(Cl, fwhm);
  GenerateFromAlm(sph,pixelSizeIndex, alm ); 
}



template <class T>
TVector<T>  SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
{
    Alm<T> alm=DecomposeToAlm( sph, 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