source: Sophya/trunk/SophyaLib/Samba/alm.cc@ 3667

#include "sopnamsp.h"
#include "alm.h"


/*!
  \class Alm 
  \ingroup Samba 
  Class for manipulating the coefficients \f$a_{lm}\f$ of the development
  in spherical harmonics of  a function efined on a sphere.
*/
 
/*!
  fwhm specifies the gaussian beam half witdh in arc.minutes
*/
template <class T>
Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm) 

{
  int_4 nlmax= clin.NElts()-1;

  //alm.ReSize(nlmax);
  this->ReSizeRow(nlmax+1); 
  GenFromCl(clin, fwhm, *(RandomGeneratorInterface::GetGlobalRandGenP()) );
}

/*!
  fwhm specifies the gaussian beam half witdh in arc.minutes
*/
template <class T>
Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm, RandomGeneratorInterface & rg)
{
  int_4 nlmax= clin.NElts()-1;

  //alm.ReSize(nlmax);
  this->ReSizeRow(nlmax+1);
  GenFromCl(clin, fwhm, rg);  
}



template <class T>
void Alm<T>::GenFromCl(const TVector<T> & clin, const r_8 fwhm, RandomGeneratorInterface & rg)
{
  /*=======================================================================
     creates the a_lm from the power spectrum, 
     assuming they are gaussian  and complex 
     with a variance given by C(l)

     the input file should contain : l and C(l) with *consecutive* l's
     (missing C(l) are put to 0.)

     because the map is real we have : a_l-m = (-)^m conjug(a_lm)
     so we actually compute them only for m >= 0

=======================================================================*/
  int_4 nlmax= clin.NElts()-1;

  r_8 sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;
  int_4 n_l = nlmax+1;


  //    --- smoothes the initial power spectrum ---
  TVector<T> cl(clin, false);
  int l;
  for (l=0;l<n_l;l++)
    {
      r_8 gauss=exp(-l*(l+1.)*sig_smooth*sig_smooth);
      cl(l)*=(T)gauss;
    }
  

  //     --- generates randomly the alm according to their power spectrum ---
  r_8 hsqrt2 = 1.0 / Rac2;

  for (l=0;l<n_l;l++)
    {
      T rms=sqrt(cl(l));
      //        ------ m = 0 ------
      complex<T> zeta1((T)rg.Gaussian() );
      (*this)(l,0)   = zeta1 * rms;

      //------ m > 0 ------
      for (int m=1;m<=l;m++)
        {
          complex<T> aux1(hsqrt2);
          complex<T> aux2((T)rg.Gaussian() , (T)rg.Gaussian() );
          zeta1=aux1*aux2;
          (*this)(l,m)=rms*zeta1;
        }
    }
}


template <class T>
TVector<T> Alm<T>::powerSpectrum() const
{
  int_4 nlmax=Lmax();

  TVector<T>  powsp(nlmax+1);

  for (int l=0; l<=nlmax;l++)
    {
      powsp(l)=( (*this)(l,0) ).real()*( (*this)(l,0) ).real();
      for (int m=1; m<=l; m++)
        {
          powsp(l)+=2.*( (*this)(l,m).real()*(*this)(l,m).real()+
             (*this)(l,m).imag()*(*this)(l,m).imag() );
        }
      powsp(l)/=(2.*l+1.);
    }
  return powsp;
}

/*!
  \class Bm
  \ingroup Samba 
  Class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ 
  (then the size of the vector will be actually \f$2m_{max}+1)\f$.
  This class is used by the spherical harmonics transform server 
*/

 
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template Alm<r_8>
#pragma define_template Alm<r_4>
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class SOPHYA::Alm<r_8>;
template class SOPHYA::Alm<r_4>;
#endif