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

Timestamp:

Aug 4, 2010, 1:59:10 PM (15 years ago)

Author:

ansari

Message:

Amelioration soustraction des avant-plans par l'introduction d'un fit polynome deg 2 sur ln(Temp)=f(ln(freq)), Reza 04/08/2010

File:

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

trunk/Cosmo/RadioBeam/fgndsub.cc

@@ -9 +9 @@
 #include "cubedef.h"
 #include "matharr.h"
+#include "poly.h"
 
 #include "ctimer.h"
@@ -90 +91 @@
 
 /* --Methode-- */
-TArray< TF >  ForegroundCleaner::extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx)
-{
-  Timer tm("ForegroundCleaner::extractLSSCube");
+TArray< TF >  ForegroundCleaner::extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx)
+{
+  Timer tm("ForegroundCleaner::extractLSSCubeP1");
 // Inputs : maplss, mapsyc, freq0, dfreq
 // Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
@@ -130 +131 @@
       bool fgnan = false;
       if (!isfinite(alpha)||(!isfinite(beta))) {
-        cout << "extractLSSCube[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+        cout << "extractLSSCubeP1[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
              << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
         alpha=beta=-999.;
@@ -140 +141 @@
         
         if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
-        cout << "extractLSSCube[" << i << "," << j << "] BAD alpha=" << alpha 
+        cout << "extractLSSCubeP1[" << i << "," << j << "] BAD alpha=" << alpha 
              << " beta=" << beta << " T0=" << T0 << " - (mapsync="
              << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
@@ -147 +148 @@
       }
       if ((i%imodprt==0)&&(j%jmodprt==0)) 
-        cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
+        cout << "extractLSSCubeP1[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
              << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
       synctemp(i,j) = T0;
@@ -161 +162 @@
       }
     }
-  cout << " ForegroundCleaner::extractLSSCube() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
+  cout << " ForegroundCleaner::extractLSSCubeP1() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
   return omap;
 }
 
+/*
+static inline val_polyn2(double alpha, double beta, double gamma, double x)
+{
+  return (beta+alpha*x+gamma*x*x);
+}
+*/
+
+/* --Methode-- */
+TArray< TF >  ForegroundCleaner::extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx)
+{
+  Timer tm("ForegroundCleaner::extractLSSCubeP2");
+// Inputs : maplss, mapsyc, freq0, dfreq
+// 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); 
+  specidx.SetSize(2, sz);
+  TArray<r_4> omap; 
+  omap.SetSize(skycube_, true);
+  Vector vlnf(skycube_.SizeZ());
+  Vector vlnT(skycube_.SizeZ());
+  int nprt = 0;
+  // 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_);
+  vlnf -= vlnf(0);
+  //  cout << " DBG*extractLSSCubeP2 vlnf(0)=" << vlnf(0) << " vlnf(1)=" << vlnf(1) 
+  //       << "vlnf(last)=" << vlnf(vlnf.Size()-1) << endl;
+  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++)  {
+      vlnT = -12.;
+      Poly polyn;
+      for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
+        double lnf = vlnf(k);
+        double ttot=(r_8)(skycube_(i,j,k));
+        if (ttot < 1.e-5) continue;
+        vlnT(k)=log(ttot);
+      }
+      polyn.Fit(vlnf,vlnT,2);
+      double beta = polyn[0];
+      double alpha = polyn[1];
+      double gamma = polyn[2];
+      double T0 = exp(polyn(vlnf(0))); // exp( val_polyn2(alpha, beta, gamma, vlnf(0)) ); 
+
+      bool fgnan = false;
+      if (!isfinite(alpha)||(!isfinite(beta))||(!isfinite(gamma))) {
+        cout << "extractLSSCubeP2[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        alpha=beta=gamma=-999.;
+        fgnan = true;  nbinfini++;
+      }
+      else {
+        double axp1 = polyn(vlnf(0));
+        double axp2 = polyn(vlnf(vlnf.Size()-1));
+        
+        if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
+          cout << "extractLSSCubeP2[" << i << "," << j << "] BAD alpha=" << alpha 
+            //         << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " axp1=" << axp1 << " axp2=" << axp2 
+               << " beta=" << beta << " gamma=" << gamma << " 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 << "extractLSSCubeP2[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha << " gamma=" << gamma  
+             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+      synctemp(i,j) = T0;
+      specidx(i,j) = alpha;
+      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(polyn(vlnf(k))) );
+          omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+        }
+      }
+    }
+  cout << " ForegroundCleaner::extractLSSCubeP2() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
+  return omap;
+}
+