Index: trunk/Cosmo/RadioBeam/README
===================================================================
--- trunk/Cosmo/RadioBeam/README	(revision 3797)
+++ trunk/Cosmo/RadioBeam/README	(revision 3825)
@@ -22,4 +22,5 @@
   - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
   - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
+  - 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 3797)
+++ trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3825)
@@ -53,5 +53,5 @@
     bool fgaddnoise=false;
     if (narg>4) {
-      pixsignoise=atof(arg[3]);
+      pixsignoise=atof(arg[4]);
       fgaddnoise=true;
     }
Index: trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- trunk/Cosmo/RadioBeam/cubedef.h	(revision 3797)
+++ trunk/Cosmo/RadioBeam/cubedef.h	(revision 3825)
@@ -16,9 +16,10 @@
 static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
 
+/* ---  Parametres avec generation Reza from HASLAM + NVSS */
 static double Freq0MHz = 875.;       // premiere frequence 
 static double FreqSizeMHz = 70.4;    // Taille de la carte en MHz (axe Z)
 
-/* z = 0.56, freq_center=910 MHz , dA_comov ~= 2060 Mpc ,  
-       5 arcmin -> 3 Mpc , 0.275 MHz -> 3 Mpc  (H(ze) ~ 94 km/s/Mpc */
+// z = 0.56, freq_center=910 MHz , dA_comov ~= 2060 Mpc ,  
+//       5 arcmin -> 3 Mpc , 0.500 MHz -> 3 Mpc  (H(ze) ~ 94 km/s/Mpc)
 // Taille de cellule pour le calcul du spectre de puissance 3D 
 static double ComovDA = 2060.;
@@ -26,4 +27,18 @@
 static double YCellSizeMpc = 3.;
 static double ZCellSizeMpc = 1.5;
+
+/* --- Parametres pour utilisation des cartes GSM  
+static double Freq0MHz = 820.;       // premiere frequence 
+static double FreqSizeMHz = 128.;    // Taille de la carte en MHz (axe Z)
+
+// z = 0.6, freq_center=884 MHz , dA_comov ~= 2185 Mpc ,  
+//       5 arcmin -> 3 Mpc , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
+// Taille de cellule pour le calcul du spectre de puissance 3D 
+static double ComovDA = 2185.;
+static double XCellSizeMpc = 3.;
+static double YCellSizeMpc = 3.;
+static double ZCellSizeMpc = 3.;
+*/
+
 
 //---  Parametres des lois de puissance en frequence 
Index: trunk/Cosmo/RadioBeam/gsm2cube.cc
===================================================================
--- trunk/Cosmo/RadioBeam/gsm2cube.cc	(revision 3825)
+++ trunk/Cosmo/RadioBeam/gsm2cube.cc	(revision 3825)
@@ -0,0 +1,152 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme d'extraction d'une partie de carte synchrotron 
+  (HASLAM @ 400 MHz) et fabrication d'un cube 3D (angles,fre) 
+    R. Ansari , C. Magneville - Juin 2010 
+
+  Usage: syncube InFitsName Out3DPPFName [Out2DMapName] 
+---------------------------------------------------------------  */
+
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <math.h>
+#include <iostream>
+#include <typeinfo>
+
+#include "tvector.h"
+#include "srandgen.h"
+#include "fioarr.h"
+#include "sopemtx.h"
+#include "pexceptions.h"
+
+#include "randr48.h"
+
+#include "tvector.h"        // Pour l'utilisation des classes TArray, TMatrix , TVector
+#include "matharr.h"
+
+/*
+#include "spherehealpix.h"   // Pour les cartes spheriques pixelisees au format HEALPix 
+#include "spherethetaphi.h"   // Pour les cartes spheriques pixelisees au format Theta-Phi 
+#include "localmap.h"         // Pour les cartes locales 
+#include "mapoperation.h"         // Pour les cartes locales 
+*/
+#include "skymap.h"
+#include "mapoperation.h"         // Pour les cartes locales 
+
+#include "fitsspherehealpix.h"   // Pour les I/O fits de HEALPix
+#include "fitsspherethetaphi.h"   // Pour les I/O fits de SphereThetaPhi
+#include "fitslocalmap.h"         // Pour les I/O fits de LocalMap<T> 
+
+#include "xastropack.h"     // Pour faire les conversions de coordonnees celestes
+
+// Pour l'initialisation des modules 
+#include "tarrinit.h"
+#include "skymapinit.h"     
+#include "fiosinit.h"     
+
+#include "timing.h"
+#include "ctimer.h"
+
+#include "cubedef.h"
+
+//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------   
+int main(int narg, char* arg[])
+{
+  // Sophya modules initialization
+  TArrayInitiator  _inia;
+  SkyMapInitiator  _inis;
+  FitsIOServerInitiator  _inif;
+  //------- AU LIEU DE ------>  SophyaInit();  
+
+  InitTim();   // Initializing the CPU timer
+  Timer tm("gsm2cube");
+
+  if (narg < 5) {
+    cout << "Usage: gsm2cube GSMPPFDirectory GSMMapNumStart GSMMapNumEnd  Out3DPPFName  \n" << endl;
+    return 1;
+  }
+
+  // decodage arguments 
+  string inpath = arg[1];
+  int numapstart = atoi(arg[2]);
+  int numapend = atoi(arg[3]);
+  string outname = arg[4];
+
+  int rc = 91;
+
+  cout << " ====== gsm2cube :   Input GSM map path= " << inpath << " NumStart=" << numapstart 
+       << " End=" << numapend << " OutName=" << outname << endl;
+  try {
+    int nfreqgsm = numapend-numapstart+1;
+    if (nfreqgsm!=NFreq) {
+      cout << " gsm2cube/ERROR: (nfreqgsm=" << nfreqgsm << ") <> NFreq (=" << NFreq << ") -> exit 9" << endl;
+      return 9;
+    }
+
+    double tet0 = Angle(Theta0Degre,Angle::Degree).ToRadian();
+    double phi0 = Angle(Phi0Degre,Angle::Degree).ToRadian();
+    double dtet = Angle(ThetaSizeDegre,Angle::Degree).ToRadian()/(double)NTheta;
+    double dphi = Angle(PhiSizeDegre,Angle::Degree).ToRadian()/(double)NPhi;
+
+    TArray<r_4> ocube(NPhi,NTheta,NFreq);
+    double mjd2000 = MJDfrDate(1,1,2000);   // Modified Julian Date pour 1 Janvier 2000
+
+    cout << "gsm2cube[1] Loop over input GSM spherical maps ... " << endl;
+
+    for(sa_size_t kf=0; kf<NFreq; kf++) {
+      char mapname[512];
+      sprintf(mapname,"%s/gsm820948_%d.ppf",inpath.c_str(),(int)(kf+numapstart));
+      cout << kf << "- Reading GSM map from file " << mapname << " ... " ;
+      PInPersist pin(mapname);
+      SphereHEALPix<r_4> gsm;
+      pin >> gsm;
+      cout << " Freq=" << gsm.Info()["FMHz"] << " MHz " << endl;
+      for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
+	for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
+	  double theta = j*dtet+tet0;
+	  double phi = i*dphi+phi0;
+	  double ra, dec;  // Angle horaire (right ascension) [0...24], dec [-90 90]
+	  ra=Angle(phi).ToDegree()/360.*24.;
+	  dec=90-Angle(theta).ToDegree();
+	  double glong, glat;  // Longitude et latitude galactique - en degres [0...360], [-90...90]
+	  EqtoGal(mjd2000,ra,dec,&glong,&glat);
+	  phi=Angle(glong,Angle::Degree).ToRadian();
+	  theta=Angle(90-glat,Angle::Degree).ToRadian();
+	  ocube(i,j,kf)=gsm(theta,phi);
+	}
+      }
+    }
+
+    double mean, sigma;
+    MeanSigma(ocube, mean, sigma);
+    cout << "gsm2cube[2] ocube : Mean=" << mean << " Sigma=" << sigma << " Sizes: " << endl;
+    ocube.Show(); 
+
+    
+    // On sauve le cube de sortie
+    {
+      cout << " gsm2cube[3]: Saving output cube to -> " << outname << endl;
+      POutPersist poc(outname);
+      poc << ocube;
+    }
+    rc = 0;
+  }
+  catch (PThrowable& exc) {
+    cerr << " gsm2cube.cc catched Exception " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n gsm2cube.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+  }
+  catch (...) {
+    cerr << " gsm2cube.cc catched unknown (...) exception  " << endl; 
+    rc = 78; 
+  } 
+
+  cout << ">>>> gsm2cube ------- FIN ----------- Rc=" << rc << endl;
+  return rc;
+}
+
+
Index: trunk/Cosmo/RadioBeam/makefile
===================================================================
--- trunk/Cosmo/RadioBeam/makefile	(revision 3797)
+++ trunk/Cosmo/RadioBeam/makefile	(revision 3825)
@@ -6,5 +6,5 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe 
+all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube
 
 clean : 
@@ -29,4 +29,7 @@
 syncube : Objs/syncube
 	echo 'makefile : syncube made'
+
+gsm2cube : Objs/gsm2cube
+	echo 'makefile : gsm2cube made'
 
 srcat2cube : Objs/srcat2cube
@@ -95,4 +98,10 @@
 	$(CXXCOMPILE) -o Objs/applobe.o applobe.cc
 
+### programme gsm2cube
+Objs/gsm2cube : Objs/gsm2cube.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/gsm2cube Objs/gsm2cube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
+
+Objs/gsm2cube.o : gsm2cube.cc  $(PKGHLIST)
+	$(CXXCOMPILE) -o Objs/gsm2cube.o gsm2cube.cc
 
 ### les classes / fonctions 
Index: trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.cc	(revision 3797)
+++ trunk/Cosmo/RadioBeam/specpk.cc	(revision 3825)
@@ -46,4 +46,5 @@
 {
   typ_=typ;
+  SetRenormFac();
 }
 
@@ -52,17 +53,18 @@
 {
   wk/=DeuxPI;
+  double retv=1.;
   switch (typ_) 
     {
     case 1:
-      return Pnu1(wk);
+      retv=Pnu1(wk);
       break;
     case 2:
-      return Pnu2(wk);
+      retv=Pnu2(wk);
       break;
     case 3:
-      return Pnu3(wk);
+      retv=Pnu3(wk);
       break;
     case 4:
-      return Pnu4(wk);
+      retv=Pnu4(wk);
       break;
     default :
@@ -83,8 +85,9 @@
   	}
       }
-      return csp;
-      }
-      break;
-    }
+      retv=csp;
+      }
+      break;
+    }
+  return retv*renorm_fac;
 }
 // Return a vector representing the power spectrum (for checking) 
@@ -97,4 +100,14 @@
 }
 
+double SpectralShape::Sommek2Pk(double kmax, int n)
+{
+  double dk=kmax/(double)n;
+  double s=0.;
+  for(int i=1; i<=n; i++) {
+    double ck=(double)i*dk;
+    s += Value(ck)*ck*ck;
+  }
+  return s*dk*4.*M_PI;
+}
 //--------------------------------------------------
 // -- Four2DResponse class : test P(k) class
Index: trunk/Cosmo/RadioBeam/specpk.h
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.h	(revision 3797)
+++ trunk/Cosmo/RadioBeam/specpk.h	(revision 3825)
@@ -30,6 +30,9 @@
   inline  double Value(double wk) { return((*this)(wk)); }
 // Return a vector representing the power spectrum (for checking) 
-  Histo GetPk(int n=256);      	
+  Histo GetPk(int n=256);
+  double Sommek2Pk(double kmax=1000., int n=5000);
+  inline void SetRenormFac(double f=1.) { renorm_fac=f; } 
   int typ_;
+  double renorm_fac;
 };
 
Index: trunk/Cosmo/RadioBeam/srcat2cube.cc
===================================================================
--- trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3797)
+++ trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3825)
@@ -72,5 +72,9 @@
     sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
     sa_size_t idxf = nvss.IndexNom("S1_4");
-    cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf << endl;
+    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;
 
     TArray<r_4> omap(NPhi,NTheta);
@@ -83,4 +87,6 @@
 
     sa_size_t srccnt=0;
+    sa_size_t extendedsrccnt=0;
+
     double meanflx=0.;
     double flxmin=9.e99;
@@ -89,4 +95,6 @@
     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);
@@ -94,4 +102,6 @@
       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;
@@ -100,5 +110,21 @@
       if ((i<0)||(i>=omap.SizeX()))  continue;
       if ((j<0)||(j>=omap.SizeY()))  continue;
-      omap(i,j) += flx;
+      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;
@@ -108,8 +134,11 @@
     cout << "srccat2cube[3]: Output rectangular map computed " << endl;
     meanflx /= (double)srccnt;
-    cout << " SrcCount in map: " << srccnt << " -> meanFlx=" << meanflx << " min=" << flxmin 
+    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;
Index: trunk/Cosmo/RadioBeam/subtractradsrc.cmd
===================================================================
--- trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3797)
+++ trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3825)
@@ -9,5 +9,8 @@
 # 1.a/ Run SimLSS
 csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,1.5 -Z 0.56 -8 1. -n 10000 -O 0,2 -o lssz056 -T 2
-# 1.b/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention 
+# 1.b/ To run SimLSS with GSM map parameters (DeltaFreq=500 MHz)
+csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,3 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
+
+# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention 
 #  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray) 
 #  RadioBeam cubes : the radial (redshift) direction along Z axis of the cube (TArray) 
@@ -27,8 +30,11 @@
   KeepObj(omap); 
 
-print omap
-saveppf omap lsscube.ppf 
+rename omap lsscube
+print lsscube
+# expmeansig lsscube val
+saveppf lsscube lsscube.ppf 
 
 csh> spiapp -term -exec racube.pic 
+
 
 ## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
@@ -69,4 +75,5 @@
 #  with the lobe effect 
 csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2  
+csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.2 0.35  
 csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2  
 
@@ -78,21 +85,44 @@
 # 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
 csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0.
-# 4.d / 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 subpklobe.ppf 0.35 50. 130. 3.
+# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction  
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. 
+# 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.
 
 ### Step 5 / Check the results using spiapp 
+setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
 delobjs *
 openppf fgndpk.ppf 
-disp fgndpk 'logx logy'
-# Adjust the limits
 openppf fgndpklobe.ppf 
-disp fgndpklobe 'same red'
-disp fgndpklobe 'same red nsta'
 openppf fgndpkcorlobe.ppf 
-disp fgndpkcorlobe 'same magenta nsta'
 openppf lsspk.ppf 
-disp lsspk 'same gold nsta'
 openppf lsspklobe.ppf 
-disp lsspklobe 'same green nsta'
+openppf lsspklobewn.ppf 
 openppf subpklobe.ppf 
-disp subpklobe 'same blue nsta' 
+openppf subpkcorlobe.ppf 
+
+openppf subpknolss.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'
+
+disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
+disp fgndpklobe 'same nsta blue'
+disp fgndpkcorlobe 'same nsta skyblue'
+disp lsspk 'same nsta gold'
+disp lsspklobewn 'same nsta siennared' 
+settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
+set lines ( 'Pk[Foreground]'  'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
+set cols ( navyblue blue skyblue gold siennared ) 
+textdrawer lines cols 'font=helvetica,bold,16 frame'
+
+disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
+disp lsspklobewn 'same nsta siennared' 
+disp subpkcorlobe 'same nsta red'
+disp subpknolsslobe '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 ) 
+textdrawer lines cols 'font=helvetica,bold,16 frame'
