#include <stdlib.h>
#include <iostream>
#include <algorithm> //copy, sort
#include <iterator> // for ostream_iterator

#include "lagSphericTransform.h"
#include "lagsht_exceptions.h"

#include "walltimer.h"    //timing



#define DEBUG 0

namespace LagSHT {


  void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& flmn, 
					     std::vector<r_8>& fijk,
					     std::vector< std::complex<r_8> >& flmk) {
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
    
  fijk.resize(Ntot3DPix);

  //Perform the Inverse Laguerre Transform first

  flmk.resize(NtotCoeff);

#if DEBUG >= 1
  cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
  cout << "Start   Inverse Laguerre...." << endl;
#endif
  
  tstack_push("Laguerre MultiSynthesis");
  lagTrans_->MultiSynthesis(flmn, flmk, Nalm_);
  tstack_pop("Laguerre MultiSynthesis");


#if DEBUG >= 1
  cout << "End   Inverse Laguerre...." << endl;
  cout << "Start Inverse H.Spherical Transform " << endl;
#endif

  tstack_push("SHT Synthesis");
  int off7n, off7k;
  //perform the synthesis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 2
    std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
    sphere_->alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation

  }//loop on shell

    tstack_pop("SHT Synthesis");

#if DEBUG >= 1
    std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis



  void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk, 
					    std::vector< std::complex<r_8> >& flmn,
					    std::vector< std::complex<r_8> >& almk
					  ) {
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(fijk.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error: fijk");
  if(almk.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Analysis size error: almk");

  //done in Laguerre routine  flmn.resize(NtotCoeff);

#if DEBUG >= 1
  std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
  //Perform the Spherical Hamonic analysis first
    
  int off7n, off7k;

  tstack_push("SHT Analysis");

  //Ylm analysis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 1
    std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif

    sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
    
  }//loop on shell

#if DEBUG >= 1
  std::cout << "End   Harmonic Spheric..." << std::endl;
  std::cout << "Start Laguerre..." << std::endl;
#endif

  tstack_pop("SHT Analysis");
  tstack_push("Laguerre MultiAnalysis");

  //Laguerre Transform
  lagTrans_->MultiAnalysis(almk, flmn, Nalm_);

#if DEBUG >= 1
  std::cout << "End Laguerre..." << std::endl;
#endif

  tstack_pop("Laguerre MultiAnalysis");


#if DEBUG >= 2
  for (int i=0;i<NtotCoeff; i++) {
    std::cout << "flmn("<<i<<"): " <<flmn[i] <<  std::endl;
  }
#endif
#if DEBUG >= 1
  std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis



  void LaguerreSphericalTransform::Synthesis(const std::vector< std::complex<r_8> >& Elmn, 
					     const std::vector< std::complex<r_8> >& Blmn,
					     int spin,
					     std::vector<r_8>& fijk_re, std::vector<r_8>& fijk_im,
					     std::vector< std::complex<r_8> >& Elmk, 
					     std::vector< std::complex<r_8> >& Blmk
					   )   {
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(Elmn.size() != (size_t)NtotCoeff || Blmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
    
  fijk_re.resize(Ntot3DPix);
  fijk_im.resize(Ntot3DPix);

  //Perform the Inverse Laguerre Transform first

  Elmk.resize(NtotCoeff);
  Blmk.resize(NtotCoeff);

#if DEBUG >= 1
  std::cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << std::endl;
  std::cout << "Start   Inverse Laguerre...." << std::endl;
#endif
  
  tstack_push("Laguerre MultiSynthesis");
  lagTrans_->MultiSynthesis(Elmn, Elmk, Nalm_);
  lagTrans_->MultiSynthesis(Blmn, Blmk, Nalm_);
  tstack_pop("Laguerre MultiSynthesis");

#if DEBUG >= 1
  std::cout << "End   Inverse Laguerre...." << std::endl;
  std::cout << "Start Inverse H.Spherical Transform " << std::endl;
#endif

  tstack_push("SHT Synthesis");
  int off7n, off7k;
  //perform the synthesis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 2
    std::cout << "Synthesis... at shell num: " << n << std::endl;
#endif
    sphere_->alm2map_spin(reinterpret_cast<const r_8*>(&(Elmk.data()[off7n])),
			  reinterpret_cast<const r_8*>(&(Blmk.data()[off7n])),
			  &fijk_re.data()[off7k],
			  &fijk_im.data()[off7k],
			  spin,
			  false); //the false is to forbidd accumulation

  }//loop on shell

    tstack_pop("SHT Synthesis");

#if DEBUG >= 1
    std::cout << "SphLagTransform::Synthesis.... Done" << std::endl;
#endif
}//Synthesis


  void LaguerreSphericalTransform::Analysis(const std::vector<r_8>& fijk_re, 
					    const std::vector<r_8>& fijk_im,
					    int spin,
					    std::vector< std::complex<r_8> >& Elmn, 
					    std::vector< std::complex<r_8> >& Blmn,
					    std::vector< std::complex<r_8> >& Elmk, 
					    std::vector< std::complex<r_8> >& Blmk
					    ) {
  int Ntot3DPix = Npix_*N_;
  int NtotCoeff = Nalm_*N_;
  if(fijk_re.size() != (size_t)Ntot3DPix || fijk_im.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
//done in Laguerre routine   Elmn.resize(NtotCoeff);
//   Blmn.resize(NtotCoeff);

  //Perform the Inverse Laguerre Transform first

  Elmk.resize(NtotCoeff);
  Blmk.resize(NtotCoeff);


#if DEBUG >= 1
  std::cout << "SphLagTransform::Analysis.... Go..." << std::endl;
#endif
  //Perform the Spherical Hamonic analysis first
    
  int off7n, off7k;

  tstack_push("SHT Analysis");

  //Ylm analysis per shell
  for(int n =0; n<N_; n++){
    off7n = n*Nalm_;
    off7k = n*Npix_; 

#if DEBUG >= 1
    std::cout << "Start Harmonic Spheric... at n = "<< n << std::endl;
#endif

    sphere_->map2alm_spin(&(fijk_re.data()[off7k]), 
			  &(fijk_im.data()[off7k]),
			  reinterpret_cast<r_8*>(&(Elmk.data()[off7n])),
			  reinterpret_cast<r_8*>(&(Blmk.data()[off7n])),
			  spin,
			  false); //the false is to forbidd accumulation


  }//loop on shell

#if DEBUG >= 1
  std::cout << "End   Harmonic Spheric..." << std::endl;
  std::cout << "Start Laguerre..." << std::endl;
#endif

  tstack_pop("SHT Analysis");
  tstack_push("Laguerre MultiAnalysis");

  //Laguerre Transform
  lagTrans_->MultiAnalysis(Elmk, Elmn, Nalm_);
  lagTrans_->MultiAnalysis(Blmk, Blmn, Nalm_);

#if DEBUG >= 1
  std::cout << "End Laguerre..." << std::endl;
#endif

  tstack_pop("Laguerre MultiAnalysis");

#if DEBUG >= 2
  for (int i=0;i<NtotCoeff; i++) {
    std::cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] <<  std::endl;
  }
#endif
#if DEBUG >= 1
  std::cout << "SphLagTransform::Analysis.... Done" << std::endl;
#endif
}//Analysis

}// Fin de namespace