#ifndef LAMBDABUILDER_SEEN
#define LAMBDABUILDER_SEEN

#include <math.h>
#include "ndatablock.h"
#include "tvector.h" 
#include "alm.h"
#include "pexceptions.h"


namespace SOPHYA {

/*!  class of Legendre polynomials  */
class Legendre {
public :
  Legendre(r_8 x, int_4 lmax);

  inline r_8 getPl(int_4 l) const
  {
  if (l>lmax_)  {
    throw RangeCheckError("Legendre::getPl(l) l>lmax -> Instanciate class with greater lmax ...");
  }
  return Pl_(l);
  }

private : 
  Legendre() {;}   // this is forbidden

  void array_init(int_4 lmax);

  r_8 x_;
  int_4 lmax_;
  TVector<r_8> Pl_;
};



class LambdaLMBuilder {

public:

LambdaLMBuilder() {}
LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax);
LambdaLMBuilder(r_8 ct, r_8 st,int_4 lmax, int_4 mmax);

virtual ~LambdaLMBuilder() {};


//--- Optimisation Reza mai 2006 Serie de fonctions statiques 
// specialisees/optimisees pour SphericalTransform
static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
			   const Alm<r_8>& alm, Bm< complex<r_8> >& bm);
static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
			   const Alm<r_4>& alm, Bm< complex<r_4> >& bm);
static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
			      TVector< complex<r_8> >& phase, Alm<r_8> & alm);
static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
			      TVector< complex<r_4> >& phase, Alm<r_4> & alm);
//---- Fin foctions specialisees

 /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
inline double lamlm(int l, int m) const 
  {
    lambda_.CheckRelativeIndices(l,m);
    return lambda_(l,m); 
  }

inline double lamlm(int index) const {
 return lambda_(index); 
}


  //Return pointer to first element address of the alm's
//inline r_8* Data()  {return lambda_.Data();}


 private:
 static void updateArrayRecurrence(int_4 lmax);

 void array_init(int lmax, int mmax);


 // static TriangularMatrix<r_8>* a_recurrence_;
 static TriangularMatrix<r_8> a_recurrence_;
 TriangularMatrix<r_8> lambda_;

 protected :

 void updateArrayLamNorm();

 // static  TriangularMatrix<r_8>* lam_fact_;
 static  TriangularMatrix<r_8> lam_fact_;
 static  TVector<r_8>*  normal_l_;
 int_4 lmax_;
 int_4 mmax_;
 r_8 cth_;
 r_8 sth_;
 
 };


class LambdaWXBuilder : public LambdaLMBuilder
{
public:


LambdaWXBuilder() {}

LambdaWXBuilder(r_8 theta, int_4 lmax, int_4 mmax);
 
 /*! return the value of the coefficients \f$ _{w}\lambda_l^m\f$ and \f$_{x}\lambda_l^m\f$ */
inline void lam_wx(int l, int m, r_8& w, r_8& x) const
   {
  if (m > l ) 
    {
      throw RangeCheckError("LambdaWXBuilder::lam_wx : l < m !" );
    }
     w=lamWlm_(l,m);
     x=lamXlm_(l,m);
   }

 private:

 void array_init();


  TriangularMatrix<r_8> lamWlm_;
  TriangularMatrix<r_8> lamXlm_;
  

};

class LambdaPMBuilder : public LambdaLMBuilder
{
 public:

LambdaPMBuilder() {}

 LambdaPMBuilder(r_8 theta, int_4 lmax, int_4 mmax);
 /*! return the value of the coefficients \f$ _{+}\lambda_l^m\f$ and \f$_{-}\lambda_l^m\f$ */
inline void lam_pm(int l, int m, r_8& lambda_plus, r_8& lambda_moins) const
   {
  if (m > l ) 
    {
      throw RangeCheckError("LambdaPMBuilder::lam_pm : l < m !" );
    }
      lambda_plus = lamPlm_(l,m);
     lambda_moins = lamMlm_(l,m);
   }

 private:
 void array_init();

  TriangularMatrix<r_8> lamPlm_;
  TriangularMatrix<r_8> lamMlm_;

};

} // Fin du namespace

#endif