Changeset 2958 in Sophya for trunk/SophyaLib/Samba/lambdaBuilder.cc

Timestamp:

Jun 1, 2006, 1:34:50 PM (19 years ago)

Author:

ansari

Message:

1/ passage en sa_size_t (au lieu de int) dans Alm<T> et Bm<T>
2/ Ajout des methodes optimisees (statiques) pour calcul transforme Ylm
ds LambdaLMBuilder et utilisation ds SphericalTransformServer

Reza , 1/06/2006

File:

• trunk/SophyaLib/Samba/lambdaBuilder.cc (modified) (3 diffs)

trunk/SophyaLib/Samba/lambdaBuilder.cc

@@ -82 +82 @@
 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 +119 @@
    }
 
+/*!
+  \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 +260 @@
 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++)