#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)
   {
     if (a_recurrence_.Size() == 0)
       {
	 //	 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_.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;
	   }
       }

   }