Changeset 1683 in Sophya for trunk/SophyaLib/Samba

Timestamp:

Oct 11, 2001, 5:22:03 PM (24 years ago)

Author:

lemeur

Message:

methode iterative pour analyse harmonique

Location:

trunk/SophyaLib/Samba

Files:

• alm.cc (modified) (1 diff)
• alm.h (modified) (1 diff)
• lambdaBuilder.cc (modified) (1 diff)
• lambdaBuilder.h (modified) (1 diff)
• sphericaltransformserver.cc (modified) (15 diffs)
• sphericaltransformserver.h (modified) (2 diffs)

trunk/SophyaLib/Samba/alm.cc

@@ -78 +78 @@
   return powsp;
 }
-
+ 
 #ifdef __CXX_PRAGMA_TEMPLATES__
 #pragma define_template Alm<r_8>

trunk/SophyaLib/Samba/alm.h

@@ -30 +30 @@
 TVector<T> powerSpectrum() const;
 
+inline  Alm<T>& SetComplexT(complex<T> a)
+   {
+     this->SetT(a);
+     return *this;
+   }
+ inline Alm<T>& operator = (complex<T> a)
+ {
+   return SetComplexT(a);
+ } 
 
 };

trunk/SophyaLib/Samba/lambdaBuilder.cc

@@ -71 +71 @@
       cth_=cos(theta);
       sth_=sin(theta);
+      array_init(lmax, mmax);
+    }
+LambdaLMBuilder::LambdaLMBuilder(r_8 costet, r_8 sintet,int_4 lmax, int_4 mmax)
+    {
+      cth_=costet;
+      sth_=sintet;
       array_init(lmax, mmax);
     }

trunk/SophyaLib/Samba/lambdaBuilder.h

@@ -41 +41 @@
 LambdaLMBuilder() {}
 LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax);
+LambdaLMBuilder(r_8 ct, r_8 st,int_4 lmax, int_4 mmax);
 virtual ~LambdaLMBuilder() {};
 
  /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
 inline double lamlm(int l, int m) const { return lambda_(l,m); }
+
+inline double lamlm(int index) const { return lambda_(index); }
+
+
+  //Return pointer to first element address of the alm's
+//inline r_8* Data()  {return lambda_.Data();}
+
 
  private:

trunk/SophyaLib/Samba/sphericaltransformserver.cc

@@ -7 +7 @@
 #include "nbrandom.h"
 #include "nbmath.h"
+#include "timing.h"
+//#include "spherehealpix.h"
+
 
 /*! \class SOPHYA::SphericalTransformServer
@@ -184 +187 @@
         throw IOExc(" unknown type of sphere: " + (string)sphere_type );
       }
-  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
-  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
-  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
+  //  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
+  //  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
+  //  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
   if (nlmax>3*nsmax-1) 
     {
-      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
+      //     cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
       if (strncmp(sphere_type,"TETAFI",6) == 0)
         {
+          cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
           cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
         }
@@ -221 +225 @@
       for (m = 0; m <= nmmax; m++)
         {
-          //  somme sur l de m a l'infini
           b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
-          //    if (ith==0 && m==0)
-          //    {
-          //      cout << " guy: lmm= " <<  lb.lamlm(m,m) << " alm " << alm(m,m) << "b00= " <<  b_m_theta(m) << endl;
-          //    } 
           for (int l = m+1; l<= nlmax; l++)
             {
                b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
-
-
-               //     if (ith==0 && m==0)
-               //        {
-               //          cout << " guy:l=" << l << " m= " << m << " lmm= " <<  lb.lamlm(l,m) << " alm " << alm(l,m) << "b00= " <<  b_m_theta(m) << endl;
-
-               //        } 
-
             }
         }
@@ -245 +236 @@
       for (m=1;m<=nmmax;m++)
         {
-          //compiler doesn't have conj()
           b_m_theta(-m) = conj(b_m_theta(m));
         }
@@ -251 +241 @@
       //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
       // ---------------------------------------------------------------*/
-      //    TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0); 
-          TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
+      TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
       for (int i=0;i< nph;i++)
         {
-          //      map(pixNumber(i))=Temp(i).real();
-                          map(pixNumber(i))=Temp(i);
+          map(pixNumber(i))=Temp(i);
         }
     }
@@ -421 +409 @@
     }
 
-  //sortie.ReSize(nph);
   TVector<T> sortie;
 
@@ -432 +419 @@
  /*! \fn  Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
 
-\return the Alm coefficients from analysis of a temperature map.
+\return the Alm coefficients from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
 
     \param<nlmax> : maximum value of the l index
 
      \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
+
+\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
   */ 
 template<class T>
- Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
+void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
+{
+  int_4  nmmax = nlmax;
+  //  PrtTim("appel  carteVersAlm");
+  carteVersAlm(map, nlmax, cos_theta_cut, alm);
+  //  PrtTim("retour  carteVersAlm");
+  if (iterationOrder > 0)
+    {
+      TVector<int_4> fact(iterationOrder+2);
+      fact(0) = 1;
+      for (int k=1; k <= iterationOrder+1; k++)
+        {
+          fact(k) = fact(k-1)*k;
+        }
+      Alm<T> alm2(alm);
+      T Tzero = (T)0.;
+      complex<T> complexZero = complex<T>(Tzero, Tzero);
+      alm = complexZero;
+      int signe = 1;
+      int nbIteration = iterationOrder+1;
+      for (int k=1; k <= nbIteration; k++)
+        {
+          T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
+          for (int m = 0; m <= nmmax; m++)
+            {
+              for (int l = m; l<= nlmax; l++)
+                {
+                  alm(l,m) += facMult*alm2(l,m); 
+                }
+            }
+          if (k == nbIteration) break;
+          signe = -signe;
+          for (int k=0; k< map.NbPixels(); k++) map(k) = (T)0.;
+          //        synthetize a map from the estimated alm
+          //      PrtTim("appel  GenerateFromAlm");
+          GenerateFromAlm( map, map.SizeIndex(), alm2);
+          //      PrtTim("retour  GenerateFromAlm");
+          alm2 = complexZero;
+          //        analyse the new map
+          //      PrtTim("appel  carteVersAlm");
+          carteVersAlm(map, nlmax, cos_theta_cut, alm2);
+          //      PrtTim("retour  carteVersAlm");
+        }
+    }
+}
+
+template<class T>
+ void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, Alm<T>& alm) const
 {
   
@@ -450 +486 @@
   int_4  nmmax = nlmax;
   TVector< complex<T> > phase(nmmax+1);
-  Alm<T> alm;
+  
   alm.ReSizeToLmax(nlmax);
   for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
@@ -456 +492 @@
       r_8 phi0;
       r_8 theta;
+      //  PrtTim("debut 1ere tranche ");
       map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
-      for (int i=0;i< nmmax+1;i++)
-        {
-          phase(i)=0; 
-        }
+      phase = complex<T>((T)0.,(T)0.);
       double cth = cos(theta);
       
       //part of the sky out of the symetric cut
       bool keep_it = (fabs(cth) >= cos_theta_cut); 
+
+      //    PrtTim("fin 1ere tranche ");
+  
+      if (keep_it)
+        {
+          //      phase = CFromFourierAnalysis(nmmax,data,phi0);
+          //      PrtTim("avant Fourier ");
+          CFromFourierAnalysis(nmmax,data,phase, phi0);
+          //      PrtTim("apres Fourier ");
+
+        }
       
-      if (keep_it)
-        {
-          // tableau datain a supprimer
-          //  TVector<complex<T> > datain(pixNumber.NElts());
-          // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);
-
-          //  phase = CFromFourierAnalysis(nmmax,datain,phi0);
-          phase = CFromFourierAnalysis(nmmax,data,phi0);
-
+//      ---------------------------------------------------------------------
+//      computes the a_lm by integrating over theta
+//      lambda_lm(theta) * phas_m(theta)
+//      for each m and l
+//      -----------------------------------------------------------------------
+      //        PrtTim("avant instanciation LM ");
+      LambdaLMBuilder lb(theta,nlmax,nmmax);
+      //        PrtTim("apres instanciation LM ");
+      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
+
+      //   PrtTim("avant mise a jour Alm ");
+      complex<T> fi;
+      T facteur;
+      int index;
+      for (int m = 0; m <= nmmax; m++)
+        {
+          fi = phase(m);
+          for (int l = m; l<= nlmax; l++)
+            {
+                 index = alm.indexOfElement(l,m);
+                 //  facteur = (T)(lb.lamlm(l,m) * domega);
+                    facteur = (T)(lb.lamlm(index) * domega);
+                  // alm(l,m) += facteur * fi ; 
+                      alm(index) += facteur * fi ; 
+            }
         }
       
-      /*-----------------------------------------------------------------------
-        computes the a_lm by integrating over theta
-        lambda_lm(theta) * phas_m(theta)
-        for each m and l
-        -----------------------------------------------------------------------*/
-      //      LambdaBuilder lb(theta,nlmax,nmmax);
-      LambdaLMBuilder lb(theta,nlmax,nmmax);
-      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
-      for (int m = 0; m <= nmmax; m++)
-        {
-          alm(m,m) += (T)lb.lamlm(m,m) * phase(m) * (T)domega; //m,m even
-          for (int l = m+1; l<= nlmax; l++)
-            {
-              alm(l,m) += (T)lb.lamlm(l,m) * phase(m)*(T)domega; 
-            }
-        }
-    }
-  return alm;
+
+      
+      //
+      //
+      //       PrtTim("apres mise a jour Alm ");
+    }
 }
   /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
@@ -518 +567 @@
   fftIntfPtr_-> FFTForward(datain, transformedData);
 
-  //dataout.ReSize(nmmax+1);
   TVector< complex<T> > dataout(nmmax+1);
 
   int im_max=min(nph,nmmax+1);
   int i;
-  for (i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
+  dataout = complex<T>((T)0.,(T)0.);
+  //  for (i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
   for (i=0;i<im_max;i++) dataout(i)=transformedData(i);
 
 
-  //  for (int i = 0;i <im_max;i++){
-  //    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
-  //  }
   for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
   for (i = 0;i <dataout.NElts();i++){
@@ -542 +588 @@
 same as previous one, but with a "datain" which is real (not complex) */
 template<class T>
-TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const
+void SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, TVector< complex<T> >& dataout, r_8 phi0) const
 {
   //=======================================================================
@@ -559 +605 @@
       throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
     }
+  if (nph%2 != 0 )
+    {
+      throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
+    }
   TVector<complex<T> > transformedData;
-    // a remodifier
-  //FFTPackServer ffts;
-  //ffts.setNormalize(false);
-  //ffts.FFTForward(datain, transformedData);
-
+
+  // la taille du vecteur complexe retourne est nph/2+1 (si la taille 
+  // du vecteur reel entre est nph)
   fftIntfPtr_-> FFTForward(datain, transformedData);
-    //
-
-  //dataout.ReSize(nmmax+1);
-  TVector< complex<T> > dataout(nmmax+1);
+
+  //  TVector< complex<T> > dataout(nmmax+1);
+  dataout.ReSize(nmmax+1);
 
   // on transfere le resultat de la fft dans dataout.
-  // on s'assure que ca ne depasse pas la taille de dataout
-  int sizeOfTransformToGet = min(transformedData.NElts(),nmmax+1);
-  //  int im_max=min(transformedData.NElts()-1,nmmax);
+
+  int maxFreqAccessiblesParFFT = min(nph/2,nmmax);
   int i;
-  for (i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);
-
-
-  // si dataout n'est pas plein, on complete jusqu'a  nph valeurs (a moins 
+  for (i=0;i<=maxFreqAccessiblesParFFT;i++) dataout(i)=transformedData(i);
+
+
+  // si dataout n'est pas plein, on complete jusqu'a  nph+1 valeurs (a moins 
   // que dataout ne soit plein avant d'atteindre nph)
-  if (sizeOfTransformToGet == (transformedData.NElts()))
-    {
-      for (i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
-        {
-
-          //      dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
+  if (maxFreqAccessiblesParFFT != nmmax )
+    {
+      int maxMfft = min(nph,nmmax);
+      for (i=maxFreqAccessiblesParFFT+1; i<=maxMfft; i++)
+        {
           dataout(i) = conj(dataout(nph-i) );
         }
       // on conplete, si necessaire, par periodicite
-      for (int kk=nph; kk<dataout.NElts(); kk++) 
-        {
-          dataout(kk)=dataout(kk%nph);
-        }
-    }
-  for (i = 0;i <dataout.NElts();i++){
-    dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
-  }
-  return dataout;
+      if ( maxMfft != nmmax )
+        {
+          for (int kk=nph+1; kk <= nmmax; kk++) 
+            {
+              dataout(kk)=dataout(kk%nph);
+            }
+        }
+    }
+  for (i = 0;i <dataout.NElts();i++)
+    {
+      dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
+    }
+  //  return dataout;
 }
 
@@ -701 +750 @@
  */ 
 template<class T>
-void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
+void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& mapq,
+                                              SphericalMap<T>& mapu,
+                                              Alm<T>& alme,
+                                              Alm<T>& almb,
+                                              int_4 nlmax,
+                                              r_8 cos_theta_cut, 
+                                              int iterationOrder) const
+{
+  int_4  nmmax = nlmax;
+  carteVersAlm(mapq, mapu, alme, almb, nlmax, cos_theta_cut);
+  if (iterationOrder > 0)
+    {
+      TVector<int_4> fact(iterationOrder+2);
+      fact(0) = 1;
+      for (int k=1; k <= iterationOrder+1; k++)
+        {
+          fact(k) = fact(k-1)*k;
+        }
+      Alm<T> alme2(alme);
+      Alm<T> almb2(almb);
+      T Tzero = (T)0.;
+      complex<T> complexZero = complex<T>(Tzero, Tzero);
+      alme = complexZero;
+      almb = complexZero;
+      int signe = 1;
+      int nbIteration = iterationOrder+1;
+      for (int k=1; k <= nbIteration; k++)
+        {
+          T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
+          for (int m = 0; m <= nmmax; m++)
+            {
+              for (int l = m; l<= nlmax; l++)
+                {
+                  alme(l,m) += facMult*alme2(l,m); 
+                  almb(l,m) += facMult*almb2(l,m); 
+                }
+            }
+          if (k == nbIteration) break;
+          signe = -signe;
+          for (int k=0; k< mapq.NbPixels(); k++)
+            {
+              mapq(k) = (T)0.;
+              mapu(k) = (T)0.;
+            }
+          //        synthetize a map from the estimated alm
+          GenerateFromAlm(mapq,mapu,mapq.SizeIndex(),alme2,almb2); 
+          alme2 = complexZero;
+          almb2 = complexZero;
+          //        analyse the new map
+          carteVersAlm(mapq, mapu, alme2, almb2, nlmax, cos_theta_cut);
+        }
+    }
+}
+
+template<class T>
+void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& mapq,
                                               const SphericalMap<T>& mapu,
                                               Alm<T>& alme,
@@ -797 +901 @@
   //  for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
 
-  // phaseq = CFromFourierAnalysis(nmmax,datain,phi0); 
-  phaseq = CFromFourierAnalysis(nmmax,dataq,phi0); 
-
-  // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(datau(kk),0.);
-
-  // phaseu=  CFromFourierAnalysis(nlmax,nmmax,datain,phi0); 
-  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 
+  //  phaseq = CFromFourierAnalysis(nmmax,dataq,phi0); 
+  CFromFourierAnalysis(nmmax,dataq,phaseq, phi0); 
+
+  //  phaseu=  CFromFourierAnalysis(nmmax,datau,phi0); 
+  CFromFourierAnalysis(nmmax,datau,phaseu, phi0); 
 
   LambdaWXBuilder lwxb(theta,nlmax,nmmax);
@@ -1141 +1243 @@
 /*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
 
-\return power spectrum from analysis of a temperature map.
+\return power spectrum from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0) 
 
      \param<nlmax> : maximum value of the l index
 
      \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
+
+\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
+
   */ 
 template <class T>
-TVector<T>  SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
-{
-    Alm<T> alm=DecomposeToAlm( sph, nlmax, cos_theta_cut);
+TVector<T>  SphericalTransformServer<T>::DecomposeToCl(SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
+{
+  Alm<T> alm;
+  DecomposeToAlm( sph, alm, nlmax, cos_theta_cut, iterationOrder);
+    //    cout << " SphericalTransformServer : impression alm " << endl;
+    //     alm.Print();
   // power spectrum
      return  alm.powerSpectrum();

trunk/SophyaLib/Samba/sphericaltransformserver.h

@@ -46 +46 @@
 
 
-  Alm<T> DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const;
+void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
 
-  void DecomposeToAlm(const SphericalMap<T>& mapq, const SphericalMap<T>& mapu,
-                     Alm<T>& a2lme, Alm<T>& a2lmb,
-                     int_4 nlmax, r_8 cos_theta_cut) const;
+void DecomposeToAlm(SphericalMap<T>& mapq,
+                    SphericalMap<T>& mapu,
+                    Alm<T>& alme,
+                    Alm<T>& almb,
+                    int_4 nlmax,
+                    r_8 cos_theta_cut, 
+                    int iterationOrder = 0) const;
 
-  TVector<T>  DecomposeToCl(const SphericalMap<T>& sph,  
-                           int_4 nlmax, r_8 cos_theta_cut) const;
+TVector<T>  DecomposeToCl(SphericalMap<T>& sph,  
+                           int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
 
 
  private:
-  TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+
+void carteVersAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, Alm<T>& alm) const;
+
+void carteVersAlm(const SphericalMap<T>& mapq,
+                  const SphericalMap<T>& mapu,
+                  Alm<T>& alme,
+                  Alm<T>& almb,
+                  int_4 nlmax,
+                  r_8 cos_theta_cut) const;
+
+
+TVector< complex<T> >  fourierSynthesisFromB(const Bm<complex<T> >& b_m,  
                                   int_4 nph, r_8 phi0) const;
-  TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
+TVector<T>  RfourierSynthesisFromB(const Bm<complex<T> >& b_m,  
                                   int_4 nph, r_8 phi0) const;
-
-  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
+TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
                                   const TVector<complex<T> > datain,
                                   r_8 phi0) const;
-  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
-                         const TVector<T> datain,  
-                                  r_8 phi0) const;
+  //  TVector< complex<T> > CFromFourierAnalysis(int_4 mmax, 
+  //                     const TVector<T> datain,  
+  //                              r_8 phi0) const;
+  void CFromFourierAnalysis(int_4 mmax, const TVector<T> datain,
+                            TVector< complex<T> >& dataout,
+                            r_8 phi0) const;
 
   void almFromWX(int_4 nlmax, int_4 nmmax, r_8 phi0, 
@@ -90 +107 @@
 
  FFTServerInterface* fftIntfPtr_;
+
 };
 } // Fin du namespace