Changeset 3789 in Sophya for trunk/Cosmo/RadioBeam/fgndsub.cc

Timestamp:

Jun 27, 2010, 4:57:36 PM (15 years ago)

Author:

ansari

Message:

Corrections diverses, Reza 27/06/2010

File:

• trunk/Cosmo/RadioBeam/fgndsub.cc (modified) (5 diffs)

trunk/Cosmo/RadioBeam/fgndsub.cc

@@ -33 +33 @@
   beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_);
   cout << " ForegroundCleaner::BeamCorrections() done " << endl;
+}
+
+/* --Methode-- */
+int ForegroundCleaner::CleanNegatives(TF seuil)
+{
+  sa_size_t nneg = 0.;
+  for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++) 
+    for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++) 
+      for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) 
+        if (skycube_(kx, ky, kz) < seuil)  {
+          nneg++; skycube_(kx, ky, kz)=seuil;
+        }
+  cout << " ForegroundCleaner::CleanNegatives " << nneg << " sky-pixels <" << seuil << " changed to" << seuil << endl;
+  return (int)nneg;
 }
 
@@ -82 +96 @@
 // Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
 // Return_Array : foreground subtracted LSS signal 
-{
   sa_size_t sz[5];   sz[0]=skycube_.SizeX();  sz[1]=skycube_.SizeY();
   synctemp.SetSize(2, sz); 
@@ -92 +105 @@
   // double freq0 : Frequence premier index en k (MHz)
   // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
+  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  
+    vlnf(k)=log((double)k*dfreq_+freq0_);
+
+  sa_size_t nbinfini=0;
+  sa_size_t nbbad=0;
+  sa_size_t imodprt=skycube_.SizeX()/6;
+  sa_size_t jmodprt=skycube_.SizeY()/6;
   for(sa_size_t i=0; i<skycube_.SizeX(); i++) 
     for(sa_size_t j=0; j<skycube_.SizeY(); j++)  {
@@ -97 +117 @@
       s1=sx=sx2=sy=sxy=0.;
       for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
-        double lnf=log((double)k*dfreq_+freq0_);
-        vlnf(k)=lnf;
+        double lnf = vlnf(k);
         double ttot=(r_8)(skycube_(i,j,k));
         if (ttot < 1.e-5) continue;
@@ -108 +127 @@
       double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
       double T0 = exp(beta+alpha*vlnf(0));
-      if ((i%16==0)&&(j%27==0)) 
+
+      bool fgnan = false;
+      if (!isfinite(alpha)||(!isfinite(beta))) {
+        cout << "extractLSSCube[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        alpha=beta=-999.;
+        fgnan = true;  nbinfini++;
+      }
+      else {
+        double axp1 = beta+alpha*vlnf(0);
+        double axp2 = beta+alpha*vlnf(vlnf.Size()-1);
+        
+        if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
+        cout << "extractLSSCube[" << i << "," << j << "] BAD alpha=" << alpha 
+             << " beta=" << beta << " T0=" << T0 << " - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        fgnan = true;  nbbad++;
+        }
+      }
+      if ((i%imodprt==0)&&(j%jmodprt==0)) 
         cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
-             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,125) << ")" << endl; 
+             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
       synctemp(i,j) = T0;
       specidx(i,j) = alpha;
-      for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
-        r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
-        omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+      if (fgnan) {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++)   omap(i,j,k) = 0.;
+      }
+      else {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
+          r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
+          omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+        }
       }
     }
+  cout << " ForegroundCleaner::extractLSSCube() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
   return omap;
 }
 
-}