(JEC) 7/10/15 add Spherical Harmonic synthesis in Laguerre2Bessel::Synthesis

src/lagsht_testsuite.cc

@@ -349,7 +349,7 @@ void TestLaguerre2Bessel() {
 
   cout << "tstack 2 Start" <<endl;
   tstack_push("processing");
-  int Pmax = 32*Nmax; //To see 
+  int Pmax = 64*Nmax; //To see 
   cout << "tstack 2a Start" <<endl;
   tstack_push("Init Laguerre 2 Bessel");
   Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax);
@@ -368,7 +368,7 @@ void TestLaguerre2Bessel() {
   cout << "tstack 2c Start" <<endl;
   tstack_push("Compute Alm(rk) from FB coeffs.");
   vector< complex<r_8> > FBalmk;
-  vector<r_8> fFBijk; //Not yet computed  JEC 22/9/15
+  vector<r_8> fFBijk;
   lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
   cout << "tstack 2c End" <<endl;
   tstack_pop("Compute Alm(rk) from FB coeffs.");
@@ -382,60 +382,52 @@ void TestLaguerre2Bessel() {
   cout << "tstack 2d End" <<endl;
   tstack_pop("Compute Alm(rk) from FB coeffs.");
   
-   cout << "tstack 2 End" <<endl;
-   tstack_pop("processing");
-
-//   //FLaguerre Coeff.
-//   for(int n=0;n<Nmax;n++){
-//     for(int l=0;l<Lmax;l++){
-//       for (int m=0;m<=l;m++) {
-// 	int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
-// 	cout << "FL("
-// 	     <<l<<","
-// 	     <<m<<","
-// 	     <<n<<"): " << setprecision(12) << flmnOrig[id] << endl;
-//       }
-//     }
-//   }
-
-// //   //FBessel Coeff.
-// //   for(int p=0;p<Pmax;p++){
-// //     for(int l=0;l<Lmax;l++){
-// //       for (int m=0;m<=l;m++) {
-// // 	int id= p*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
-// // 	cout << "FB("
-// // 	     <<l<<","
-// // 	     <<m<<","
-// // 	     <<p<<"): " << setprecision(12) << FBlmp[id] << endl;
-// //       }
-// //     }
-// //   }
-
-  cout << "Dump FL or FB  reconstructed Alm(r_k)" <<endl;
-
-   std::ofstream ofs;
-   ofs.open ("Almi.txt", std::ofstream::out);
-
-  for(int n=0;n<Nmax;n++){
-    for(int l=0;l<Lmax;l++){
-      for (int m=0;m<=l;m++) {
-	int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
-	cout << "Almk ("
-	     <<l<<","
-	     <<m<<","
-	     <<n<< "): FL : " << setprecision(12) << FLalmk[id] 
-	     << " FB : " << FBalmk[id]  << endl;
-	ofs <<l<<" "
-	     <<m<<" "
-	    <<n<< " " << setprecision(12) << FLalmk[id].real() << " " <<   FLalmk[id].imag()
-	    << "  " << FBalmk[id].real() << " " << FBalmk[id].imag()  << endl;
+  cout << "tstack 2 End" <<endl;
+  tstack_pop("processing");
+  
+  
+  {//dump 1
+    cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
 
+    std::ofstream ofs;
+    ofs.open ("Almi.txt", std::ofstream::out);
+    
+    for(int n=0;n<Nmax;n++){
+      for(int l=0;l<Lmax;l++){
+	for (int m=0;m<=l;m++) {
+	  int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
+	  cout << "Almk ("
+	       <<l<<","
+	       <<m<<","
+	       <<n<< "): FL : " << setprecision(12) << FLalmk[id] 
+	       << " FB : " << FBalmk[id]  << endl;
+	  ofs <<l<<" "
+	      <<m<<" "
+	      <<n<< " " << setprecision(12) << FLalmk[id].real() << " " <<   FLalmk[id].imag()
+	      << "  " << FBalmk[id].real() << " " << FBalmk[id].imag()  << endl;
+	  
+	}
+      }
+    }
+
+    ofs.close();
+  }//dump 1
+
+  {//dump 2
+    std::ofstream ofs;
+    ofs.open ("fijk.txt", std::ofstream::out);
+
+    for (int ish=0;ish<Nshell; ish++){
+      for (int ip=0; ip<Npix; ip++) {
+	int id = ish*Npix+ip;
+	cout << "Shell("<<ish<<") Pix["<<ip<<"] FL fijk = " << fFLijk[id] << " FB fijk = " << fFBijk[id] << endl;
+	ofs << ish << " " << ip << " " << setprecision(12) <<  fFLijk[id] << " " << fFBijk[id] << endl;
       }
     }
-  }
 
-  ofs.close();
 
+    ofs.close();
+  }//dump 2
 
   
 

src/laguerre2bessel.cc

@@ -145,7 +145,7 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
 
 void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp, 
 				vector< complex<r_8> >& FBalmk,
-				vector<r_8>& fijk //not computed for the moment
+				vector<r_8>& fijk
 				) {
 
   int NtotCoeffFB = Nalm_*Pmax_; 
@@ -153,8 +153,6 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 
   cout << "Dump  Laguerre2Bessel::Synthesis START "<<endl;
     
-  //  int Ntot3DPix = Npix_*Nmax_;
-  //  fijk.resize(Ntot3DPix); //Not used for the moment
 
   int NtotCoeffFL = Nalm_*Nmax_; 
   FBalmk.resize(NtotCoeffFL); 
@@ -214,6 +212,21 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
     }//loop on l
   }//loop on k
 
+
+  int Ntot3DPix = Npix_*Nmax_;
+  fijk.resize(Ntot3DPix); //Not used for the moment
+
+  int off7n, off7k;
+  //perform the synthesis per shell
+  for(int n =0; n<Nmax_; n++){
+    off7n = n*Nalm_;
+    off7k = n*Npix_; 
+
+    sphere_->alm2map(reinterpret_cast<const r_8*>(&(FBalmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
+
+  }//loop on shell
+  
+
 }//Synthesis