Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3829)
@@ -21,5 +21,5 @@
   - applobe.cc : Antenna/array beam effect on 3D sky cube 
   - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
-  - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
+  - srcat2cube.cc : radio source 3D temperature cube from NVSS or North20cm catalogs
   - gsm2cube.cc : Computes 3D temperature cube from a set of GSM (Global Sky Model) healpix PPF maps  
   - tjyk.cc : test de la classe H21Conversions 
Index: /trunk/Cosmo/RadioBeam/calcpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3829)
@@ -36,5 +36,5 @@
 {
   if (narg<3) {
-    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
+    cout << " Usage: calcpk In3DMap_PPFName OutPk_PPFName  [InRenormFactor] [PixNoiseLevel] [Four3D_PPFName]" << endl;
     return 1;
   }
@@ -54,6 +54,13 @@
     if (narg>4) {
       pixsignoise=atof(arg[4]);
-      fgaddnoise=true;
+      if (pixsignoise>1.e-9)  fgaddnoise=true;
     }
+    bool fgsavef3d=false;
+    string outf3dname="four3d.ppf";
+    if (narg>5) {
+      outf3dname=arg[5];
+      fgsavef3d=true;
+    }
+
 
     cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
@@ -82,4 +89,6 @@
     TArray< complex<r_4> > four3d;
     ffts.FFTForward(inmap, four3d);
+    cout << " calcpk[2.b] inmap.Norm2()=" << inmap.Norm2() << " four3d.Norm2()=" << four3d.Norm2() << endl;
+
     tm.Split(" After FFTForward ");
 
@@ -92,4 +101,5 @@
     pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
 
+
     HProf hp = pkc.ComputePk(0.,256);
     {
@@ -97,7 +107,8 @@
       POutPersist po(outname);
       po << hp;
-      outname = "f3d_" + outname;
-      cout << "calcpk[4.b] : writing four3d to PPF file  " << outname << endl;
-      POutPersist pof(outname);
+    }
+    if (fgsavef3d)   {
+      cout << "calcpk[4.b] : writing four3d to PPF file  " << outf3dname << endl;
+      POutPersist pof(outf3dname);
       pof << four3d;
     }
Index: /trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3829)
@@ -11,5 +11,5 @@
 static sa_size_t NFreq = 256;
 
-static double Theta0Degre = 75.;  // -> Delta = +30 deg
+static double Theta0Degre = 65.;  // -> Delta = +30 deg
 static double Phi0Degre = 135.;   // -> alpha = 9h00
 static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
@@ -45,5 +45,5 @@
 //---  Parametres des lois de puissance en frequence 
 static double AmpPL1 = 1.;   // amp max PowerLaw 1 (synchrotron  
-static double PLidx1 = -2.5;  // index de la loi de puissance synchrotron 
+static double PLidx1 = -2.8;  // index de la loi de puissance synchrotron 
 static double sigPLidx1 = 0.1;  // Sigma de la variation (gaussienne) de index1 
 // Amplitude max de la 2eme composante en loi de puissance (tirage plat 0 ... AmpPL2)
Index: /trunk/Cosmo/RadioBeam/radutil.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/radutil.cc	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/radutil.cc	(revision 3829)
@@ -10,4 +10,5 @@
   setFrequency(freq);
   setOmegaPix(opix);
+  setCosmoParam();
 }
 
@@ -28,2 +29,21 @@
   return jy*1e-26 / omegapix_ * lambda_ * lambda_ / 2. / k_Boltzman_Cst ;
 }
+
+
+double H21Conversions::setCosmoParam(double omegamatter, double omegabaryon, double h100, double fracHI)
+{
+  OmegaMatter_=omegamatter;
+  OmegaBaryons_=omegabaryon;
+  h100_=h100;
+  fracHI_=fracHI;
+  OmegaLambda_=1.-OmegaMatter_;  // Je neglige OmegaRadiation
+  return OmegaLambda_;
+}
+
+double H21Conversions::Mean21cmTemperature_mK()
+{
+  double zz = 1.+z_;
+  double cc=zz*zz/sqrt(OmegaMatter_*zz*zz*zz+OmegaLambda_);
+  cc *= ((h100_/0.7)*(OmegaBaryons_/0.04)*(fracHI_/0.1));
+  return (cc*0.57);
+}
Index: /trunk/Cosmo/RadioBeam/radutil.h
===================================================================
--- /trunk/Cosmo/RadioBeam/radutil.h	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/radutil.h	(revision 3829)
@@ -18,4 +18,8 @@
   double toKelvin(double jy);    // Conversion de Jansky en  temperature (Kelvin) 
 
+  double Mean21cmTemperature_mK();      // Temperature moyenne de l'emission a 21 cm en mK
+  inline double T21cm_Kelvin() { return Mean21cmTemperature_mK()/1000.;  }
+  inline double T21cm_mK() { return Mean21cmTemperature_mK();  }
+
   void setFrequency(double nu);  // on definit la frequence en MHz 
   inline void setRedshift(double z)    // on definit le redshift
@@ -31,4 +35,9 @@
   { double cf=Angle(1.,Angle::ArcMin).ToRadian(); omegapix_ = opix*cf*cf; }
 
+  // Definition des parametres cosmologiques utiles pour le calcul de la temperature d'emission a 21 cm
+  // retourne la valeur de OmegaLambda (univers plat)
+  double setCosmoParam(double omegamatter=0.02581, double omegabaryon=0.0441, double h100=0.719, double fracHI=0.02);
+  inline void setFracHI(double fracHI=0.02) { fracHI_=fracHI; }
+
   inline double getRedshift() { return z_; }
   inline double getFrequency() { return freq_; }
@@ -42,4 +51,5 @@
   { double cf=Angle(1.,Angle::Degree).ToRadian(); return omegapix_/cf/cf; }
 
+
   static double SpeedOfLight_Cst;       // Speed of light  m/sec
   static double Freq021cm_Cst;       // Speed of light  m/sec
@@ -50,4 +60,11 @@
   double lambda_;
   double omegapix_;
+
+  // Parametres cosmologiques pour calcul du coefficient de conversion Mass to T21 
+  double OmegaMatter_;
+  double OmegaBaryons_;
+  double OmegaLambda_;
+  double h100_; 
+  double fracHI_; 
 };
 
Index: /trunk/Cosmo/RadioBeam/srcat2cube.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3829)
@@ -34,4 +34,8 @@
 
 #include "cubedef.h"
+
+//----------------
+int nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& cube);
+int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& cube);
 
 //----------------------------------------------------------------------------
@@ -48,6 +52,6 @@
   Timer tm("srcat2cube");
 
-  if (narg < 3) {
-    cout << "Usage: srccat2cube NVSS_CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
+  if (narg < 4) {
+    cout << "Usage: srccat2cube -nvss/-north20 CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
     return 1;
   }
@@ -55,134 +59,37 @@
 
   // decodage arguments 
-  string outname = arg[2];
-  string inname = arg[1];
+  
+  string copt=arg[1];
+  string outname=arg[3];
+  string inname=arg[2];
   int rc = 91;
 
-  cout << " ====== srccat2cube :   Input NVSS catalog name= " << inname << " OutName=" << outname;
+  cout << " ====== srccat2cube : Input catalog name= " << inname << " OutName=" << outname;
   bool fginmap=true;
   try {
-    DataTable nvss;
-    cout << "srccat2cube[1]: reading NVSS catalog from " << inname << endl;
+    DataTable cat;
+    cout << "srccat2cube[1]: reading source catalog from " << inname << endl;
     {
       FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
-      fis >> nvss;
-    }
-    cout << nvss;
-    sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
-    sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
-    sa_size_t idxf = nvss.IndexNom("S1_4");
-    sa_size_t idxmajax = nvss.IndexNom("MajAxis");
-    sa_size_t idxminax = nvss.IndexNom("MinAxis");
-
-    cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
-	 << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
-
+      fis >> cat;
+    }
+    cout << cat ;
     TArray<r_4> omap(NPhi,NTheta);
-    double tet0 = Theta0Degre; 
-    double phi0 = Phi0Degre; 
-    double tetmax = tet0+ThetaSizeDegre;
-    double phimax = phi0+PhiSizeDegre;
-
-    cout << "srccat2cube[2]: projecting sources to map ..." << endl;
-
-    sa_size_t srccnt=0;
-    sa_size_t extendedsrccnt=0;
-
-    double meanflx=0.;
-    double flxmin=9.e99;
-    double flxmax=-9.e99;
-
-    double dtet = ThetaSizeDegre/(double)NTheta;
-    double dphi = PhiSizeDegre/(double)NPhi;
-    double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
-
-    for (sa_size_t n=0; n<nvss.NRows(); n++)  {
-      r_8* pline=nvss.GetLineD(n);
-      double alpha=pline[idxa];  // alpha en degre
-      double delta=pline[idxd];  // delta en degre 
-      double flx=pline[idxf]*1.e-3;  // flux en Jy
-      double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.;  // taille (extension de la source en arcmin 
-      if (srcszarcmin<1.) srcszarcmin=1.;
-      double tet = 90.-delta;
-      double phi = alpha;
-      sa_size_t i = (phi-phi0)/dphi;
-      sa_size_t j = (tet-tet0)/dtet;
-      if ((i<0)||(i>=omap.SizeX()))  continue;
-      if ((j<0)||(j>=omap.SizeY()))  continue;
-      double srat = (4.*srcszarcmin*srcszarcmin)/(mpixsizarcmin*mpixsizarcmin);
-      if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
-	omap(i,j) += flx*srat;
-      }
-      else {  // on repartit l'energie de la source dans plusieurs pixels
-	extendedsrccnt++;
-	for(int bi=-1;bi<=1;bi++) {
-	  for(sa_size_t bj=-1; bj<=1; bj++) {
-            sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
-	    sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
-	    if ((ii<0)||(ii>=omap.SizeX()))  continue;
-	    if ((jj<0)||(jj>=omap.SizeY()))  continue;	    
-	    if ((bi==0)&&(bj==0))  omap(ii,jj) += flx*srat*0.3;
-	    else omap(ii,jj) += flx*srat*0.7/8.;
-	  }
-	}
-      }
-      srccnt++;   meanflx+=flx;
-      if (flx<flxmin) flxmin=flx;
-      if (flx>flxmax) flxmax=flx;
-    }
-
-    cout << "srccat2cube[3]: Output rectangular map computed " << endl;
-    meanflx /= (double)srccnt;
-    cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt 
-	 << " -> meanFlx=" << meanflx << " min=" << flxmin 
-	 << " max=" << flxmax << " Jy" << endl;
-
-    double mean, sigma;
-    r_4 mintemp, maxtemp;
-    omap.MinMax(mintemp, maxtemp);
-    MeanSigma(omap, mean, sigma);
-    cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
-    omap.Show(); 
-
-    H21Conversions conv;
-    conv.setRedshift(0.);
-    conv.setOmegaPixDeg2(dphi*dtet);
-    cout << "srccat2cube[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
-	 << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
-    omap *= (r_4)conv.toKelvin(1.);
-    MeanSigma(omap, mean, sigma);
-    cout << " After conversion  : Mean=" << mean << " Sigma=" << sigma << "  Kelvin " << endl;
+    TArray<r_4> ocube(NPhi,NTheta,NFreq);
+    if (copt=="-nvss")  nvssTocube(cat, omap, ocube);
+    else north20Tocube(cat, omap, ocube);
     
-    if (narg > 3) {
-      string ppfname = arg[3];
-      cout << " srccat2cube[4]: Saving inmap/outmap tp PPF file-> " << ppfname << endl;
-      POutPersist po(ppfname);
-      po << PPFNameTag("omap") << omap;
-    }
-
-    TArray<r_4> ocube(NPhi,NTheta,NFreq);
-
-    double infreq = 1420.; //  frequence de reference du flux des sources 
-    double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
-    double dfreq = FreqSizeMHz/(double)NFreq;  
-
-    ThSDR48RandGen rg;
-    for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
-      for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
-	double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
-	for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
-	  double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
-	  ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
-	}
-      }
-    }
-    
-    // On sauve le cube de sortie
-    {
-      cout << " srccat2cube[5]: Saving output cube to -> " << outname << endl;
+    {  // On sauve le cube de sortie
+      cout << " srccat2cube[7]: Saving output cube to -> " << outname << endl;
       POutPersist poc(outname);
       poc << ocube;
     }
-
+    
+    if (narg > 4) {
+      string ppfname = arg[4];
+      cout << " srccat2cube[8]: saving 2D source map to PPF file-> " << ppfname << endl;
+      POutPersist po(ppfname);
+      po << omap;
+    }
     rc = 0;
   }
@@ -206,2 +113,211 @@
 
 
+
+/* -- Fonction -- */
+int nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& ocube)
+{
+  sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
+  sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
+  sa_size_t idxf = nvss.IndexNom("S1_4");
+  sa_size_t idxmajax = nvss.IndexNom("MajAxis");
+  sa_size_t idxminax = nvss.IndexNom("MinAxis");
+  
+  cout << " NVSS catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
+       << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
+  
+  double tet0 = Theta0Degre; 
+  double phi0 = Phi0Degre; 
+  double tetmax = tet0+ThetaSizeDegre;
+  double phimax = phi0+PhiSizeDegre;
+  
+  cout << "srccat2cube/NVSS[2]: projecting sources to map ..." << endl;
+  
+  sa_size_t srccnt=0;
+  sa_size_t extendedsrccnt=0;
+  
+  double meanflx=0.;
+  double flxmin=9.e99;
+  double flxmax=-9.e99;
+  
+  double dtet = ThetaSizeDegre/(double)NTheta;
+  double dphi = PhiSizeDegre/(double)NPhi;
+  double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
+  
+  for (sa_size_t n=0; n<nvss.NRows(); n++)  {
+    r_8* pline=nvss.GetLineD(n);
+    double alpha=pline[idxa];  // alpha en degre
+    double delta=pline[idxd];  // delta en degre 
+    double flx=pline[idxf]*1.e-3;  // flux en Jy
+    double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.;  // taille (extension de la source en arcmin 
+    if (srcszarcmin<1.) srcszarcmin=1.;
+    double tet = 90.-delta;
+    double phi = alpha;
+    sa_size_t i = (phi-phi0)/dphi;
+    sa_size_t j = (tet-tet0)/dtet;
+    if ((i<0)||(i>=omap.SizeX()))  continue;
+    if ((j<0)||(j>=omap.SizeY()))  continue;
+    if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
+      omap(i,j) += flx;
+    }
+    else {  // on repartit l'energie de la source dans plusieurs pixels
+      extendedsrccnt++;
+      for(int bi=-1;bi<=1;bi++) {
+	for(sa_size_t bj=-1; bj<=1; bj++) {
+	  sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
+	  sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
+	  if ((ii<0)||(ii>=omap.SizeX()))  continue;
+	  if ((jj<0)||(jj>=omap.SizeY()))  continue;	    
+	  if ((bi==0)&&(bj==0))  omap(ii,jj) += flx*0.3;
+	  else omap(ii,jj) += flx*0.7/8.;
+	}
+      }
+    }
+    srccnt++;   meanflx+=flx;
+    if (flx<flxmin) flxmin=flx;
+    if (flx>flxmax) flxmax=flx;
+  }
+  
+  cout << "srccat2cube/NVSS[3]: Output rectangular map computed " << endl;
+  meanflx /= (double)srccnt;
+  cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt 
+       << " -> meanFlx=" << meanflx << " min=" << flxmin 
+       << " max=" << flxmax << " Jy" << endl;
+  
+  double mean, sigma;
+  r_4 minjy, maxjy;
+  omap.MinMax(minjy, maxjy);
+  MeanSigma(omap, mean, sigma);
+  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky ;  Sizes:" << endl;
+  omap.Show(); 
+  
+  H21Conversions conv;
+  conv.setRedshift(0.);
+  conv.setOmegaPixDeg2(dphi*dtet);
+  cout << "srccat2cube/NVSS[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
+       << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
+  omap *= (r_4)conv.toKelvin(1.);
+  MeanSigma(omap, mean, sigma);
+  r_4 minT, maxT;
+  omap.MinMax(minT, maxT);
+  cout << " NVSS/ After conversion  : Mean=" << mean << " Sigma=" << sigma << " Min=" << minT 
+       << " Max=" << maxT << "  Kelvin " << endl;
+
+  double infreq = 1420.; //  frequence de reference du flux des sources 
+  double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
+  double dfreq = FreqSizeMHz/(double)NFreq;  
+
+  ThSDR48RandGen rg;
+  for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
+    for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
+      double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
+      for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
+	double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
+	ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
+      }
+    }
+  }
+  cout << "srccat2cube/NVSS[5] data cube created from sources " << endl;
+  ocube.MinMax(minT, maxT);
+  MeanSigma(ocube, mean, sigma);
+  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
+ 
+  return srccnt;
+}
+
+/* -- Fonction -- */
+int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& ocube)
+{
+
+  sa_size_t idxa = nor.IndexNom("ra");
+  sa_size_t idxd = nor.IndexNom("dec");
+  sa_size_t idxf = nor.IndexNom("flux");
+  sa_size_t idxslo = nor.IndexNom("SpLO");
+  sa_size_t idxshi = nor.IndexNom("SpHI");
+  
+  cout << " North20cm catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf 
+       << " SpLO: " << idxslo << " SpHI: " << idxshi << endl;
+  
+  double tet0 = Theta0Degre; 
+  double phi0 = Phi0Degre; 
+  double tetmax = tet0+ThetaSizeDegre;
+  double phimax = phi0+PhiSizeDegre;
+  
+  cout << "srccat2cube/North20[2]: projecting sources to map ..." << endl;
+  
+  sa_size_t srccnt=0;
+  sa_size_t lowoksrccnt=0;
+  
+  double meanflx=0.;
+  double flxmin=9.e99;
+  double flxmax=-9.e99;
+  
+  double dtet = ThetaSizeDegre/(double)NTheta;
+  double dphi = PhiSizeDegre/(double)NPhi;
+
+  double infreq = 1420.; //  frequence de reference du flux des sources 
+  double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
+  double dfreq = FreqSizeMHz/(double)NFreq;  
+  
+  for (sa_size_t n=0; n<nor.NRows(); n++)  {
+    r_8* pline=nor.GetLineD(n);
+    double alpha=pline[idxa];  // alpha en degre
+    double delta=pline[idxd];  // delta en degre 
+    double flx=pline[idxf]*1.e-3;  // flux en Jy
+    double tet = 90.-delta;
+    double phi = alpha*360./24.;
+    sa_size_t i = (phi-phi0)/dphi;
+    sa_size_t j = (tet-tet0)/dtet;
+    if ((i<0)||(i>=omap.SizeX()))  continue;
+    if ((j<0)||(j>=omap.SizeY()))  continue;
+    omap(i,j) += flx;
+    srccnt++;   meanflx+=flx;
+    if (flx<flxmin) flxmin=flx;
+    if (flx>flxmax) flxmax=flx;
+    double slo=pline[idxslo];
+    if (slo<9.) {  // source detected at 80 cm 
+      lowoksrccnt++;
+    }
+    else slo=5.;
+    for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
+      double rapfreq = pow((freq0+k*dfreq)/infreq, slo);
+      ocube(i,j,k) += flx*rapfreq;
+    }
+  }
+  
+  cout << "srccat2cube/North20[3]: Output rectangular map computed " << endl;
+  meanflx /= (double)srccnt;
+  cout << " SrcCount in map: " << srccnt << " SpLowOK=" << lowoksrccnt 
+       << " -> meanFlx=" << meanflx << " min=" << flxmin 
+       << " max=" << flxmax << " Jy" << endl;
+  
+  double mean, sigma;
+  r_4 minjy, maxjy;
+  omap.MinMax(minjy, maxjy);
+  MeanSigma(omap, mean, sigma);
+  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
+
+  ocube.MinMax(minjy, maxjy);
+  MeanSigma(ocube, mean, sigma);
+  cout << " Cube : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
+
+  H21Conversions conv;
+  conv.setRedshift(0.);
+  conv.setOmegaPixDeg2(dphi*dtet);
+  cout << "srccat2cube/North20[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad"  
+       << " @1400MHz  toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
+  
+  // Jansky to Kelvin conversion 
+  for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
+    conv.setFrequency(freq0+k*dfreq);
+    //    cout << " DBG* Freq= " << freq0+k*dfreq << " -> toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
+    ocube(Range::all(), Range::all(), Range(k)) *= (r_4)conv.toKelvin(1.);
+  }
+  cout << "srccat2cube/North20[5] data cube in Kelvin computed " << endl;
+
+  r_4 minT, maxT;
+  ocube.MinMax(minT, maxT);
+  MeanSigma(ocube, mean, sigma);
+  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
+  
+  return srccnt;
+}
Index: /trunk/Cosmo/RadioBeam/subtractradsrc.cmd
===================================================================
--- /trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3829)
@@ -46,8 +46,8 @@
 csh> cat > sumcubes.pic
 openppf syncube.ppf
-openppf nvsscube.ppf
+openppf srcnor3d.ppf
 # expmeansig syncube val 
-# expmeansig nvsscube val
-c++exec TArray<r_4> fgndcube = syncube+nvsscube; KeepObj(fgndcube);
+# expmeansig srcnor3d val
+c++exec TArray<r_4> fgndcube = syncube+srcnor3d; KeepObj(fgndcube);
 print fgndcube
 # expmeansig fgndcube val
@@ -89,4 +89,7 @@
 # 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130 
 csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3.
+# 4.f / Estimate residual noise from Foreground removal : 
+csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3.
+csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.
 
 ### Step 5 / Check the results using spiapp 
@@ -103,8 +106,15 @@
 
 openppf subpknolss.ppf
+openppf subpknolssnocor.ppf
 
 disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-11,8e-6 gold'
 disp lsspklobe 'same nsta orange'
 disp lsspklobewn 'same nsta siennared'
+settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black' 
+
+#  Calcul du volume total en Mpc^3
+set VOL 3*3*3*360*360*256
+plot2d lsspklobe x val*$VOL 1 'same nsta orange'
+plot2d lsspklobe x val*$VOL 1 'same nsta orange'
 
 disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
@@ -121,8 +131,28 @@
 disp lsspklobewn 'same nsta siennared' 
 disp subpkcorlobe 'same nsta red'
-disp subpknolsslobe 'same nsta green'
+disp subpknolss 'same nsta green'
 
-settitle 'Recovered Pk[LSS] from LSS+Foreground(GSM) (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
-set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[residual,NoLSS]' )
-set cols ( gold siennared red green ) 
+#  Calcul du volume total en Mpc^3
+set VOL 3*3*3*360*360*256
+plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
+plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared' 
+plot2d subpkcorlobe x val*$VOL 1 'same nsta cpts marker=box,5 red' 
+plot2d subpklobe x val*$VOL 1 'same nsta cpts marker=box,5 blueviolet' 
+plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 
+
+settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
+setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
+set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[ExtLSS,NoBeamCor]' 'Pk[residual,NoLSS]' )
+set cols ( gold siennared red blueviolet green ) 
 textdrawer lines cols 'font=helvetica,bold,16 frame'
+
+
+plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
+plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 red' 
+plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 
+plot2d subpknolssnocor x val*$VOL 1 'same nsta cpts marker=box,5 magenta' 
+setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
+settitle 'Recovered Pk[LSS] and residual systematics' ' ' 'font=helvetica,bold,18'
+set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[residual,NoLSS]' 'Pk[residual,NoLSS,NoBeamCorrection]' )
+set cols ( gold red green magenta ) 
+textdrawer lines cols 'font=helvetica,bold,16 frame'
Index: /trunk/Cosmo/RadioBeam/tjyk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3828)
+++ /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3829)
@@ -25,4 +25,5 @@
        << " MHz , Lambda=" << conv.getLambda() << " m , OmegaPix=" << conv.getOmegaPix() << " srad" << endl;
   cout << "  toKelvin(1 Jy)= " << conv.toKelvin(1.) << "    toJansky(1 K)=" << conv.toJansky(1.) << endl;
+  cout << "  Mean 21cm emission temperature= " << conv.T21cm_mK() << " milliKelvin " << endl;
   cout << " --------------------------------------------------------------- " << endl;
 
