Index: /trunk/SophyaLib/Samba/alm.h
===================================================================
--- /trunk/SophyaLib/Samba/alm.h	(revision 2957)
+++ /trunk/SophyaLib/Samba/alm.h	(revision 2958)
@@ -21,11 +21,11 @@
 Alm() : TriangularMatrix<complex<T> >()  {;};
     /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
-Alm(int nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
+Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
 Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {;}
 
 Alm(const TVector<T>& clin, const r_8 fwhm) ;
 /*! resize with a new lmax */
-inline void ReSizeToLmax(int_4 nlmax) {this->ReSizeRow(nlmax+1);}
-inline int_4 Lmax() const {return this->rowNumber()-1;}
+inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);}
+inline sa_size_t Lmax() const {return this->rowNumber()-1;}
 TVector<T> powerSpectrum() const;
 
@@ -48,19 +48,19 @@
 Bm(): nmmax_(0) {;};
     /* instanciate from the maximum absolute value of the index */
-Bm(int mmax) : nmmax_((uint_4)mmax)   { bm_.ReSize( (uint_4)(2*mmax+1) );}
+Bm(sa_size_t mmax) : nmmax_(mmax)   { bm_.ReSize( (2*mmax+1) );}
  Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;}
  /*! resize with a new value of mmax */
-inline void ReSizeToMmax(int_4 mmax) 
+inline void ReSizeToMmax(sa_size_t mmax) 
   {
-    nmmax_=(uint_4)mmax;
+    nmmax_= mmax;
     bm_.ReSize(2* nmmax_+1);
   }
-//inline int_4 Size() const {return 2*nmmax_+1;};
-inline T& operator()(int m) {return bm_( adr_i(m));};
-inline T const& operator()(int m) const {return *(bm_.Begin()+ adr_i(m));};
+//inline sa_size_t Size() const {return 2*nmmax_+1;};
+inline T& operator()(sa_size_t m) {return bm_( adr_i(m));};
+inline T const& operator()(sa_size_t m) const {return *(bm_.Begin()+ adr_i(m));};
 /*! return the current value of the maximum absolute value of the index */
-inline int_4 Mmax() const 
+inline sa_size_t Mmax() const 
    {
-     return (int_4)nmmax_;
+     return nmmax_;
    }
 inline Bm<T>& operator = (const Bm<T>& b)
@@ -76,12 +76,12 @@
   private:
 /*! return the address in the array representing the vector */
-inline uint_4 adr_i(int i) const 
+inline sa_size_t adr_i(sa_size_t i) const 
 {
-  return(uint_4)(i+nmmax_);
+  return (i+nmmax_);
 }
 
 
 
- uint_4 nmmax_;
+ sa_size_t nmmax_;
  NDataBlock<T> bm_;
 
Index: /trunk/SophyaLib/Samba/lambdaBuilder.cc
===================================================================
--- /trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 2957)
+++ /trunk/SophyaLib/Samba/lambdaBuilder.cc	(revision 2958)
@@ -82,16 +82,6 @@
 void LambdaLMBuilder::array_init(int lmax, int mmax)
    {
-     if (a_recurrence_.Size() == 0)
-       {
-	 //	 a_recurrence_ = new TriangularMatrix<r_8>;
-	 updateArrayRecurrence(lmax);
-       }
-     else
-       if ( lmax > (a_recurrence_).rowNumber()-1     )
-	 {
-	   cout << " WARNING : The classes LambdaXXBuilder will be more efficient if instanciated with parameter lmax = maximum value of l index which will be needed in the whole application (arrays not recomputed) " << endl;
-	   cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (a_recurrence_).rowNumber() << endl;
-	 updateArrayRecurrence(lmax);
-	 }
+     updateArrayRecurrence(lmax);
+
      lmax_=lmax;
      mmax_=mmax;
@@ -129,4 +119,138 @@
    }
 
+/*!
+  \brief : Specialized/optimized static function for fast spherical harmonic transform.
+  Computes bm(m) = Sum_l>=m [ lambda(l,m) * alm(l,m) ] 
+  See SphericalTransformServer<T>::GenerateFromAlm(map, pixsize, alm)
+*/
+void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
+				     const Alm<r_8>& alm, Bm< complex<r_8> >& bm)
+{
+  updateArrayRecurrence(lmax);
+  r_8 cth = cos(theta);
+  r_8 sth = sin(theta);
+  
+  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
+  
+  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
+  r_8 lam_0, lam_1, lam_2;
+  register r_8 a_rec, b_rec;
+
+  const complex<r_8>* almp = alm.columnData(0); 
+  r_8* arecp = a_recurrence_.columnData(0); 
+
+  int_4 l, m;
+
+  for (m=0; m<=mmax;m++)   {
+    //MOV^     const complex<r_8>* almp = alm.columnData(m); 
+    complex<r_8>* bmp = &(bm(m));
+    
+    *bmp = (lam_mm / bignorm2)*(*almp);  almp++;
+
+    lam_0=0.;
+    lam_1=1. /bignorm2 ;
+    
+    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
+    b_rec = 0.;
+    for (l=m+1; l<=lmax; l++) {
+      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
+      
+      //DEL	   lambda_(l,m) = lam_2*lam_mm;
+      *bmp += lam_2*lam_mm*(*almp);   almp++;
+      b_rec=1./a_rec;
+      //	   a_rec= LWK->a_recurr(l,m);
+      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
+      lam_0 = lam_1;
+      lam_1 = lam_2;
+    }
+    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
+  }
+}
+
+/*!
+  \brief : Specialized/optimized static function for fast spherical harmonic transform.
+*/
+void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax,  
+				     const Alm<r_4>& alm, Bm< complex<r_4> >& bm)
+{
+  Alm<r_8> alm8(alm.Lmax());
+  for(sa_size_t k=0; k<alm8.Size(); k++) 
+    alm8(k) = complex<r_8>((r_8)alm(k).real() , (r_8)alm(k).imag());
+  Bm< complex<r_8> > bm8(bm.Mmax());
+  ComputeBmFrAlm(theta, lmax, mmax, alm8, bm8);
+  for(sa_size_t kk=-bm.Mmax(); kk<=bm.Mmax(); kk++) 
+    bm(kk)= complex<r_4>((r_4)bm8(kk).real() , (r_4)bm8(kk).imag());
+  return;
+}
+
+
+/*!
+  \brief : Specialized/optimized static function for fast spherical harmonic transform.
+  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ] 
+  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm) 
+*/
+void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
+					TVector< complex<r_8> >& phase, Alm<r_8> & alm)
+{
+  updateArrayRecurrence(lmax);
+  r_8 cth = cos(theta);
+  r_8 sth = sin(theta);
+  
+  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
+  
+  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
+  r_8 lam_0, lam_1, lam_2;
+  register r_8 a_rec, b_rec;
+  
+  complex<r_8>* almp = alm.columnData(0); 
+  r_8* arecp = a_recurrence_.columnData(0); 
+
+  int_4 l, m;
+  for (m=0; m<=mmax;m++)  {
+    //MOV^ complex<r_8>* almp = alm.columnData(m); 
+    complex<r_8>  phi = phase(m);
+    
+    *almp += (lam_mm / bignorm2)*phi;    almp++;
+
+    lam_0=0.;
+    lam_1=1. /bignorm2 ;
+    
+    //MOV^ r_8* arecp = a_recurrence_.columnData(m); 
+    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
+    b_rec = 0.;
+    for (l=m+1; l<=lmax; l++)   {
+      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
+      
+      //DEL	   lambda_(l,m) = lam_2*lam_mm;
+      *almp += (lam_2*lam_mm) * phi;
+      almp++;
+      
+      b_rec=1./a_rec;
+      //	   a_rec= LWK->a_recurr(l,m);
+      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
+      lam_0 = lam_1;
+      lam_1 = lam_2;
+    }
+    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
+  }
+}
+
+/*!
+  \brief : Specialized/optimized static function for fast spherical harmonic transform.
+  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ] 
+  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm) 
+*/
+void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
+					TVector< complex<r_4> >& phase, Alm<r_4> & alm) 
+{
+  TVector< complex<r_8> > phase8;
+  phase8 = phase;
+  Alm<r_8> alm8(alm.Lmax());
+  ComputeAlmFrPhase(theta, lmax, mmax, phase8, alm8);
+  for(sa_size_t k=0; k<alm.Size(); k++) 
+    alm(k) = complex<r_4>((r_4)alm8(k).real() , (r_4)alm8(k).imag());
+  return;
+}
+
 
 /*! \fn void  SOPHYA::LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
@@ -136,4 +260,15 @@
 void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
    {
+     if ( (a_recurrence_.Size() > 0) && (lmax < a_recurrence_.rowNumber()) ) 
+       return;  // Pas besoin de recalculer le tableau de recurrence
+     if (a_recurrence_.Size() > 0)  {
+       cout << " WARNING : The classes LambdaXXBuilder will be more efficient "
+	    << "if instanciated with parameter lmax = maximum value of l index "
+	    << "which will be needed in the whole application (arrays not "
+	    << "recomputed) " << endl;
+       cout << "   lmax= " << lmax << " previous instanciation with lmax= " 
+	    <<  a_recurrence_.rowNumber()-1 << endl;
+	 }
+
      a_recurrence_.ReSizeRow(lmax+1);
      for (int m=0; m<=lmax;m++)
Index: /trunk/SophyaLib/Samba/lambdaBuilder.h
===================================================================
--- /trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 2957)
+++ /trunk/SophyaLib/Samba/lambdaBuilder.h	(revision 2958)
@@ -46,5 +46,19 @@
 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() {};
+
+
+//--- Optimisation Reza mai 2006 Serie de fonctions statiques 
+// specialisees/optimisees pour SphericalTransform
+static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
+			   const Alm<r_8>& alm, Bm< complex<r_8> >& bm);
+static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
+			   const Alm<r_4>& alm, Bm< complex<r_4> >& bm);
+static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
+			      TVector< complex<r_8> >& phase, Alm<r_8> & alm);
+static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax, 
+			      TVector< complex<r_4> >& phase, Alm<r_4> & alm);
+//---- Fin foctions specialisees
 
  /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
@@ -65,5 +79,5 @@
 
  private:
- void updateArrayRecurrence(int_4 lmax);
+ static void updateArrayRecurrence(int_4 lmax);
 
  void array_init(int lmax, int mmax);
Index: /trunk/SophyaLib/Samba/sphericaltransformserver.cc
===================================================================
--- /trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 2957)
+++ /trunk/SophyaLib/Samba/sphericaltransformserver.cc	(revision 2958)
@@ -232,8 +232,10 @@
       //	      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
-      int m;
-      for (m = 0; m <= nmmax; m++)
+      for (int_4 m = 0; m <= nmmax; m++)
 	{
 	  b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
@@ -243,7 +245,11 @@
 	    }
 	}
+      ------- Fin version PRE-Mai2006 */
+      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 (m=1;m<=nmmax;m++)
+      for (int_4 m=1;m<=nmmax;m++)
 	{
 	  b_m_theta(-m) = conj(b_m_theta(m));
@@ -555,4 +561,8 @@
 //	for each m and l
 //	-----------------------------------------------------------------------
+
+      // ===> Optimisation Reza, Mai 2006 
+      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
+        qui calcule la somme au vol 
       //        PrtTim("avant instanciation LM ");
       LambdaLMBuilder lb(theta,nlmax,nmmax);
@@ -576,4 +586,9 @@
 	    }
 	}
+      ------- Fin version PRE-Mai2006 */
+      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
+      phase *= complex<T>((T)domega, 0.);
+      LambdaLMBuilder::ComputeAlmFrPhase(theta,nlmax,nmmax, phase, alm);
+      //Fin Optimisation Reza, Mai 2006 <==== 
       
 
