Index: unk/SophyaLib/Samba/syngen.cc
===================================================================
--- /trunk/SophyaLib/Samba/syngen.cc	(revision 744)
+++ 	(revision )
@@ -1,271 +1,0 @@
-#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;
-    }
-  }
-}
Index: unk/SophyaLib/Samba/syngen.h
===================================================================
--- /trunk/SophyaLib/Samba/syngen.h	(revision 744)
+++ 	(revision )
@@ -1,27 +1,0 @@
-#ifndef SYNGEN_H_SEEN
-#define SYNGEN_H_SEEN
-
-#include <complex>
-#include <vector>
-#include "tvector.h"
-#include "sphericalmap.h"
-
-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);
-
-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);
-
-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);
-#endif
