#ifndef LAGSPHERICTRANSFORM_SEEN
#define LAGSPHERICTRANSFORM_SEEN

// #include <utility>      // std::pair, std::make_pair
// using namespace std;

//Sharp
// #include "sharp_lowlevel.h"
// #include "sharp_geomhelpers.h"
// #include "sharp_almhelpers.h"
// #include "sharp_cxx.h"

#include "lagsht_spheregeom.h"
#include "laguerreTransform.h"
//#include "lagsht_geom.h"
#include "lagsht_exceptions.h"


namespace LagSHT {


////////////////////////////////////////////////////////////////
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////

class LaguerreSphericalTransform {
 public:


  LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax, 
			     r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2): 
    R_(R) {
    
    //Factory
    if (spheregeom == "ECP" || spheregeom == "Fejer1") {
      sphere_ = new ECPGeometry(Lmax, Mmax, Nrings, Nphi);
    } else if (spheregeom == "Gauss") {
      sphere_ = new GaussGeometry(Lmax, Mmax, Nrings, Nphi);      
    } else if (spheregeom == "Healpix") {
      sphere_ = new HealpixGeometry(Lmax, Mmax, Nrings); //here Nsides = Nrings
    } else {
      cout << "FAILURE: LaguerreSphericalTransform::Ctor wrong geometry type" << endl;
      sphere_ = 0;
      throw LagSHTError("LaguerreSphericalTransform::Ctor wrong geometry type");
    }

    lagTrans_ = new LaguerreTransform(Nmax, R, alpha); 

    
    //do the map pixelization
    if(sphere_) {
      //store localy usefull parameters for the algorithms Analysis/Synthesis
      Npix_ = sphere_->Npix();
      Nalm_ = sphere_->Nalm();
    } else {
      throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
    }


    if(lagTrans_) {
      N_    = lagTrans_->GetOrder();
      alpha_= lagTrans_->GetAlpha();
    } else {
      throw LagSHTError("LaguerreSphericalTransform::Ctor lagtrans_ ptr NULL");
    }


  }//Ctor

    
  //! Destructor
  virtual ~LaguerreSphericalTransform() {
    if(sphere_) delete sphere_;  sphere_ = 0;
    if(lagTrans_) delete lagTrans_; lagTrans_ = 0;
  }
  
  //! Geometry helper
  BaseGeometry* GetSphericalGeometry() const { return sphere_;}

  //! Laguerre Transform helper
  LaguerreTransform* GetLagTransform() const { return lagTrans_;}

   //! Getters
  int GetOrder() { return N_; }
  int GetAlpha() { return alpha_;}


 //! Synthesis 
 /*! \brief Coeffs -> Pixels with Input/Output using T floating representation
   \input flmn: 3D complex spherical-laguerre coefficients
   \output fijk: 3D real function values
 */
  void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk);

  //! Analysis
  /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
    \input fijk: 3D real function values
    \output flmn: 3D complex spherical-laguerre coefficients
  */
  void Analysis(const vector<r_8>& fijk, vector< complex<r_8> >& flmn){
    int NtotCoeff = Nalm_*N_;
    vector< complex<r_8> > almk(NtotCoeff);
    Analysis(fijk,flmn,almk);
  }



  //! Analysis
  /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
    \input fijk: 3D real function values
    \output flmn: 3D complex spherical-laguerre coefficients
    \output almk: 2D alm complex spherical coefficients for each subshell k
  */
  void Analysis(const vector<r_8>& fijk, 
		vector< complex<r_8> >& flmn,
		vector< complex<r_8> >& almk
		);


protected:
  r_8 R_; //!< Radial maximal value (not yet used)

  BaseGeometry* sphere_; //<! the 2D pixelization of the shells 

  LaguerreTransform* lagTrans_; //!< the Laguerre Transform
  
  int Npix_; //!< Total number of 2D pixels
  int Nalm_; //!< Total number of Alm coefficients

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

}; //class


}//end namespace
#endif