source: Sophya/trunk/SophyaLib/Samba/lambdaBuilder.cc@ 2959

#include "sopnamsp.h"
#include "lambdaBuilder.h"
#include "nbconst.h"


/*!
  \class SOPHYA::Legendre
generate Legendre polynomials : use in two steps : 

a) instanciate Legendre(\f$x\f$, \f$lmax\f$) ; \f$x\f$ is the value for wich Legendre polynomials will be required (usually equal to \f$\cos \theta\f$) and \f$lmax\f$ is the MAXIMUM value of the order of polynomials wich will be required in the following code (all polynomials, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula).

b) get the value of Legendre polynomial for a particular value of \f$l\f$ by calling the method getPl.

*/


Legendre::Legendre(r_8 x, int_4 lmax) 
{
  if (fabs(x) > 1. ) 
    {
      throw RangeCheckError("variable for Legendre polynomials must have modules inferior to 1" );
    }
  x_ = x;
  array_init(lmax);
}
   
/*! \fn void SOPHYA::Legendre::array_init(int_4 lmax)

compute all \f$P_l(x,l_{max})\f$ for \f$l=1,l_{max}\f$ 
*/
void Legendre::array_init(int_4 lmax)
{
  lmax_ = lmax;
  Pl_.ReSize(lmax_+1);
  Pl_(0)=1.;
  if (lmax>0) Pl_(1)=x_;
  for (int k=2; k<Pl_.NElts(); k++)
    {
      Pl_(k) = ( (2.*k-1)*x_*Pl_(k-1)-(k-1)*Pl_(k-2) )/k;
    }
}

TriangularMatrix<r_8> LambdaLMBuilder::a_recurrence_ = TriangularMatrix<r_8>();
TriangularMatrix<r_8> LambdaLMBuilder::lam_fact_     = TriangularMatrix<r_8>();
TVector<r_8>* LambdaLMBuilder::normal_l_     = NULL;



/*! \class SOPHYA::LambdaLMBuilder


This class generate the coefficients :
\f[
            \lambda_l^m=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}
            P_l^m(\cos{\theta})
\f]
where \f$P_l^m\f$ are the associated Legendre polynomials.  The above coefficients contain the theta-dependance of spheric harmonics : 
\f[
            Y_{lm}(\cos{\theta})=\lambda_l^m(\cos{\theta}) e^{im\phi}.
\f]

Each object has a fixed theta (radians), and maximum l and m to be calculated
(lmax and mmax).
 use the class in two steps : 
a) instanciate  LambdaLMBuilder(\f$\theta\f$, \f$lmax\f$, \f$mmax\f$) ;  \f$lmax\f$ and \f$mmax\f$ are  MAXIMUM values for which \f$\lambda_l^m\f$ will be required in the following code (all coefficients, from \f$l=0 to lmax\f$, are computed once for all by an iterative formula).
b) get the values of coefficients for  particular values of \f$l\f$ and \f$m\f$ by calling the method lamlm.
*/


LambdaLMBuilder::LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax)
    {
      cth_=cos(theta);
      sth_=sin(theta);
      array_init(lmax, mmax);
    }
LambdaLMBuilder::LambdaLMBuilder(r_8 costet, r_8 sintet,int_4 lmax, int_4 mmax)
    {
      cth_=costet;
      sth_=sintet;
      array_init(lmax, mmax);
    }
void LambdaLMBuilder::array_init(int lmax, int mmax)
   {
     updateArrayRecurrence(lmax);

     lmax_=lmax;
     mmax_=mmax;
     r_8 bignorm2 = 1.e268; // = 1e-20*1.d288

     lambda_.ReSizeRow(lmax_+1);

     r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
          
     for (int m=0; m<=mmax_;m++)
       {


         lambda_(m,m)= lam_mm / bignorm2; 

         r_8 lam_0=0.;
         r_8 lam_1=1. /bignorm2 ;
         //      r_8 a_rec = LWK->a_recurr(m,m);
         r_8 a_rec = a_recurrence_(m,m);
         r_8 b_rec = 0.;
         for (int l=m+1; l<=lmax_; l++)
         {
           r_8 lam_2 = (cth_*lam_1-b_rec*lam_0)*a_rec;
           lambda_(l,m) = lam_2*lam_mm;
           b_rec=1./a_rec;
           //      a_rec= LWK->a_recurr(l,m);
           a_rec= a_recurrence_(l,m);
           lam_0 = lam_1;
           lam_1 = lam_2;
         }

         lam_mm = -lam_mm*sth_* sqrt( (2.*m+3.)/ (2.*m+2.) );

       }
   }

/*!
  \brief : Specialized/optimized static function for fast spherical harmonic transform.
  Computes bm(m) = Sum_l>=m [ lambda(l,m) * alm(l,m) ] 
  See SphericalTransformServer<T>::GenerateFromAlm(map, pixsize, alm)
*/
void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
                                     const Alm<r_8>& alm, Bm< complex<r_8> >& bm)
{
  updateArrayRecurrence(lmax);
  r_8 cth = cos(theta);
  r_8 sth = sin(theta);
  
  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
  
  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
  r_8 lam_0, lam_1, lam_2;
  register r_8 a_rec, b_rec;

  const complex<r_8>* almp = alm.columnData(0); 
  r_8* arecp = a_recurrence_.columnData(0); 

  int_4 l, m;

  for (m=0; m<=mmax;m++)   {
    //MOV^     const complex<r_8>* almp = alm.columnData(m); 
    complex<r_8>* bmp = &(bm(m));
    
    *bmp = (lam_mm / bignorm2)*(*almp);  almp++;

    lam_0=0.;
    lam_1=1. /bignorm2 ;
    
    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
    b_rec = 0.;
    for (l=m+1; l<=lmax; l++) {
      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
      
      //DEL        lambda_(l,m) = lam_2*lam_mm;
      *bmp += lam_2*lam_mm*(*almp);   almp++;
      b_rec=1./a_rec;
      //           a_rec= LWK->a_recurr(l,m);
      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
      lam_0 = lam_1;
      lam_1 = lam_2;
    }
    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
  }
}

/*!
  \brief : Specialized/optimized static function for fast spherical harmonic transform.
*/
void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax,  
                                     const Alm<r_4>& alm, Bm< complex<r_4> >& bm)
{
  Alm<r_8> alm8(alm.Lmax());
  for(sa_size_t k=0; k<alm8.Size(); k++) 
    alm8(k) = complex<r_8>((r_8)alm(k).real() , (r_8)alm(k).imag());
  Bm< complex<r_8> > bm8(bm.Mmax());
  ComputeBmFrAlm(theta, lmax, mmax, alm8, bm8);
  for(sa_size_t kk=-bm.Mmax(); kk<=bm.Mmax(); kk++) 
    bm(kk)= complex<r_4>((r_4)bm8(kk).real() , (r_4)bm8(kk).imag());
  return;
}


/*!
  \brief : Specialized/optimized static function for fast spherical harmonic transform.
  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ] 
  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm) 
*/
void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
                                        TVector< complex<r_8> >& phase, Alm<r_8> & alm)
{
  updateArrayRecurrence(lmax);
  r_8 cth = cos(theta);
  r_8 sth = sin(theta);
  
  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
  
  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
  r_8 lam_0, lam_1, lam_2;
  register r_8 a_rec, b_rec;
  
  complex<r_8>* almp = alm.columnData(0); 
  r_8* arecp = a_recurrence_.columnData(0); 

  int_4 l, m;
  for (m=0; m<=mmax;m++)  {
    //MOV^ complex<r_8>* almp = alm.columnData(m); 
    complex<r_8>  phi = phase(m);
    
    *almp += (lam_mm / bignorm2)*phi;    almp++;

    lam_0=0.;
    lam_1=1. /bignorm2 ;
    
    //MOV^ r_8* arecp = a_recurrence_.columnData(m); 
    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
    b_rec = 0.;
    for (l=m+1; l<=lmax; l++)   {
      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
      
      //DEL        lambda_(l,m) = lam_2*lam_mm;
      *almp += (lam_2*lam_mm) * phi;    almp++;
      
      b_rec=1./a_rec;
      //           a_rec= LWK->a_recurr(l,m);
      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
      lam_0 = lam_1;
      lam_1 = lam_2;
    }
    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
  }
}

/*!
  \brief : Specialized/optimized static function for fast spherical harmonic transform.
  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ] 
  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm) 
*/
void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
                                        TVector< complex<r_4> >& phase, Alm<r_4> & alm) 
{
  updateArrayRecurrence(lmax);
  r_8 cth = cos(theta);
  r_8 sth = sin(theta);
  
  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
  
  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
  r_8 lam_0, lam_1, lam_2;
  register r_8 a_rec, b_rec;
  register r_4 cff;
  complex<r_4>* almp = alm.columnData(0); 
  r_8* arecp = a_recurrence_.columnData(0); 

  int_4 l, m;
  for (m=0; m<=mmax;m++)  {
    //MOV^ complex<r_4>* almp = alm.columnData(m); 
    complex<r_4>  phi = phase(m);
    
    cff = lam_mm / bignorm2;
    *almp += cff*phi;  almp++;

    lam_0=0.;
    lam_1=1. /bignorm2 ;
    
    //MOV^ r_8* arecp = a_recurrence_.columnData(m); 
    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
    b_rec = 0.;
    for (l=m+1; l<=lmax; l++)   {
      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
      
      //DEL        lambda_(l,m) = lam_2*lam_mm;
      cff = (lam_2*lam_mm);
      *almp += cff * phi;    almp++;
      
      b_rec=1./a_rec;
      //           a_rec= LWK->a_recurr(l,m);
      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
      lam_0 = lam_1;
      lam_1 = lam_2;
    }
    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
  }
}


/*! \fn void  SOPHYA::LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)

 compute a static array of coefficients independant from theta (common to all instances of the LambdaBuilder Class 
*/
void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
   {
     if ( (a_recurrence_.Size() > 0) && (lmax < a_recurrence_.rowNumber()) ) 
       return;  // Pas besoin de recalculer le tableau de recurrence
     if (a_recurrence_.Size() > 0)  {
       cout << " WARNING : The classes LambdaXXBuilder will be more efficient "
            << "if instanciated with parameter lmax = maximum value of l index "
            << "which will be needed in the whole application (arrays not "
            << "recomputed) " << endl;
       cout << "   lmax= " << lmax << " previous instanciation with lmax= " 
            <<  a_recurrence_.rowNumber()-1 << endl;
         }

     a_recurrence_.ReSizeRow(lmax+1);
     for (int m=0; m<=lmax;m++)
     {
       a_recurrence_(m,m) = sqrt( 2.*m +3.);
       for (int l=m+1; l<=lmax; l++)
         {
           r_8 fl2 = (l+1.)*(l+1.);
           a_recurrence_(l,m)=sqrt( (4.*fl2-1.)/(fl2-m*m) );
         }
     }
   }

/*! \fn void  SOPHYA::LambdaLMBuilder::updateArrayLamNorm()

 compute  static arrays of coefficients independant from theta (common to all instances of the derived  classes 
*/
void  LambdaLMBuilder::updateArrayLamNorm()
     {
       lam_fact_.ReSizeRow(lmax_+1);
       for(int m = 0;m<= lmax_; m++)
         {
           for (int l=m; l<=lmax_; l++)
             {
               lam_fact_(l,m) =2.*(r_8)sqrt( (2.*l+1)*(l+m)*(l-m)/(2.*l-1)  );
             }
         }
       (*normal_l_).ReSize(lmax_+1);
       (*normal_l_)(0)=0.;
       (*normal_l_)(1)=0.;
       for (int l=2; l< (*normal_l_).NElts(); l++)
         {
           (*normal_l_)(l) =(r_8)sqrt( 2./( (l+2)*(l+1)*l*(l-1) ) );    
         }
     }


/*! \class SOPHYA::LambdaWXBuilder

This class generates the coefficients :
\f[
            _{w}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}} G^+_{lm}
\f]
\f[
            _{x}\lambda_l^m=-2\sqrt{\frac{2(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}G^-_{lm}
\f]
where
\f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})
\f]
and
\f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)
\f]
 \f$P_l^m\f$ are the associated Legendre polynomials.

The coefficients express the theta-dependance of the \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ functions :
\f[W_{lm}(\cos{\theta}) = \sqrt{\frac{(l+2)!}{2(l-2)!}}_w\lambda_l^m(\cos{\theta})e^{im\phi}
\f]
\f[X_{lm}(\cos{\theta}) = -i\sqrt{\frac{(l+2)!}{2(l-2)!}}_x\lambda_l^m(\cos{\theta})e^{im\phi}
\f]
 where \f$W_{lm}(\cos{\theta})\f$ and \f$X_{lm}(\cos{\theta})\f$ are defined as :

\f[
W_{lm}(\cos{\theta})=-\frac{1}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(
_{+2}Y_l^m(\cos{\theta})+_{-2}Y_l^m(\cos{\theta})\right)
\f]
\f[X_{lm}(\cos{\theta})=-\frac{i}{2}\sqrt{\frac{(l+2)!}{(l-2)!}}\left(
_{+2}Y_l^m(\cos{\theta})-_{-2}Y_l^m(\cos{\theta})\right)
\f]

*/


LambdaWXBuilder::LambdaWXBuilder(r_8 theta, int_4 lmax, int_4 mmax) : LambdaLMBuilder(theta, lmax, mmax)
  {
    array_init();
  }


void LambdaWXBuilder::array_init()
   {
     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
       {
         //      lam_fact_ = new  TriangularMatrix<r_8>;
         normal_l_ = new TVector<r_8>;
         updateArrayLamNorm();
       }
     else
       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
         {
           updateArrayLamNorm();
         }

     r_8 one_on_s2  = 1. / (sth_*sth_) ;    // 1/sin^2
     r_8 c_on_s2    = cth_*one_on_s2;
     lamWlm_.ReSizeRow(lmax_+1);
     lamXlm_.ReSizeRow(lmax_+1);

     // calcul des lambda
     for(int m = 0;m<= mmax_; m++)
       {
         for (int l=m; l<=lmax_; l++)
           {
             lamWlm_(l,m) =  0.;
             lamXlm_(l,m) =  0.;                
           }
       }
     for(int l = 2;l<= lmax_; l++)
       {
         r_8 normal_l =  (*normal_l_)(l);       
         for (int m=0; m<=l; m++)
        {
          r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
          r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
          r_8 lam_fact_l_m = lam_fact_(l,m);
          r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
          r_8  b_w =  c_on_s2 * lam_fact_l_m;
          r_8  a_x =  2. * cth_ * (l-1.);
          lamWlm_(l,m) =   normal_l * ( a_w * lam_lm - b_w * lam_lm1m );
          lamXlm_(l,m) = - normal_l * m* one_on_s2* ( a_x * lam_lm - lam_fact_l_m * lam_lm1m );
        }
       }

   }

/*!   \class SOPHYA::LambdaPMBuilder

This class generates the coefficients
\f[
            _{\pm}\lambda_l^m=2\sqrt{\frac{(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\left( G^+_{lm} \mp G^-_{lm}\right)
\f]
where
\f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})
\f]
and
\f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)
\f]
and \f$P_l^m\f$ are the associated Legendre polynomials.
The coefficients express the theta-dependance of the  spin-2 spherical harmonics :
\f[_{\pm2}Y_l^m(\cos{\theta})=_\pm\lambda_l^m(\cos{\theta})e^{im\phi}
\f]
*/

LambdaPMBuilder::LambdaPMBuilder(r_8 theta, int_4 lmax, int_4 mmax) : LambdaLMBuilder(theta, lmax, mmax)
  {
    array_init();
  }


void LambdaPMBuilder::array_init()
   {
     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
       {
         //      lam_fact_ = new  TriangularMatrix<r_8>;
         normal_l_ = new TVector<r_8>;
         updateArrayLamNorm();
       }
     else
       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
         {
           updateArrayLamNorm();
         }

     r_8 one_on_s2 = 1. / (sth_*sth_) ;
     r_8 c_on_s2 = cth_*one_on_s2;
     lamPlm_.ReSizeRow(lmax_+1);
     lamMlm_.ReSizeRow(lmax_+1);

     // calcul des lambda
     for(int m = 0;m<= mmax_; m++)
       {
         for (int l=m; l<=lmax_; l++)
           {
             lamPlm_(l,m) =  0.;
             lamMlm_(l,m) =  0.;                
           }
       }

     for(int l = 2;l<= lmax_; l++)
       {
         r_8 normal_l =  (*normal_l_)(l);       
         for (int m=0; m<=l; m++)
           {
             r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
             r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
             r_8 lam_fact_l_m = lam_fact_(l,m);
             r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
             r_8 f_w =  lam_fact_l_m/(sth_*sth_);
             r_8  c_w =  2*m*(l-1.) * c_on_s2;

             lamPlm_(l,m)  = normal_l * ( -(a_w+c_w) * lam_lm + f_w*( cth_ + m) * lam_lm1m )/Rac2;
             lamMlm_(l,m) = normal_l * ( -(a_w-c_w) * lam_lm + f_w*( cth_ - m) * lam_lm1m )/Rac2;
           }
       }

   }