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

#include "alm.h"
template <class T>
Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm) 

{


  /*=======================================================================
     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;

  //alm.ReSize(nlmax);
  ReSizeRow(nlmax+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;
  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(NorRand());
      (*this)(l,0)   = zeta1 * rms;

      //------ m > 0 ------
      for (int m=1;m<=l;m++)
        {
          complex<T> aux1(hsqrt2);
          complex<T> aux2(NorRand(),NorRand());
          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;
}

#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 Alm<r_8>;
template class Alm<r_4>;
#endif