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

extern "C" {
  void fft_gpd_(long double* ,int& ,int& ,int& ,int& ,long double*);
  float gasdev2_(int& idum);
	   }

void alm2mapgen(const int nsmax,const int nlmax,const int nmmax,
	     const vector< vector< complex<double> > >&alm,
	     SphericalMap<double>& map,
	     vector< complex<long double> > b_north,
	     vector< complex<long double> > b_south,
	     vector< complex<long double> > bw,
	     vector< complex<long double> > data,
	     vector< complex<long double> > work){
  /*=======================================================================
    computes a map form its alm for the HEALPIX pixelisation
    map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
    = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
    
    where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
    
    * the recurrence of Ylm is the standard one (cf Num Rec)
    * the sum over m is done by FFT
    
    =======================================================================*/
  
  //  map.resize(12*nsmax*nsmax);                       
  b_north.resize(2*nmmax+1); //index m corresponds to nmmax+m
  b_south.resize(2*nmmax+1);
  bw.resize(4*nsmax);
  data.resize(4*nsmax);
  work.resize(4*nsmax);


  /*      INTEGER l, indl */

  for (int ith = 0; ith <= (map.NbThetaSlices()+1)/2-1;ith++){
    int nph;
    double phi0;
    double theta,thetas;
    Vector phin, datan,datas;
    TVector<double> phis;
    map.GetThetaSlice(map.NbThetaSlices()-ith-1,thetas,phis,datas);
    map.GetThetaSlice(ith,theta,phin,datan);
    nph = phin.NElts();
    phi0 = phin(0); 
    double cth = cos(theta);

    /*        -----------------------------------------------------
	      for each theta, and each m, computes
	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
	      ------------------------------------------------------
	      lambda_mm tends to go down when m increases (risk of underflow)
	      lambda_lm tends to go up   when l increases (risk of overflow)*/
     // lambda_0_0 * big number --> to avoid underflow
    LambdaBuilder lb(acos(cth),nlmax,nmmax);
    for (int m = 0; m <= nmmax; m++){
      complex<double> b_n = lb.lamlm(m,m) * alm[m][m];
      complex<double>  b_s = lb.lamlm(m,m) * alm[m][m]; //m,m even
      for (int l = m+1; l<= nlmax; l++){
	b_n += lb.lamlm(l,m)*alm[l][m];
	b_s += lb.lamlm(l,m,-1)*alm[l][m];
      }
      b_north[m+nmmax] = b_n;
      b_south[m+nmmax] = b_s;
    }

    //        obtains the negative m of b(m,theta) (= complex conjugate)

    for (int m=-nmmax;m<=-1;m++){
      //compiler doesn't have conj()
      complex<long double> 
	shit(b_north[-m+nmmax].real(),-b_north[-m+nmmax].imag());
      complex<long double> 
	shit2(b_south[-m+nmmax].real(),-b_south[-m+nmmax].imag());
      b_north[m+nmmax] = shit;
      b_south[m+nmmax] = shit2;
    }
    /*        ---------------------------------------------------------------
	      sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
        ---------------------------------------------------------------*/
    /*cout <<"gen ";
    for (int i=0;i<2*nmmax+1;i++){
      cout << b_north[i]<<" ";}
    cout<<endl;*/

    syn_phasgen(nsmax,nlmax,nmmax,b_north,nph,datan,phi0,bw); // north hemisph. + equator
    if (ith != map.NbThetaSlices()/2){
      syn_phasgen(nsmax,nlmax,nmmax,b_south,nph,datas,phi0,bw); // south hemisph. w/o equat
      for (int i=0;i< nph;i++){
	map.PixVal(map.PixIndexSph(thetas,phis(i)))=datas(2*i);
      }
    }
    //cout <<"gen "<<datan<<endl;
    for (int i=0;i< nph;i++){
      map.PixVal(map.PixIndexSph(theta,phin(i)))=datan(2*i);
    }
  }
}

void syn_phasgen(const int nsmax,const int nlmax,const int nmmax,
	      const vector< complex<long double> >& datain,
	      int nph, Vector& dataout, double phi0,
	      vector< complex<long double> >& bw){

  /*=======================================================================
     dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 
     with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1

     as the set of frequencies {m} is larger than nph, 
     we wrap frequencies within {0..nph-1}
     ie  m = k*nph + m' with m' in {0..nph-1}
     then
     noting bw(m') = exp(i*m'*phi0) 
                   * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
        with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
     dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
                = Fourier Transform of bw
        is real

         NB nph is not necessarily a power of 2

=======================================================================*/

  int ksign =  1;
  //long double * data = new long double[nph*2];
  complex<double>* data = new complex<double>[nph];
  long double * work = new long double[nph*2];

  int iw;
  for (iw=0;iw<nph;iw++){bw[iw]=0;}

  long double shit=phi0*nph;
  //try complex<long double> since generally this will not be an int
  //int kshift = (int) pow(-1,shit);
  //cout<<shit<<" "<<phi0*nph/M_PI<<" "<<kshift<<endl;
  //  int kshift = (int) pow(-1,kphi0);//  ! either 1 or -1
  //     all frequencies [-m,m] are wrapped in [0,nph-1]
  for (int i=1;i<=2*nmmax+1;i++){
    int m=i-nmmax-1; //in -nmmax, nmmax
    iw=((m % nph) +nph) % nph; //between 0 and nph = m'
    int k=(m-iw)/nph; //number of 'turns'
    //  cout <<(complex<long double>)pow(kshift,k)<<" ";
    //    bw[iw]+=datain[i-1]* (complex<long double>)(pow(kshift,k));//complex number
    bw[iw]+=datain[i-1]* complex<long double>(cos(shit*k),sin(shit*k));//complex number
  }


  //     applies the shift in position <-> phase factor in Fourier space
  for (int iw=1;iw<=nph;iw++){
    int m=ksign*(iw-1);
    complex<long double> shit(cos(m*phi0),sin(m*phi0));
    data[iw-1]=bw[iw-1]*shit;
    //    data[2*(iw-1)]=(bw[iw-1]*shit).real();
    //data[2*(iw-1)+1]=(bw[iw-1]*shit).imag();
  }

  //FFTVector inv((double*) data,nph,true);
  //FFTVector outv=FFTServer::solve(inv,1);
  //outv.Vreal(dataout);

  FFTServer fft;
  fft.fftf(nph, data);
  dataout.Realloc(2*nph);
  for (int i=0;i<nph;i++) {
    dataout(2*i)=data[i].real();dataout(2*i+1)=data[i].imag();
  }
  delete[] data;
  delete[] work;

  /*
  int dum0=1;
  int dum1=1;
    
  fft_gpd_(data,nph,dum0,ksign,dum1,work); // complex to complex 

  for (int iw=0;iw<nph;iw++){dataout(iw) = data[2*iw];}*/
}

void create_almgen(const int nsmax,const int nlmax,const int nmmax,
		vector< vector< complex<double> > >& alm,
		vector<float>& cls_tt,int& iseed,const float fwhm,
		vector<float> lread){


  /*=======================================================================
     creates the a_lm from the power spectrum, 
     assuming they are gaussian  and complex 
     with a variance given by C(l)

     the input file should contain : l and C(l) with *consecutive* l's
     (missing C(l) are put to 0.)

     because the map is real we have : a_l-m = (-)^m conjug(a_lm)
     so we actually compute them only for m >= 0

 modifie G. Le Meur (13/01/98) :
 on ne lit plus les C(l) sur un fichier. Le tableau est entre en argument
 (cls_tt). Ce tableau doit contenir les valeur de C(l) par ordre 
 SEQUENTIEL de l (de l=0 a l=nlmax)
=======================================================================*/
/*      CHARACTER*128 filename

      integer unit,n_l,j,il,im
      real    fs,quadrupole,xxxtmp,correct
      real    fs_tt
      real    zeta1_r, zeta1_i
      LOGICAL ok
      character*20 string_quad

      real*4    prefact(0:2)
      logical   prefact_bool

      integer lneffct
      real    gasdev2
      external lneffct, gasdev2

c-----------------------------------------------------------------------*/
  alm.resize(nlmax+1);
  for (int i=0; i< (signed) alm.size();i++){
    alm[i].resize(nmmax+1);
  }
  //cout << "cing createalm"<<endl;
  iseed = -abs(iseed);
  if (iseed == 0){iseed=-1;}
  int idum = iseed;

  float sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;

  for (int i=0;i<nlmax+1;i++){lread[i]=0;}

  for (int i=0;i <= nlmax;i++){lread[i]  = i;}

  int n_l = nlmax+1;

  cout<<lread[0]<<" <= l <= "<<lread[n_l-1]<<endl;

  //    --- smoothes the initial power spectrum ---

  for (int i=0;i<n_l;i++){
    int l= (int) lread[i];
    float gauss=exp(-l*(l+1.)*sig_smooth*sig_smooth);
    cls_tt[i]*=gauss;
  }
  

  //     --- generates randomly the alm according to their power spectrum ---
  float hsqrt2 = 1.0 / sqrt(2.0);

  for (int i=0;i<n_l;i++){
    int l=(int) lread[i];
    float rms_tt=sqrt(cls_tt[i]);
    //        ------ m = 0 ------
    complex<float> zeta1(gasdev2_(idum));
    //cout << "c2 "<<idum << " "<<zeta1 <<endl;
    alm[l][0]   = zeta1 * rms_tt;

    //------ m > 0 ------
    for (int m=1;m<=l;m++){
      complex<float> shit1(hsqrt2);
      //    cout << "c "<<idum << endl;
      complex<float> shit2(gasdev2_(idum),gasdev2_(idum));
      zeta1=shit1*shit2;
      alm[l][m]=rms_tt*zeta1;
    }
  }
}