Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3784)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3785)
@@ -19,4 +19,5 @@
   - calcpk2.cc : LSS deltaT/T P(k) computation from LSS+foreground cubes after cleaning
   - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
+  - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
   - tjyk.cc : test de la classe H21Conversions 
 
@@ -25,4 +26,5 @@
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
   - radutil.h radutil.cc : utilitaire de conversion a 21 cm
+  - cubedef.h : size and limits for 3D temperature cubes generated by  syncube.cc , srcat2cube.cc
  
 E/ P(k) in temperature plot + "cosmic variance" + noise ...
Index: /trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3785)
+++ /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3785)
@@ -0,0 +1,33 @@
+//  Definition de la carte synchrotron/radio sources ...
+//  Projet BAORadio , Juin 2010 
+
+#ifndef CUBEDEF_SEEN
+#define CUBEDEF_SEEN
+
+// Definition tailles et bornes du cube angles (X,Y), frequence (Z)
+static sa_size_t NTheta=256;
+static sa_size_t NPhi=256;
+static sa_size_t NFreq = 128;
+
+static double Theta0Degre = 60.;  // -> Delta = +30 deg
+static double Phi0Degre = 120.;   // -> alpha = 8h00
+static double ThetaSizeDegre = 60.;  // Taille de la carte en degre selon alpha (axe X)
+static double PhiSizeDegre = 60.;    // Taille de la carte en degre selon delta (axe Y)
+
+static double Freq0MHz = 840.;       // premiere frequence 
+static double FreqSizeMHz = 128.;    // Taille de la carte en MHz (axe Z)
+
+
+//---  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 sigPLidx1 = 0.1;  // Sigma de la variation (gaussienne) de index1 
+// Amplitude max de la 2eme composante en loi de puissance (tirage plat 0 ... AmpPL2)
+static double AmpPL2 = 0.1;   // amp max PowerLaw 2 
+static double PLidx2 = -3.2; 
+static double sigPLidx2 = 0.15;
+
+static double PLidxSrc = -2.2;  // index de la loi de puissance des sources  
+static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
+
+#endif
Index: /trunk/Cosmo/RadioBeam/makefile
===================================================================
--- /trunk/Cosmo/RadioBeam/makefile	(revision 3784)
+++ /trunk/Cosmo/RadioBeam/makefile	(revision 3785)
@@ -6,5 +6,5 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise calcpk calcpk2 syncube tjyk
+all : pknoise calcpk calcpk2 syncube srcat2cube tjyk
 
 clean : 
@@ -23,4 +23,7 @@
 syncube : Objs/syncube
 	echo 'makefile : syncube made'
+
+srcat2cube : Objs/srcat2cube
+	echo 'makefile : srcat2cube made'
 
 tjyk : Objs/tjyk
@@ -55,5 +58,12 @@
 	$(CXXCOMPILE) -o Objs/syncube.o syncube.cc
 
-### programme syncube
+### programme srcat2cube
+Objs/srcat2cube : Objs/srcat2cube.o 
+	$(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o  Objs/radutil.o $(SOPHYAEXTSLBLIST)
+
+Objs/srcat2cube.o : srcat2cube.cc radutil.h 
+	$(CXXCOMPILE) -o Objs/srcat2cube.o srcat2cube.cc
+
+### programme tjyk
 Objs/tjyk : Objs/tjyk.o Objs/radutil.o 
 	$(CXXLINK) -o Objs/tjyk Objs/tjyk.o Objs/radutil.o $(SOPHYAEXTSLBLIST)
Index: /trunk/Cosmo/RadioBeam/srcat2cube.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3785)
+++ /trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3785)
@@ -0,0 +1,178 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de fabrication d'un cube 3D (angles,fre) 
+  a partir du catalogue de source radio (NVSS) 
+    R. Ansari , C. Magneville - Juin 2010 
+
+  Usage: srccat2cube CatalogFitsName Out3DPPFName [Out2DMapName]
+---------------------------------------------------------------  */
+
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <math.h>
+#include <iostream>
+#include <typeinfo>
+
+#include "array.h"
+#include "histats.h"
+
+#include "swfitsdtable.h"
+#include "fitshdtable.h"
+
+#include "randr48.h"
+
+#include "xastropack.h"     // Pour faire les conversions de coordonnees celestes
+
+#include "radutil.h"
+
+// Pour l'initialisation des modules 
+#include "tarrinit.h"
+#include "histinit.h"
+#include "fiosinit.h"     
+
+#include "timing.h"
+#include "ctimer.h"
+
+#include "cubedef.h"
+
+//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------   
+int main(int narg, char* arg[])
+{
+  // Sophya modules initialization
+  TArrayInitiator  _inia;
+  HiStatsInitiator  _inih;
+  FitsIOServerInitiator  _inif;
+  //------- AU LIEU DE ------>  SophyaInit();  
+
+  InitTim();   // Initializing the CPU timer
+  Timer tm("srcat2cube");
+
+  if (narg < 3) {
+    cout << "Usage: srccat2cube NVSS_CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
+    return 1;
+  }
+
+
+  // decodage arguments 
+  string outname = arg[2];
+  string inname = arg[1];
+  int rc = 91;
+
+  cout << " ====== srccat2cube :   Input NVSS catalog name= " << inname << " OutName=" << outname;
+  bool fginmap=true;
+  try {
+    DataTable nvss;
+    cout << "srccat2cube[1]: reading NVSS 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");
+    cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf << endl;
+
+    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;
+    double meanflx=0.;
+    double flxmin=9.e99;
+    double flxmax=-9.e99;
+
+    double dtet = ThetaSizeDegre/(double)NTheta;
+    double dphi = PhiSizeDegre/(double)NPhi;
+    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 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;
+      omap(i,j) += flx;
+      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 << " -> meanFlx=" << meanflx << " min=" << flxmin 
+	 << " max=" << flxmax << " Jy" << endl;
+
+    double mean, sigma;
+    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;
+    
+    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;
+      POutPersist poc(outname);
+      poc << ocube;
+    }
+
+    rc = 0;
+  }
+  catch (PThrowable& exc) {
+    cerr << " srccat2cube.cc catched Exception " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n srccat2cube.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+  }
+  catch (...) {
+    cerr << " srccat2cube.cc catched unknown (...) exception  " << endl; 
+    rc = 78; 
+  } 
+
+  cout << ">>>> srccat2cube[9] ------- FIN ----------- Rc=" << rc << endl;
+  return rc;
+}
+
+
Index: /trunk/Cosmo/RadioBeam/syncube.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3784)
+++ /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3785)
@@ -44,7 +44,8 @@
 #include "fiosinit.h"     
 
-
 #include "timing.h"
 #include "ctimer.h"
+
+#include "cubedef.h"
 
 //----------------------------------------------------------------------------
@@ -65,13 +66,4 @@
     return 1;
   }
-
-  //---  Parametres des lois de puissance en frequence 
-  double AmpPL1 = 1.;   // amp max PowerLaw 2 
-  double PLidx1 = -2.5;  // index de la loi de puissance synchrotron 
-  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)
-  double AmpPL2 = 0.1;   // amp max PowerLaw 2 
-  double PLidx2 = -3.2; 
-  double sigPLidx2 = 0.15;
 
   // decodage arguments 
@@ -101,7 +93,5 @@
     cout << "syncube[1]: Input resolution doubled : " << inmap << endl;
 
-    sa_size_t NTet=256;
-    sa_size_t NPhi=256;
-    LocalMap<r_4> outmap(NPhi,NTet,60.,60.,110.,150.);
+   LocalMap<r_4> outmap(NPhi,NTheta,60.,60.,110.,150.);
     for(int_4 kk=0; kk<outmap.NbPixels(); kk++) {
       double theta, phi;  // Theta, Phi en radians
@@ -110,9 +100,9 @@
     }    
 
-    TArray<r_4> omap(NPhi,NTet);
-    double tet0 = Angle(60.,Angle::Degree).ToRadian();
-    double phi0 = Angle(120.,Angle::Degree).ToRadian();
-    double dtet = Angle(60.,Angle::Degree).ToRadian()/(double)NTet;
-    double dphi = Angle(60.,Angle::Degree).ToRadian()/(double)NPhi;
+    TArray<r_4> omap(NPhi,NTheta);
+    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;
     for (sa_size_t j=0; j<omap.SizeY(); j++)  {
       double theta = j*dtet+tet0;
@@ -122,7 +112,12 @@
       }
     }
+
+    double mean, sigma;
+    MeanSigma(omap, mean, sigma);
+    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
+    omap.Show(); 
+
     //    Sph2Sph(inmap,outmap);
     
-    cout << "syncube[2]: Output rectangular map computed " << outmap << endl;
     if (narg > 3) {
       string ppfname = arg[3];
@@ -134,10 +129,9 @@
     }
 
-    sa_size_t NFreq = 128;
-    TArray<r_4> ocube(NPhi,NTet,NFreq);
+    TArray<r_4> ocube(NPhi,NTheta,NFreq);
 
     double infreq = 400.; // Frequence carte input en MHz
-    double freq0 = 840.;  // Freq0 du cube de sortie 
-    double dfreq = 1.;  
+    double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
+    double dfreq = FreqSizeMHz/(double)NFreq;  
 
     ThSDR48RandGen rg;
@@ -162,5 +156,5 @@
     // On sauve le cube de sortie
     {
-      cout << " syncube[2]: Saving output cube to -> " << outname << endl;
+      cout << " syncube[4]: Saving output cube to -> " << outname << endl;
       POutPersist poc(outname);
       poc << ocube;
@@ -170,5 +164,5 @@
   }
   catch (PThrowable& exc) {
-    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
+    cerr << " syncube.cc catched Exception " << exc.Msg() << endl;
     rc = 77;
   }  
Index: /trunk/Cosmo/RadioBeam/tjyk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3784)
+++ /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3785)
@@ -2,4 +2,6 @@
 #include <iostream>
 
+// Appel : tjyk [z=redshift] [resol_arcmin]
+ 
 int main(int narg, char* arg[])
 {
@@ -7,7 +9,13 @@
   H21Conversions conv;
   
-  conv.setRedshift(0.6);
-  conv.setOmegaPixArcmin2(10.*10.);
-  cout << " --------------------------------------------------------------- " << endl;
+  double z = 0.;
+  double resol = 10.;
+  if (narg>1) z = atof(arg[1]);
+  if (narg>2) resol = atof(arg[2]);
+
+  conv.setRedshift(z);
+  conv.setOmegaPixArcmin2(resol*resol);
+  cout << " ------------------------------ H21Conversions ------------------------- " << endl;
+  cout << "  Redshift=z= " << z << " Resol(arcmin)= " << resol << endl;
   cout << " H21Conversions: z=" << conv.getRedshift() << " Freq=" << conv.getFrequency()
        << " MHz , Lambda=" << conv.getLambda() << " m , OmegaPix=" << conv.getOmegaPix() << " srad" << endl;
