#ifndef LAGSPHERICTRANSFORM_SEEN
#define LAGSPHERICTRANSFORM_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 "lagsht_spheregeom.h"
#include "laguerreTransform.h"
#include "lagsht_exceptions.h"


namespace LagSHT {


////////////////////////////////////////////////////////////////
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////
/* Example to access the Laguerre Coeff flmn and the intermediate alm(r_k) Coeff.
  for(int n=0;n<Nmax;n++){
    for(int l=0;l<Lmax;l++){
      for (int m=0;m<=l;m++) {
	int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
      }
    }
  }
 */
class LaguerreSphericalTransform {
 public:


  LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax, 
			     r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2): 
    Lmax_(Lmax), Mmax_(Mmax), 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()   const { return N_; }
  int GetAlpha()   const { return alpha_; }
  int GetNPixels() const { return Npix_; }
  int GetNAlm()    const { return Nalm_; }


  //! Get the f(l,m,n) coefficient in the vector collection
  //make a Flmn Class
  inline complex<r_8> GetFlmn(const vector< complex<r_8> >& flmn, int l , int m, int n) const {
    int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
    return flmn[id];
 }
  //! Set the f(l,m,n) coefficient in the vector collection
  inline void SetFlmn(vector< complex<r_8> >& flmn, int l , int m, int n, complex<r_8> value) const {
     int id= n*Nalm_ + l +m*Lmax_ -m*(m+1)/2;
     flmn[id] = value;
  }
  


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

  //backward compatibility
  void Synthesis(const vector< complex<r_8> >& flmn, vector<r_8>& fijk) {
    int NtotCoeff = Nalm_*N_;
    vector< complex<r_8> > almk(NtotCoeff);
    Synthesis(flmn,fijk,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
  */
  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:

  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 Lmax_; //!< Spherical harmonic limit such that l:0,...,Lmax-1   m:0,...,l
  int Mmax_; //!< maximum value of m


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

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


}; //class


}//end namespace
#endif