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

#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 (abs(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.;
  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_ = NULL;
TriangularMatrix<r_8>* LambdaLMBuilder::lam_fact_     = NULL;
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);
    }
void LambdaLMBuilder::array_init(int lmax, int mmax)
   {
     if (a_recurrence_ == NULL)
       {
         a_recurrence_ = new TriangularMatrix<r_8>;
         updateArrayRecurrence(lmax);
       }
     else
       if ( lmax > (*a_recurrence_).rowNumber()-1     )
         {
           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() << endl;
         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.) );

       }
   }


/*! \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)
   {
     (*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_ == NULL || 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_ == NULL || 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;
           }
       }

   }