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

Timestamp:

Sep 28, 2011, 5:13:51 PM (14 years ago)

Author:

ansari

Message:

modifs pour pouvoir imposer la moyenne en temp des plans X,Y des cubes lors de l'extraction du signal HI, Reza 28/9/2011

File:

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

trunk/Cosmo/RadioBeam/fgndsub.cc

@@ -12 +12 @@
 
 #include "ctimer.h"
+
+/* --Methode-- */
+PowerLawChecker::PowerLawChecker(TArray< TF >& skycube)
+  : skycube_(skycube)
+{
+}
+/* --Methode-- */
+void PowerLawChecker::CheckXYMean()
+{
+  // double freq0 : Frequence premier index en k (MHz)
+  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
+  double freq0_ = Freq0MHz;
+  double dfreq_ = FreqSizeMHz/(double)NFreq;
+
+  double tempfirst,templast;
+  r_8 s1, sx, sx2, sy, sxy;
+  s1=sx=sx2=sy=sxy=0.;
+  double lnf0=0.;
+  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
+    double lnf=log((double)k*dfreq_+freq0_);
+    double ttot = Mean(skycube_(Range::all(), Range::all(), Range(k)));
+    if (k==0) { tempfirst=ttot;  lnf0=lnf; }
+    if (k==skycube_.SizeZ()-1) templast=ttot;
+    if (ttot < 1.e-5) continue;
+    double lntt=log(ttot);
+    s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
+    sy+=lntt;    sxy+=(lnf*lntt);
+  }
+  double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
+  double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
+  double T0 = exp(beta+alpha*lnf0);
+  
+  cout << " PowerLawChecker::CheckMean() meanTemp(0 ...last) " << tempfirst << " ... " 
+       << templast << endl;
+  bool fgnan = false;
+  if (!isfinite(alpha)||(!isfinite(beta))) {
+    cout << "ePowerLawChecker::CheckMean()  Not finite alpha, beta " << endl;
+    fgnan = true; 
+  }
+  cout << "PowerLawChecker::CheckMean() - T0=" << T0 << " alpha=" << alpha << "(beta=" << beta << ")" << endl; 
+  return;
+}
 
 /* --Methode-- */
@@ -34 +76 @@
   beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_, maxratio_);
   cout << " ForegroundCleaner::BeamCorrections() done " << endl;
+}
+
+/* --Methode-- */
+int ForegroundCleaner::FixMeanXYTemp(double T0, double alpha)
+{
+  cout << "ForegroundCleaner::FixMeanXYTemp(T0=" << T0 << ",alpha=" << alpha << ")" << endl; 
+  double lnf0=log(freq0_);
+  sa_size_t modprt=skycube_.SizeZ()/12;
+  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
+    double lnf=log((double)k*dfreq_+freq0_);
+    double fittedtemp = T0*exp(alpha*(lnf-lnf0));
+    TArray<TF> slice = skycube_(Range::all(), Range::all(), Range(k));
+    TF deltatemp = (TF)(fittedtemp-(double)Mean(slice));
+    slice += deltatemp;
+    if (k%modprt == 0) 
+      cout << "FixMeanXYTemp[k=" << k << " MeanXYTemp=" << fittedtemp-deltatemp << " -> " << fittedtemp 
+           << " (DeltaTemp=" << deltatemp << ")" << endl;
+  }
+  cout << "ForegroundCleaner::FixMeanXYTemp done" << endl; 
+  return 0;
 }