source: Sophya/trunk/SophyaLib/Samba/sphericaltransformserver.h@ 988

#ifndef SPHERICALTRANFORMSERVER_SEEN
#define SPHERICALTRANFORMSERVER_SEEN

#include "sphericalmap.h"
#include "fftservintf.h"
#include "fftpserver.h"
#include "alm.h"
#include "lambdaBuilder.h"


namespace SOPHYA {

//
/*! 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]

 */
template <class T>
class SphericalTransformServer
{

 public:

 SphericalTransformServer()  
{
    fftIntfPtr_=new FFTPackServer;
    fftIntfPtr_->setNormalize(false);
};
 ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;};
/*!
 Set a fft server. The constructor sets a default fft server (fft-pack). So it is not necessary to call this method for a standard use.
*/
 void SetFFTServer(FFTServerInterface* srv=NULL) 
{
  if (fftIntfPtr_!=NULL) delete fftIntfPtr_;
  fftIntfPtr_=srv;
  fftIntfPtr_->setNormalize(false);
}
 /*! synthesis of a temperature  map from  Alm coefficients */
 void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; 
 /*! synthesis of a polarization map from  Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */
 void GenerateFromAlm(SphericalMap<T>& mapq, SphericalMap<T>& mapu, int_4 pixelSizeIndex, const Alm<T>& alme,    const Alm<T>& almb) const;

 /*! synthesis of a temperature  map from  power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */
 void GenerateFromCl(SphericalMap<T>& sph, int_4 pixelSizeIndex,
                     const TVector<T>& Cl, const r_8 fwhm) const; 
  /*! synthesis of a polarization  map from  power spectra electric-Cl and magnetic-Cl (Alm's are generated randomly, following a gaussian distribution). 
  \param fwhm FWHM in arcmin for random generation of Alm's (eg. 5) 

*/
 void GenerateFromCl(SphericalMap<T>& sphq, SphericalMap<T>& sphu, 
                     int_4 pixelSizeIndex,
                     const TVector<T>& Cle, const TVector<T>& Clb, 
                    const r_8 fwhm) const; 
 /*!return the Alm coefficients from analysis of a temperature map.

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


Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const;
 /*!analysis of a polarization map into Alm coefficients.

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

 void DecomposeToAlm(const SphericalMap<T>& mapq, const SphericalMap<T>& mapu,
                     Alm<T>& a2lme, Alm<T>& a2lmb,
                     int_4 nlmax, r_8 cos_theta_cut) const;

/*!return power spectrum from analysis of a temperature map.

     \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.
  */ 
 TVector<T>  DecomposeToCl(const SphericalMap<T>& sph,  
                           int_4 nlmax, r_8 cos_theta_cut) const;


  private:
  /*! return a vector with nph elements  which are  sums :\f$\sum_{m=-mmax}^{mmax}b_m(\theta)e^{im\varphi}\f$ for nph values of \f$\varphi\f$ regularly distributed in \f$[0,\pi]\f$ ( calculated by FFT)

  The object b_m (\f$b_m\f$) of the class Bm is a special vector which index goes from -mmax to mmax. 
  */
 TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
                                  int_4 nph, r_8 phi0) const;
/*! same as fourierSynthesisFromB, but return a real vector, taking into account the fact that b(-m) is conjugate of b(m) */
 TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
                                  int_4 nph, r_8 phi0) const;

  /*! return a vector with mmax elements  which are  sums :
\f$\sum_{k=0}^{nphi}datain(\theta,\varphi_k)e^{im\varphi_k}\f$ for (mmax+1) values of \f$m\f$ from 0 to mmax.
   */
 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
                                  const TVector<complex<T> > datain,
                                  r_8 phi0) const;
/* same as previous one, but with a "datain" which is real (not complex) */
 TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
                         const TVector<T> datain,  
                                  r_8 phi0) const;

 /*! 
Compute polarized Alm's as : 
\f[
a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(\,_{w}\lambda_l^m\tilde{Q}-i\,_{x}\lambda_l^m\tilde{U}\right)}
\f]
\f[
a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(i\,_{x}\lambda_l^m\tilde{Q}+\,_{w}\lambda_l^m\tilde{U}\right)}
\f]

where \f$\tilde{Q}\f$ and \f$\tilde{U}\f$ are C-coefficients computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters.

\f$\omega_{pix}\f$ are solid angle of each pixel.

dataq, datau : Stokes parameters.

  */
void almFromWX(int_4 nlmax, int_4 nmmax, r_8 phi0, 
               r_8 domega, r_8 theta, 
               const TVector<T>& dataq, const TVector<T>& datau,
               Alm<T>& alme, Alm<T>& almb) const;
 /*! 
Compute polarized Alm's as : 
\f[
a_{lm}^E=-\frac{1}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}+\,_{-}\lambda_l^m\tilde{P^-}\right)}
\f]
\f[
a_{lm}^B=\frac{i}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}-\,_{-}\lambda_l^m\tilde{P^-}\right)}
\f]

where \f$\tilde{P^{\pm}}=\tilde{Q}\pm\tilde{U}\f$  computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters,\f$Q\f$ and \f$U\f$ .

\f$\omega_{pix}\f$ are solid angle of each pixel.

dataq, datau : Stokes parameters.

  */
void almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax, 
               r_8 phi0, r_8 domega, r_8 theta, 
               const TVector<T>& dataq, const TVector<T>& datau,
               Alm<T>& alme, Alm<T>& almb) const;

/*! synthesis of Stokes parameters following formulae : 

\f[
Q=\sum_{m=-mmax}^{mmax}b_m^qe^{im\varphi}
\f]
\f[
U=\sum_{m=-mmax}^{mmax}b_m^ue^{im\varphi}
\f]

computed by FFT (method fourierSynthesisFromB called by the present one)

with :

\f[
b_m^q=-\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(\,_{w}\lambda_l^ma_{lm}^E-i\,_{x}\lambda_l^ma_{lm}^B\right) }
\f]
\f[
b_m^u=\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(i\,_{x}\lambda_l^ma_{lm}^E+\,_{w}\lambda_l^ma_{lm}^B\right) }
\f]
 */
void mapFromWX(int_4 nlmax, int_4 nmmax,
               SphericalMap<T>& mapq, SphericalMap<T>& mapu, 
               const Alm<T>& alme, const Alm<T>& almb) const;

/*! synthesis of polarizations following formulae : 

\f[
P^+ = \sum_{m=-mmax}^{mmax} {b_m^+e^{im\varphi} }
\f]
\f[
P^- = \sum_{m=-mmax}^{mmax} {b_m^-e^{im\varphi} }
\f]

computed by FFT (method fourierSynthesisFromB called by the present one)

with :

\f[
b_m^+=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E+ia_{lm}^B \right) }
\f]
\f[
b_m^-=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E-ia_{lm}^B \right) }
\f]
 */

void mapFromPM(int_4 nlmax, int_4 nmmax,
               SphericalMap<T>& mapq, SphericalMap<T>& mapu, 
               const Alm<T>& alme, const Alm<T>& almb) const;



 FFTServerInterface* fftIntfPtr_;
};
} // Fin du namespace


#endif