source: Sophya/trunk/SophyaLib/Samba/ana2fast.cc@ 502

#include <math.h>
#include <vector>
#include <fftserver.h>
#include <complex>
#include "ana2fast.h"
#include "lambuilder.h"

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

void map2a2lm(int nsmax,int nlmax,int nmmax,const vector<float>& mapq,
              const vector<float>& mapu,
             vector< vector< complex<double> > >& a2lme,
              vector< vector< complex<double> > >& a2lmb,
              double cos_theta_cut){
  
  //       REAL*4 powspec(0:nlmax)
    
  //       integer npmiss,npmt,id_miss(10000)
  
  //create the maps for which there are nice basis functions
  
  vector< complex<float> > mapp(mapq.size());
  vector< complex<float> > mapm(mapq.size());
  for (int i=0;i< (signed) mapq.size();i++){
    mapp[i]=complex<float>(mapq[i],mapu[i]);
    mapm[i]=complex<float>(mapq[i],-mapu[i]);
    //cout <<"the maps"<< mapp[i]<<" "<<mapm[i]<<endl;
  }

  vector< vector< complex<double> > > a2lmp;
  vector< vector< complex<double> > > a2lmm;
  a2lmp.resize(nlmax+1);
  for (int i=0; i< (signed) a2lmp.size();i++){
    a2lmp[i].resize(nmmax+1);
    for (int j=0; j< (signed) a2lmp[i].size();j++)a2lmp[i][j]=0;
  }
  a2lmm.resize(nlmax+1);
  for (int i=0; i< (signed) a2lmm.size();i++){
    a2lmm[i].resize(nmmax+1);
    for (int j=0; j< (signed) a2lmm[i].size();j++)a2lmm[i][j]=0;
  }

  /*-----------------------------------------------------------------------
    computes the integral in phi : phas_m(theta)
    for each parallele from north to south pole
    -----------------------------------------------------------------------*/
  
  int istart_north = 0;
  int  istart_south = 12*nsmax*nsmax;
    
  double dth1 = 1. / (3.*nsmax*nsmax);
  double  dth2 = 2. / (3.*nsmax);
  double  dst1 = 1. / (sqrt(6.) * nsmax);

  vector< complex<double> > phas_np(nmmax+1), phas_sp(nmmax+1),
    phas_nm(nmmax+1),phas_sm(nmmax+1);

  for (int ith = 1; ith <= 2*nsmax;ith++){
    int nph, kphi0;
    double cth, sth, sth2;
    //assign doesn't seem to exist in our compiler
    //phas_n.assign(nmmax+1,(complex<float>) 0);
    //phas_s.assign(nmmax+1,(complex<float>) 0);
    for (int i=0;i< nmmax+1;i++){
      phas_np[i]=0; phas_sp[i]=0;phas_nm[i]=0;phas_sm[i]=0;
    }

    if (ith <= nsmax-1){      /* north polar cap */
      nph = 4*ith;
      kphi0 = 1; 
      cth = 1.  - dth1*ith*ith; /* cos(theta) */
      sth = sin( 2. * asin( ith * dst1 ) ) ;  /* sin(theta) */
      sth2 = sth*sth;
    } else { /* tropical band + equat. */
      nph = 4*nsmax;
      kphi0 = (ith+1-nsmax) % 2;
      cth = (2.*nsmax-ith) * dth2;
      sth = sqrt((1.-cth)*(1.+cth)); /* ! sin(theta)*/
      sth2=(1.-cth)*(1.+cth);
    }

    //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_phas2_2(nsmax,nlmax,nmmax,mapp,mapm,istart_north,nph,phas_np,
                   phas_nm,kphi0);
    }
    istart_north = istart_north + nph;
            
    istart_south = istart_south - nph;
    if (ith < 2*nsmax && keep_it){
      comp_phas2_2(nsmax,nlmax,nmmax,mapp,mapm,istart_south,nph,phas_sp,
                   phas_sm,kphi0);
    }
    /*-----------------------------------------------------------------------
      computes the a_lm by integrating over theta
      lambda_lm(theta) * phas_m(theta)
      for each m and l
      -----------------------------------------------------------------------*/
    Lambda2Builder l2b(acos(cth),nlmax,nmmax);
    //cout << "fft:"<<phas_np[0]<<" "<<phas_sp[0]<<" "<<phas_nm[0]<<" "<<phas_sm[0]<<endl;
    for (int m = 0; m <= nmmax; m++){
      //      cout << phas_np[m]<<" "<<phas_sp[m]<<" "<<phas_nm[m]<<" "<<phas_sm[m]<<endl;
      a2lmp[m][m]+=l2b.lam2lmp(m,m)*phas_np[m]+l2b.lam2lmp(m,m,-1)*phas_sp[m];
      a2lmm[m][m]+=l2b.lam2lmm(m,m)*phas_nm[m]+l2b.lam2lmm(m,m,-1)*phas_sm[m];
      for (int l = m+1; l<= nlmax; l++){
        a2lmp[l][m]+=
          l2b.lam2lmp(l,m)*phas_np[m]+l2b.lam2lmp(l,m,-1)*phas_sp[m];
        a2lmm[l][m]+=
          l2b.lam2lmm(l,m)*phas_nm[m]+l2b.lam2lmm(l,m,-1)*phas_sm[m];
      }
    }
  }
  complex<double> im(0,1);
  a2lme.resize(nlmax+1);
  for (int i=0; i< (signed) a2lme.size();i++){
    a2lme[i].resize(nmmax+1);
  }
  a2lmb.resize(nlmax+1);
  for (int i=0; i< (signed) a2lmb.size();i++){
    a2lmb[i].resize(nmmax+1);
  }
  float domega=(4.*M_PI)/(12.*nsmax*nsmax);
  for (int m = 0; m <= nmmax; m++){
    for (int l = m+1; l<= nlmax; l++){
      a2lme[l][m]=-(a2lmp[l][m]+a2lmm[l][m])/2.*(double)domega;
      a2lmb[l][m]=im*(a2lmp[l][m]-a2lmm[l][m])/2.*(double)domega;
    }
  }
  //for (int l = 2; l<= nlmax; l++){
  //cout << "calc almp,m"<<a2lmp[l][0]<<" "<<a2lmm[l][0]<<endl;}
}

void comp_phas2_2(int nsmax,int nlmax,int nmmax,
                  const vector< complex<float> >& datain,
                  const vector< complex<float> >& datain2,
                int start,int nph,vector< complex<double> >& dataout,
                 vector< complex<double> >& dataout2, int kphi0){
  /*=======================================================================
    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
    =======================================================================*/

  int ksign = -1;
  double phi0 = kphi0*M_PI/nph;

  complex<double>* data= new complex<double>[4*nsmax];
  complex<double>* data2= new complex<double>[4*nsmax];
  for (int i = 0; i< nph;i++){
    data[i] = datain[i+start];
    data2[i] = datain2[i+start];
  }
  for (int i = nph; i< 4*nsmax;i++){
    data[i] = 0;
    data2[i] = 0;
  }

  FFTServer fft;
  fft.fftb(nph,data);
  fft.fftb(nph,data2);

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