#include <math.h>
#include "anagen.h"
#include "lambuilder.h"
#include "fftserver.h"

/* assumes the map is symmetric about the equator and that
each strip is divided equally*/
 
extern "C" {
  void fft_gpd_(long double* ,int& ,int& ,int& ,int& ,long double*);
	   }

void map2almgen(int nsmax,int nlmax,int nmmax,const SphericalMap<double>& map,
	     vector< vector< complex<double> > >& alm, double cos_theta_cut){
  
  //       REAL*4  map(12*nsmax**2)               ! 4*12*nsmax**2
  
  //       REAL*4 powspec(0:nlmax)
    
  //       integer npmiss,npmt,id_miss(10000)
  
  alm.resize(nlmax+1);
  for (int i=0; i< (signed) alm.size();i++)
    {
      alm[i].resize(nmmax+1);
      for (int j=0; j< (signed) alm[i].size();j++)alm[i][j]=0;
    }
  
  /*-----------------------------------------------------------------------
    computes the integral in phi : phas_m(theta)
    for each parallele from north to south pole
    -----------------------------------------------------------------------*/
  
  vector< complex<double> > phas_n(nmmax+1), phas_s(nmmax+1);
  for (int ith = 0; ith <= (map.NbThetaSlices()+1)/2-1;ith++){
    int nph;
    double phi0;
    double theta;
    Vector phin, datan, datas;
    TVector<double> phis;
    map.GetThetaSlice(map.NbThetaSlices()-ith-1,theta,phis,datas);
    map.GetThetaSlice(ith,theta,phin,datan);
    for (int i=0;i< nmmax+1;i++){
      phas_n[i]=0; phas_s[i]=0;
    }
    nph = phin.NElts();
    phi0 = phin(0); 
    double cth = cos(theta);

    //part of the sky out of the symetric cut
    bool keep_it = (abs(cth) >= cos_theta_cut); 

    //make sure that map is well defined
    if (keep_it){
      comp_phas2gen(nsmax,nlmax,nmmax,datan,nph,phas_n,phi0);
    }

    if (ith != map.NbThetaSlices()/2 && keep_it){
      comp_phas2gen(nsmax,nlmax,nmmax,datas,nph,phas_s,phi0);
    }
    /*-----------------------------------------------------------------------
      computes the a_lm by integrating over theta
      lambda_lm(theta) * phas_m(theta)
      for each m and l
      -----------------------------------------------------------------------*/
    LambdaBuilder lb(acos(cth),nlmax,nmmax);
    double domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
    for (int m = 0; m <= nmmax; m++)
      {
	alm[m][m] += (lb.lamlm(m,m) * (phas_n[m] + phas_s[m])) * domega; //m,m even
	for (int l = m+1; l<= nlmax; l++)
	  {
	    alm[l][m] += (lb.lamlm(l,m) * phas_n[m] + lb.lamlm(l,m,-1)*phas_s[m])*domega; //assuming the map is symmetric about the equator
	  }
      }
  }
}


void comp_phas2gen(int nsmax,int nlmax,int nmmax, Vector datain,
	       	int nph,vector< complex<double> >& dataout,double phi0){
  /*=======================================================================
    integrates (data * phi-dependence-of-Ylm) over phi
    --> function of m can be computed by FFT
    with  0<= m <= npoints/2 (: Nyquist)
    because the data is real the negative m are the conjugate of the 
    positive ones
    
    arguments d'appels : GLM
    =======================================================================*/

  //FFTVector inf(datain);
  //FFTVector outf=FFTServer::solve(inf,-1);
  FFTServer fft;
  Vector outf;
  //cout<<"in :"<<datain<<endl;
  fft.fftb(datain,outf);
  //  cout<<outf<<endl;
  long double * data = new long double[nph*2];
  //outf.d((double*)data, nph);
  data[0]=outf(0); data[1]=0;
  for (int i=2;i<=nph;i++) data[i]=outf(i-1);
  for (int i=nph+1;i<2*nph;i++) data[i]=0;

  int ksign = -1;
//   long double data[nph*2];

//   for (int i = 0; i< nph;i++){
//     data[2*i] = datain(i);
//     data[2*i+1]=0;
//   }

//   long double work[nph*2];
//   int dum1=1;
//   int dum0=0;
//   fft_gpd_(data,nph,dum1,ksign,dum0,work); /* real to complex 
// 				      for any nph (not necessary power of 2)*/
//     //in the output the frequencies are respectively 0,1,2,..,nph/2,-nph/2+1,..,-2,-1
//     //     only the first nph/2+1 (positive freq.) are interesting
//   for (int i=0;i<nph*2;i++){cout << "heh"<<test[i]<<" "<<data[i]<<endl;}*/

  int im_max = min(nph/2,nmmax);
  dataout.resize(nmmax+1);
  for (int i = 1;i <= im_max + 1;i++){
    int  m = ksign*(i-1);
    complex<double> shit(data[2*(i-1)],data[2*(i-1)+1]);
    complex<double> fuck(cos(m*phi0),sin(m*phi0));
    dataout[i-1]=shit*fuck;
  }
  for (int i = im_max + 2;i <= nmmax + 1;i++){
         dataout[i-1] = 0;
  }
  delete data;
}