Index: trunk/Cosmo/RadioBeam/applobe.cc
===================================================================
--- trunk/Cosmo/RadioBeam/applobe.cc	(revision 4014)
+++ trunk/Cosmo/RadioBeam/applobe.cc	(revision 4022)
@@ -58,9 +58,10 @@
 
   if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
-    cout << "Usage: applobe [-t -g -fib -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n" 
+    cout << "Usage: applobe [-t -g -fib -kzf/-nzf -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n" 
 	 << "               [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
     if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
       cout << "   -t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
 	   << "   -fib :  Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n" 
+	   << "   -kzf -nzf :  Keep (default) or Not zero space frequency when applying lobes (BeamEffect class) \n" 
 	   << "   -mxr val: Max beam correction factor (default=10.) \n"
            << "   Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n" 
@@ -79,4 +80,5 @@
   bool fixedbeam=false;  // true -> apply freq. independent beam 
   double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe 
+  bool preservefreq0=true; // true -> keep zero frequency when appyling lobe 
 
   // decodage argument optionnel 
@@ -87,4 +89,6 @@
     else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
     else if (fbo=="-fib")  { fixedbeam=true; arg++; narg--; }
+    else if (fbo=="-kzf")  { preservefreq0=true; arg++; narg--; }
+    else if (fbo=="-nzf")  { preservefreq0=false; arg++; narg--; }
     else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
     else fgoptarg=false;
@@ -164,5 +168,5 @@
     else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
 	      << " DoL=" << DIAMETRE/lambda << " ) " << endl;
-    BeamEffect beam(*fresp_p);
+    BeamEffect beam(*fresp_p,preservefreq0);
 
     if (fgcorrbeam) {
Index: trunk/Cosmo/RadioBeam/calcpk2.cc
===================================================================
--- trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 4014)
+++ trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 4022)
@@ -44,9 +44,10 @@
 {
   if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
-    cout << " Usage: [-t -g -mxr val] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
+    cout << " Usage: [-t -g -mxr val -fmt T0,alpha] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
          << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
 	 << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
     if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
       cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
+	   << "-fmt T0,alpha: fix mean slice temp. according to power law T0,alpha \n "
 	   << "-mxr val: Max beam correction factor (default=10.) \n "
 	   << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
@@ -79,4 +80,7 @@
     double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe
 
+    bool fgfix_mXYtemp=false;
+    double T0_pl=1.; 
+    double alpha_pl=0.;
     // decodage argument optionnel 
     bool fgoptarg=true;
@@ -86,4 +90,6 @@
       else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
       else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
+      else if (fbo=="-fmt")  
+	{ fgfix_mXYtemp=true; arg++; sscanf(arg[1],"%lg,%lg",&T0_pl,&alpha_pl); arg++; narg-=2; }
       else fgoptarg=false;
     }
@@ -202,13 +208,17 @@
 
     ForegroundCleaner  cleaner(*arep_p, tbeam, skycube, maxratio); 
+    if (fgfix_mXYtemp) {
+      cout << "calcpk2[3.b] : calling cleaner.FixMeanXYTemp(T0=" << T0_pl << ",alpha=" << alpha_pl << ")" << endl;
+      cleaner.FixMeanXYTemp(T0_pl,alpha_pl);
+    }
     if (fgcorrbeam) {
-      cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam  
+      cout << "calcpk2[3.c] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam  
 	   << " D/Lambda=" << DoL  << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
       cleaner.BeamCorrections();
     }
-    cout << " calcpk2[3.c] : calling cleaner.CleanNegatives() ... " << endl;
+    cout << " calcpk2[3.d] : calling cleaner.CleanNegatives() ... " << endl;
     cleaner.CleanNegatives();
     if (fgclnsrc) {
-      cout << "calcpk2[3.d] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
+      cout << "calcpk2[3.e] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
       cleaner.CleanPointSources(nsigsrc);
     }
Index: trunk/Cosmo/RadioBeam/ckpowl.cc
===================================================================
--- trunk/Cosmo/RadioBeam/ckpowl.cc	(revision 4022)
+++ trunk/Cosmo/RadioBeam/ckpowl.cc	(revision 4022)
@@ -0,0 +1,100 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de verification comportement en loi de puissance 
+    R. Ansari , Sep 2011
+
+  Usage: ckpowl In3DMap  [InRenormFactor] 
+---------------------------------------------------------------  */
+
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include <iostream>
+#include <string>
+#include <math.h>
+
+#include <typeinfo>
+
+#include "specpk.h"
+#include "histats.h"
+
+#include "qhist.h" 
+#include "lobe.h" 
+#include "cubedef.h" 
+#include "fgndsub.h" 
+
+#include "histinit.h"
+#include "fftwserver.h"
+#include "randr48.h"      
+
+#include "ctimer.h"
+
+typedef ThSDR48RandGen RandomGenerator ;
+
+//-------------------------------------------------------------------------
+//      ------------------ MAIN PROGRAM ------------------------------
+//-------------------------------------------------------------------------
+int main(int narg, const char* arg[])
+{
+  if (narg<3) {
+    cout << " Usage: calcpk In3DMap_PPFName OutName [InRenormFactor] [PixNoiseLevel] " << endl;
+    return 1;
+  }
+  Timer tm("calcpk");
+  int rc = 0;
+  try {
+    string inppfname = arg[1];
+    string outname = arg[2];
+    double renfact=1.;
+    bool fgrenfact=false;
+    if (narg>3) {
+      renfact=atof(arg[3]);
+      fgrenfact=true;
+    }
+    double pixsignoise = 0.;
+    bool fgaddnoise=false;
+    if (narg>4) {
+      pixsignoise=atof(arg[4]);
+      if (pixsignoise>1.e-9)  fgaddnoise=true;
+    }
+
+
+    cout << "ckpowl[1] : reading 3D map from file " << inppfname << endl;
+    TArray<r_4> inmap;
+    {
+      PInPersist pin(inppfname);
+      pin >> inmap;
+    }
+    if (fgrenfact) { 
+      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
+      inmap *= renfact;
+    }
+    double mean, sigma;
+    MeanSigma(inmap, mean, sigma);
+    cout << " InMap sizes " << inmap.InfoString() << endl;
+    inmap.Show(); 
+    cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
+    if (fgaddnoise) { 
+      BeamEffect::AddNoise(inmap, pixsignoise);
+    }
+
+    tm.Split(" After read ");
+    cout << "ckpowl[2] : Checking power law behaviour " << endl;
+    PowerLawChecker plck(inmap);
+    plck.CheckXYMean();
+    tm.Split(" Done PowerLawChecker::CheckXYMean() ");    
+  }  // End of try bloc 
+  catch (PThrowable & exc) {  // catching SOPHYA exceptions
+    cerr << " ckpowl.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
+         << "\n...exc.Msg= " << exc.Msg() << endl;
+    rc = 99;
+  }
+  catch (std::exception & e) {  // catching standard C++ exceptions
+    cerr << " ckpowl.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
+    rc = 98;
+  }
+  catch (...) {  // catching other exceptions
+    cerr << " ckpowlx.cc: some other exception (...) was caught ! " << endl;
+    rc = 97;
+  }
+  cout << " ==== End of ckpowl.cc program  Rc= " << rc << endl;
+  return rc;    
+}
Index: trunk/Cosmo/RadioBeam/fgndsub.cc
===================================================================
--- trunk/Cosmo/RadioBeam/fgndsub.cc	(revision 4014)
+++ trunk/Cosmo/RadioBeam/fgndsub.cc	(revision 4022)
@@ -12,4 +12,46 @@
 
 #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,4 +76,24 @@
   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;
 }
 
Index: trunk/Cosmo/RadioBeam/fgndsub.h
===================================================================
--- trunk/Cosmo/RadioBeam/fgndsub.h	(revision 4014)
+++ trunk/Cosmo/RadioBeam/fgndsub.h	(revision 4022)
@@ -23,8 +23,21 @@
 #endif
 
+//------------------------------------------------------------------------------------
+// Cacarterisation du comportement en loi de puissance des cubes de temperature 
+class PowerLawChecker {
+public:
+  PowerLawChecker(TArray< TF >& skycube);
+  void CheckXYMean();
+
+  TArray< TF > skycube_;
+};
+
+//------------------------------------------------------------------------------------
+// Classe implementant la soustraction des avant-plans sous forme de loi de puissance 
 class ForegroundCleaner {
 public:
   ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube, double maxratio=10.);
   void BeamCorrections();
+  int FixMeanXYTemp(double T0, double alpha);
   int CleanNegatives(TF seuil=1.e-6);
   int CleanPointSources(double nsigmas=5.);
@@ -39,7 +52,7 @@
   double maxratio_;
   double dx_, dy_;   // taille des pixels (radians) de skycube
-  double freq0_, dfreq_;    // 1ere frequence et bin en frequence de skycube_;
-  
+  double freq0_, dfreq_;    // 1ere frequence et bin en frequence de skycube_;  
 };
 
+
 #endif
Index: trunk/Cosmo/RadioBeam/makefile
===================================================================
--- trunk/Cosmo/RadioBeam/makefile	(revision 4014)
+++ trunk/Cosmo/RadioBeam/makefile	(revision 4022)
@@ -1,10 +1,10 @@
 ############################################################ 
 ##  makefile for interferometer response 
-##  R.Ansari - Avril 2010  
+##  R.Ansari - Avril 2010 - Sep 2011
 ############################################################ 
 
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube
+all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube ckpowl
 
 clean : 
@@ -41,4 +41,7 @@
 applobe : Objs/applobe
 	echo 'makefile : applobe made'
+
+ckpowl : Objs/ckpowl
+	echo 'makefile : ckpowl made'
 
 ### programme pknoise
@@ -105,4 +108,11 @@
 	$(CXXCOMPILE) -o Objs/gsm2cube.o gsm2cube.cc
 
+### programme ckpowl
+Objs/ckpowl : Objs/ckpowl.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/ckpowl Objs/ckpowl.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
+
+Objs/ckpowl.o : ckpowl.cc $(PKGHLIST)
+	$(CXXCOMPILE) -o Objs/ckpowl.o ckpowl.cc
+
 ### les classes / fonctions 
 Objs/fgndsub.o : fgndsub.cc $(PKGHLIST)
Index: trunk/Cosmo/RadioBeam/sensfgnd21cm.tex
===================================================================
--- trunk/Cosmo/RadioBeam/sensfgnd21cm.tex	(revision 4014)
+++ trunk/Cosmo/RadioBeam/sensfgnd21cm.tex	(revision 4022)
@@ -128,5 +128,5 @@
 cm emission. Such a 3D matter distribution map can be used to test the Cosmological model and to constrain the Dark Energy 
 properties or its equation of state. A novel approach, called intensity mapping can be used to map the \HI distribution, 
-using radio interferometers with large instanteneous field of view and waveband.}
+using radio interferometers with large instantaneous field of view and waveband.}
  % aims heading (mandatory)
   { In this paper, we study the sensitivity of different radio interferometer configurations, or multi-beam 
@@ -352,6 +352,6 @@
 \end{equation}
 Introducing the \HI mass fraction relative to the total baryon mass $\gHI$, the 
-neutral hydrogen number density relative fluctuations can be written as, and the corresponding
-21 cm emission temperature can be written as:
+neutral hydrogen number density and the corresponding 21 cm emission temperature 
+can be written as a function of \HI relative density fluctuations: 
 \begin{eqnarray}
 \etaHI (\vec{\Theta}, z(\lambda) ) & = & \gHIz \times \Omega_B  \frac{\rho_{crit}}{m_{H}}  \times 
