Index: /trunk/SophyaLib/Samba/sphericaltransformserver.cc
===================================================================
--- /trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 2990)
+++ /trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 2991)
@@ -180,73 +180,57 @@
   Bm<complex<T> > b_m_theta(nmmax);
 
-
-
   // pour chaque tranche en theta
-  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)
-    {
-      int_4 nph;
-      r_8 phi0;
-      r_8 theta;
-      TVector<int_4> pixNumber; 
-      TVector<T> datan;
-      
-      map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
-      nph = pixNumber.NElts();
-      if (nph < 2) continue;  // On laisse tomber les tranches avec un point
-      //       -----------------------------------------------------
-      //	      for each theta, and each m, computes
-      //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
-      //	      ------------------------------------------------------
-      // ===> Optimisation Reza, Mai 2006 
-      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
-        qui calcule la somme au vol 
+  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)  {
+    int_4 nph;
+    r_8 phi0;
+    r_8 theta;
+    TVector<int_4> pixNumber; 
+    TVector<T> datan;
+    
+    map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
+    nph = pixNumber.NElts();
+    if (nph < 2) continue;  // On laisse tomber les tranches avec un point
+    //       -----------------------------------------------------
+    //	      for each theta, and each m, computes
+    //	      b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m) 
+    //	      ------------------------------------------------------
+    // ===> Optimisation Reza, Mai 2006 
+    /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
+      qui calcule la somme au vol 
       LambdaLMBuilder lb(theta,nlmax,nmmax);
       //  somme sur m de 0 a l'infini
-      for (int_4 m = 0; m <= nmmax; m++)
-	{
-	  b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
-	  for (int l = m+1; l<= nlmax; l++)
-	    {
-	       b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
-	    }
-	}
+      for (int_4 m = 0; m <= nmmax; m++) {
+      b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
+      for (int l = m+1; l<= nlmax; l++)
+        b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
+      }
       ------- Fin version PRE-Mai2006 */
-      LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
-      //Fin Optimisation Reza, Mai 2006 <==== 
+    LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
+    //Fin Optimisation Reza, Mai 2006 <==== 
 
     //        obtains the negative m of b(m,theta) (= complex conjugate)
-
-      for (int_4 m=1;m<=nmmax;m++)
-	{
-	  b_m_theta(-m) = conj(b_m_theta(m));
-	}
-      // ---------------------------------------------------------------
-      //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
-      // ---------------------------------------------------------------*/
-
-      /* ----- Reza, Juin 2006 : 
-         En verifiant la difference entre deux cartes 
-         cl -> map -> alm -> map2 et mapdiff = map-map2 
-         je me suis apercu qu'il y avait des differences importantes - dans les
-         deux zones 'polar cap' de HEALPix - apres avoir pas mal chercher, il semble
-         que la routine RfourierSynthesisFromB en soit responsable 
-         Je fais donc tout passer dans fourierSynthesisFromB 
-      if ( (sphere_type == "RING") || (sphere_type == "NESTED") ) {
-	  TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
-	  for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i); 
-	}
-      else 
-      */
-	// pour des pixelisations quelconques (autres que HEALPix
-	//  nph n'est pas toujours pair
-	// ca fait des problemes pour les transformees de Fourier
-	// car le server de TF ajuste la longueur du vecteur reel 
-	// en sortie de TF, bref, la securite veut qu'on prenne une 
-	// TF complexe
-	{
-	  TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0);
-	  for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i).real();
-	}
-    }
+    for (int_4 m=1;m<=nmmax;m++)
+      b_m_theta(-m) = conj(b_m_theta(m));
+    // ---------------------------------------------------------------
+    //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
+    // ---------------------------------------------------------------*/
+
+    /* ----- Reza, Juin 2006 : 
+       En verifiant la difference entre deux cartes 
+       cl -> map -> alm -> map2 et mapdiff = map-map2 
+       je me suis apercu qu'il y avait des differences importantes - dans les
+       deux zones 'polar cap' de HEALPix - qui utilisait RfourierSynthesisFromB
+       TF complex -> reel . Le probleme venant de l'ambiguite de taille, lie 
+       a la partie imaginaire de la composante a f_nyquist , j'ai corrige et 
+       tout mis en TF complexe -> reel 
+    */
+    TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
+    // Si on peut acceder directement les pixels d'un tranche, on le fait
+    T* pix = map.GetThetaSliceDataPtr(ith);
+    if (pix != NULL) 
+      for (int_4 i=0;i< nph;i++) pix[i] = Temp(i);
+    else 
+      for (int_4 i=0;i< nph;i++) map(pixNumber(i))=Temp(i); 
+  }
 }
 
@@ -380,38 +364,30 @@
   //**********************************************************************
   TVector< complex<T> > bw(nph);
-  TVector< complex<T> > dataout(nph);
   TVector< complex<T> > data(nph/2+1);
-
 
   for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
   int m;
-  for (m=-b_m.Mmax();m<=-1;m++)
-    {
-      int maux=m;
-      while (maux<0) maux+=nph;
-      int iw=maux%nph;
-      double aux=(m-iw)*phi0;
-      bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
-    }
-  for (m=0;m<=b_m.Mmax();m++)
-    {
-      //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
-      int iw=m%nph;
-      double aux=(m-iw)*phi0;
-      bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
-    }
+  for (m=-b_m.Mmax();m<=-1;m++)  {
+    int maux=m;
+    while (maux<0) maux+=nph;
+    int iw=maux%nph;
+    double aux=(m-iw)*phi0;
+    bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) )  ;
+  }
+  for (m=0;m<=b_m.Mmax();m++) {
+    //      int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
+    int iw=m%nph;
+    double aux=(m-iw)*phi0;
+    bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
+  }
 
   //     applies the shift in position <-> phase factor in Fourier space
-  // cout << " TF : nph= " << nph << " vec. entree " << data.Size() << endl;
-  for (int mprime=0; mprime <= nph/2; mprime++)
-    {
-      complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
-      data(mprime)=bw(mprime)*
-	               (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
-    }
-
-  TVector<T> sortie;
-  fftIntfPtr_-> FFTBackward(data, sortie);
-  
+  for (int mprime=0; mprime <= nph/2; mprime++)  
+    data(mprime)=bw(mprime)*complex<T>((T)cos(mprime*phi0),(T)sin(mprime*phi0));    
+  TVector<T> sortie(nph);
+// On met la partie imaginaire du dernier element du data a zero pour nph pair
+  if (nph%2 == 0) data(nph/2) = complex<T>(data(nph/2).real(), (T)0.);
+// et on impose l'utilisation de la taille en sortie pour FFTBack (..., ..., true)
+  fftIntfPtr_-> FFTBackward(data, sortie, true); 
   return sortie;
 }
