Index: trunk/SophyaLib/Samba/ana2fast.cc
===================================================================
--- trunk/SophyaLib/Samba/ana2fast.cc	(revision 730)
+++ 	(revision )
@@ -1,193 +1,0 @@
-#include <math.h>
-#include <vector>
-#include <fftserver.h>
-#include <complex>
-#include "ana2fast.h"
-#include "lambuilder.h"
-#ifdef __MWERKS__
-#include "unixmac.h"
-#endif
-
-
-/*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++){
-      a2lme[m][m]=-(a2lmp[m][m]+a2lmm[m][m])/2.*static_cast<double>(domega);
-      a2lmb[m][m]=im*(a2lmp[m][m]-a2lmm[m][m])/2.*static_cast<double>(domega);
-    for (int l = m+1; l<= nlmax; l++){
-      a2lme[l][m]=-(a2lmp[l][m]+a2lmm[l][m])/2.*static_cast<double>(domega);
-      a2lmb[l][m]=im*(a2lmp[l][m]-a2lmm[l][m])/2.*static_cast<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;
-}
