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

#include "lagSphericTransform.h"
#include "lagsht_exceptions.h"
//#include "lagsht_healpixhelper.h"

#include "walltimer.h"    //timing



#define DEBUG 0

namespace LagSHT {


void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn, 
					   vector<r_8>& fijk) {
  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
  vector< complex<r_8> > flmk(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
    cout << "Synthesis... at shell num: " << n << 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
  cout << "SphLagTransform::Synthesis.... Done" << endl;
#endif
}//Synthesis



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

#if DEBUG >= 1
  cout << "SphLagTransform::Analysis.... Go..." << endl;
#endif
  //Perform the Spherical Hamonic analysis first
    
  //  vector< complex<r_8> > flmk(NtotCoeff);

  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
    cout << "Start Harmonic Spheric... at n = "<< n << 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
  cout << "End   Harmonic Spheric..." << endl;
  cout << "Start Laguerre..." << endl;
#endif

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

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

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

  tstack_pop("Laguerre MultiAnalysis");


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



}// Fin de namespace