#ifndef LAGUERRETRANSFORM_SEEN
#define LAGUERRETRANSFORM_SEEN
/*
 *  This file is part of LagSHT.
 *
 *  LagSHT is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  LagSHT is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libsharp; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 *  LagSHT is being developed at the Linear Accelerateur Laboratory (LAL)
 *  91898 ORSAY CEDEX - FRANCE
 *  main author: J.E Campagne
 */


#include <complex>

#include "laguerreBuilder.h"


namespace LagSHT {
class BaseLaguerreTransform {
 public:
  //! Creator
  BaseLaguerreTransform(int n, int alpha=2): N_(n), alpha_(alpha) {}
  //! Destructor
  virtual ~BaseLaguerreTransform() {}
  //! Getters
  int GetOrder() { return N_; }
  int GetAlpha() { return alpha_;}

protected:
  int N_; //!< Order of the Generalized Lagurerre polynomial
  int alpha_;  //!< second parameter of the generalized Laguerre polynomial 

};//BaseLaguerreTransform


///////////////////////////////////////////////////////////////////////
//// ------------------- Class LaguerreTransform   ----------------- //
//// Perform single or multiple 1D-Laguerre Transform and Inverse    //
//// using generalized Laguerre Polynomials except for T:double      //
//// where in this case the generalized Laguerre functions are used  //
///////////////////////////////////////////////////////////////////////
//! Laguerre Transform Syntheis/Analysis
class LaguerreTransform : public BaseLaguerreTransform {
 public:
  //! Creator
  LaguerreTransform(int N, r_8 R, int alpha=2, bool R2tau=true);
  //! Destructor
  virtual ~LaguerreTransform() {}

  //! Getters
  vector<r_8> GetNodes()   const { return nodes_; }
  vector<r_8> GetWeights() const { return weights_; }
  r_8         GetTau()     const { return tau_; }

  //! Find the closest nodes: no boundchecking
  int R2Index(r_8 r) const;
  //! return the nodes position: no boundchecking
  r_8 Index2R(int k) const { 
    //    return nodes_[k]*R_/nodes_[N_-1];
    //JEC 22/9/15
    return nodes_[k]*tau_;
  }
  
  


  //! Multi Analysis at once (Values -> Coefficients)
  /*!
    \param fi: values f(r_i) with r_i zero L^(alpha)_N
    \param fn : result
    \param stride : step between 2 conscecutive bunches of fi & fn 
  */
  void MultiAnalysis(const vector< complex<r_8> >& fi, vector< complex<r_8> >& fn, int stride);



  //! Multi Analysis at once (Values -> Coefficients)
  /*!
    \param fi: values f(r_i) with r_i zero L^(alpha)_N (Real case)
    \param fn : result
    \param stride : step between 2 conscecutive bunches of fi & fn 
  */
  void MultiAnalysisR(const vector<r_8>& fi, vector<r_8>& fn, int stride);



  //! Multiple Synthesis at once (Coefficientss -> Values)
  /*!
  \param fn: (input) complex coefficients of Laguerre Transform
  \param stride: (input) is the size of the non-radial dimension (eg. Y"lm" index range size = number of alm)
  \param fi: (output) complex values of the Inverse Laguerre Transform
  */
  void MultiSynthesis(const vector< complex<r_8> >& fn, vector< complex<r_8> >& fi, int stride);

  //! Single Analysis (Values -> Coefs)
  /*!
    \param fi: values f(r_i) with r_i zero L^(alpha)_N
    \param fn : result
  */
  void Analysis(const vector< complex<r_8> >& fi, vector< complex<r_8> >& fn) {
    MultiAnalysis(fi,fn,1);
  }//Analysis

  //! Single Inverse Laguerre Transform (Coefficientss -> Values)
  void Synthesis(const vector< complex<r_8> >& fn, vector< complex<r_8> >& fi) {
    MultiSynthesis(fn,fi,1);
  }//Synthesis

protected:

  r_8 alphaFact_; //!< sqrt(alpha!)
  r_8 R_;         //!< maximal radius

  vector<r_8> nodes_;    //!< nodes   of the Gauss-Laguerre quadrature used
  vector<r_8> weights_;  //!< weights of the Gauss-Laguerre quadrature used

  r_8 tau_;              //!< a scaling factor = R/nodes[N-1] JEC 22/9/15
  r_8 tau3sqrt_;         // tau^(3/2)
  r_8 tau3sqrtinv_;      // tau^(-3/2)
 public:
  class distance {
  public: 
    distance(r_8 r0): r0_(r0) {}
    r_8 operator()(r_8 x){ r_8 xn=x-r0_; return xn*xn;}
  private:
    r_8 r0_;
  };//distance

};//LaguerreTransform



}//end namespace
#endif