Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3782)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3783)
@@ -1,4 +1,4 @@
 ===============================================================================
-  Projet BAORadio / HSHS
+  Projet BAORadio 
   Simulation programs for radio lobe reconstruction, noise estimate etc ... 
  
@@ -9,26 +9,35 @@
 Note :  Uses SOPHYA class library  (http://www.sophya.org )
 
-A/ (u-v) plane coverage and Noise-P(k) computation for different interferometer configuration
-============================================================================================== 
-List of files :
+A/ Interferometer sensitivity/noise power spectrum and foreground/radio sources
+==============================================================================================
+
+A.1/ (u-v) plane coverage and Noise-P(k) computation for different interferometer configuration
+  - pknoise.cc : main program for Noise-P(k) computation
+  - plpkn.pic : spiapp script for plotting pknoise.cc program output 
+A.2/ power spectrum computation and foreground/radio source 
+  - calcpk.cc : power spectrum P(k) computation from 3D mass or temperature map
+  - 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
+  - tjyk.cc : test de la classe H21Conversions 
+
+A.3/ List of other common files :
   - specpk.h specpk.cc : 3D power spectrum computation classes  
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
-  - pknoise.cc : main program , makefile 
-  - plpkn.pic : spiapp script for plotting pknoise.cc program output 
-   
-B/ P(k) in temperature plot + "cosmic variance" + noise ...
+  - radutil.h radutil.cc : utilitaire de conversion a 21 cm
+ 
+E/ P(k) in temperature plot + "cosmic variance" + noise ...
 =============================================================
 Note : to be executed in spiapp 
 List of Files : anapkn.pic + cmvhpkz.ppf
 
-C/ Interferometer configuration / noise estimate 
+F/ Interferometer configuration / noise estimate 
 =================================================
 Note : to be executed in spiapp 
 
-C.1/ formnoise.pic + vnpairs.cc + vnpairs.h 
+F.1/ formnoise.pic + vnpairs.cc + vnpairs.h 
 
-C.2/ pa.pic + pa.cc + pa.h : Lobe synthesis (dipoles sum)
+F.2/ pa.pic + pa.cc + pa.h : Lobe synthesis (dipoles sum)
 
-D/ 3D reconstruction (alpha,delta,freq) simulation program 
+G/ 3D reconstruction (alpha,delta,freq) simulation program 
 ===========================================================
 
Index: /trunk/Cosmo/RadioBeam/calcpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3783)
@@ -0,0 +1,96 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de calcul du spectre de puissance (3D) a partir d'un 
+  cube de donnees (delta rho/rho ou delta T/T)
+    R. Ansari , C. Magneville - Juin 2010 
+
+  Usage: calcpk In3DMap OutPk  [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 "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 OutPk  [InRenormFactor] " << endl;
+    return 1;
+  }
+  Timer tm("calcpk");
+  int rc = 0;
+  try {
+    string inppfname = arg[1];
+    string outname = arg[2];
+    r_4 renfact=1.;
+    bool fgrenfact=false;
+    if (narg>3) {
+      renfact=atof(arg[3]);
+      fgrenfact=true;
+    }
+    cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
+    TArray<r_4> inmap;
+    {
+      PInPersist pin(inppfname);
+      pin >> inmap;
+    }
+    if (fgrenfact) inmap *= renfact;
+    double mean, sigma;
+    MeanSigma(inmap, mean, sigma);
+    cout << " InMap sizes " << inmap.InfoString() << endl;
+    inmap.Show(); 
+    cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
+    tm.Split(" After read ");
+    
+    cout << "calcpk[2] : computing 3D Fourier coefficients ... " << endl;
+    FFTWServer ffts;
+    TArray< complex<r_4> > four3d;
+    ffts.FFTForward(inmap, four3d);
+    tm.Split(" After FFTForward ");
+
+    cout << "calcpk[3] : computing power spectrum ... " << endl;
+    RandomGenerator rg;
+    Four3DPk pkc(four3d, rg);
+    
+    HProf hp = pkc.ComputePk(0.,256);
+    {
+      cout << "calcpk[4] : writing profile P(k) to PPF file  " << outname << endl;
+      POutPersist po(outname);
+      po << hp;
+    }
+    tm.Split(" Done CalcP(k) ");    
+  }  // End of try bloc 
+  catch (PThrowable & exc) {  // catching SOPHYA exceptions
+    cerr << " calcpk.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 << " calcpk.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
+    rc = 98;
+  }
+  catch (...) {  // catching other exceptions
+    cerr << " calcpk.cc: some other exception (...) was caught ! " << endl;
+    rc = 97;
+  }
+  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
+  return rc;    
+}
Index: /trunk/Cosmo/RadioBeam/makefile
===================================================================
--- /trunk/Cosmo/RadioBeam/makefile	(revision 3782)
+++ /trunk/Cosmo/RadioBeam/makefile	(revision 3783)
@@ -6,21 +6,70 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise  
+all : pknoise calcpk calcpk2 syncube tjyk
 
 clean : 
 	rm Objs/* 
 
+### les executables 
 pknoise : Objs/pknoise
 	echo 'makefile : pknoise made'
 
-Objs/pknoise : Objs/pknoise.o Objs/specpk.o Objs/mdish.o
-	$(CXXLINK) -o Objs/pknoise Objs/pknoise.o Objs/specpk.o Objs/mdish.o $(SOPHYAEXTSLBLIST)
+calcpk : Objs/calcpk
+	echo 'makefile : calcpk made'
 
-Objs/pknoise.o : pknoise.cc specpk.h  mdish.h 
+calcpk2 : Objs/calcpk2
+	echo 'makefile : calcpk2 made'
+
+syncube : Objs/syncube
+	echo 'makefile : syncube made'
+
+tjyk : Objs/tjyk
+	echo 'makefile : tjyk made'
+
+### programme pknoise
+Objs/pknoise : Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
+	$(CXXLINK) -o Objs/pknoise Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+
+Objs/pknoise.o : pknoise.cc specpk.h  mdish.h qhist.h 
 	$(CXXCOMPILE) -o Objs/pknoise.o pknoise.cc
 
-Objs/specpk.o : specpk.cc specpk.h  mdish.h 
+### programme calcpk
+Objs/calcpk : Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
+	$(CXXLINK) -o Objs/calcpk Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+
+Objs/calcpk.o : calcpk.cc specpk.h  mdish.h qhist.h 
+	$(CXXCOMPILE) -o Objs/calcpk.o calcpk.cc
+
+### programme calcpk2
+Objs/calcpk2 : Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o 
+	$(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
+
+Objs/calcpk2.o : calcpk2.cc specpk.h  mdish.h qhist.h 
+	$(CXXCOMPILE) -o Objs/calcpk2.o calcpk2.cc
+
+### programme syncube
+Objs/syncube : Objs/syncube.o 
+	$(CXXLINK) -o Objs/syncube Objs/syncube.o $(SOPHYAEXTSLBLIST)
+
+Objs/syncube.o : syncube.cc 
+	$(CXXCOMPILE) -o Objs/syncube.o syncube.cc
+
+### programme syncube
+Objs/tjyk : Objs/tjyk.o Objs/radutil.o 
+	$(CXXLINK) -o Objs/tjyk Objs/tjyk.o Objs/radutil.o $(SOPHYAEXTSLBLIST)
+
+Objs/tjyk.o : tjyk.cc radutil.h 
+	$(CXXCOMPILE) -o Objs/tjyk.o tjyk.cc
+
+### les classes / fonctions 
+Objs/specpk.o : specpk.cc specpk.h  mdish.h qhist.h 
 	$(CXXCOMPILE) -o Objs/specpk.o specpk.cc
 
-Objs/mdish.o : mdish.cc mdish.h 
+Objs/mdish.o : mdish.cc mdish.h qhist.h 
 	$(CXXCOMPILE) -o Objs/mdish.o mdish.cc
+
+Objs/qhist.o : qhist.cc qhist.h 
+	$(CXXCOMPILE) -o Objs/qhist.o qhist.cc
+
+Objs/radutil.o : radutil.cc radutil.h 
+	$(CXXCOMPILE) -o Objs/radutil.o radutil.cc
Index: /trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3782)
+++ /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3783)
@@ -102,69 +102,4 @@
 };
 
-// -----------------------------------
-// -- Classe ressemblant a un histo 2D 
-// -----------------------------------
-QHis2D::QHis2D()
-  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
-{
-  ixb0 = jyb0 = 0;
-}
-QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
-  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
-{
-  Define(xMin, xMax, nxb, yMin, yMax, nyb);
-}
-void QHis2D::Define(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
-{
-  nx=nxb; ny=nyb;
-  xmin=xMin;  xmax=xMax;  
-  ymin=yMin;  ymax=yMax;  
-  dxb=(xmax-xmin)/(double)nx;
-  dyb=(ymax-ymin)/(double)ny;
-  sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
-  aw.ReSize(2,sz);
-  SetZeroBin();
-  sumw0=0.;
-  return;
-}
-double QHis2D::Add(r_8 x, r_8 y, r_8 w, bool fgfh)
-{
-  sa_size_t ix = (sa_size_t)((x-xmin)/dxb);
-  sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
-  if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
-  double rw = ((ix==ixb0)&&(jy==jyb0)) ? w : 0.;
-  sumw0 += rw; 
-  if (fgfh) aw(ix,jy) += w;
-  return rw;
-}
-void QHis2D::SetZeroBin(r_8 x, r_8 y)
-{
-  ixb0 = (sa_size_t)((x-xmin)/dxb);
-  jyb0 = (sa_size_t)((y-ymin)/dyb);
-}
-Histo2D QHis2D::Convert()
-{
-  int_4 imn,jmn,imx,jmx;
-  r_8 min = aw(0,0);
-  r_8 max = aw(0,0);
-  imn=jmn=imx=jmx=0;
-  Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
-  for(int_4 j=0; j<ny; j++) 
-    for(int_4 i=0; i<nx; i++) {
-      h2(i,j) = aw(i,j);
-      if (aw(i,j)>max) {
-	imx=i;  jmx=j;  max=aw(i,j);
-      }
-      if (aw(i,j)<min) {
-	imn=i;  jmn=j;  min=aw(i,j);
-      }
-    }
-  cout << "QHis2D::Convert()/Info: Nx,Ny=" << nx << "," << ny << " SumW=" << sumw0 
-       << "\n ... Max:" << imx << "," << jmx << " ->" << max 
-       << " @" << imx*dxb+xmin << "," << jmx*dyb+ymin
-       << "\n ...Min:" << imn << "," << jmn << " ->" << min 
-       << " @" << imn*dxb+xmin << "," << jmn*dyb+ymin << endl;
-  return h2;
-}
 
 //----------------------------------------------------------------------
Index: /trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.h	(revision 3782)
+++ /trunk/Cosmo/RadioBeam/mdish.h	(revision 3783)
@@ -14,7 +14,9 @@
 #include "genericfunc.h"   // SOPHYA .h  
 #include "array.h"         // SOPHYA .h  
-#include "histats.h"       // SOPHYA .h 
+#include "qhist.h" 
 
+#ifndef DeuxPI
 #define DeuxPI 2.*M_PI 
+#endif
 
 // -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (angles theta,phi)
@@ -78,24 +80,4 @@
 };
 
-// -----------------------------------
-// -- Classe ressemblant a un histo 2D 
-class QHis2D {
-public:
-  QHis2D();
-  QHis2D(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
-  void Define(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
-  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
-  void SetZeroBin(r_8 x=0, r_8 y=0);
-  inline double WBinX() { return dxb; }
-  inline double WBinY() { return dyb; }
-  Histo2D Convert();
-  
-  r_8 xmin,xmax,ymin,ymax;
-  r_8 dxb,dyb;
-  sa_size_t nx,ny;
-  sa_size_t ixb0, jyb0;
-  TArray<r_8> aw;
-  double sumw0;
-};
 
 // -------------------------------------------------------------------
Index: /trunk/Cosmo/RadioBeam/qhist.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/qhist.cc	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/qhist.cc	(revision 3783)
@@ -0,0 +1,139 @@
+//  Simple specific histo-like classes 
+//  R. Ansari - Avril-Mai 2010 
+
+#include "qhist.h"
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 2D 
+// -----------------------------------
+
+/* --Methode-- */
+QHis1D::QHis1D()
+  : nx_(0),xmin_(0.),xmax_(0.),
+    nhis_(0), nunder_(0), nover_(0),
+    swunder_(0.), swover_(0.)
+{
+}
+
+/* --Methode-- */
+QHis1D::QHis1D(QHis1D const& a)
+  : nx_(a.nx_),xmin_(a.xmin_),xmax_(a.xmax_), aw_(a.aw_), 
+    nhis_(a.nhis_), nunder_(a.nunder_), nover_(a.nover_),
+    swunder_(a.swunder_), swover_(a.swover_)
+{
+}
+
+/* --Methode-- */
+QHis1D::QHis1D(r_8 xMin,r_8 xMax,int_4 nxb)
+  : nx_(0),xmin_(0.),xmax_(0.),
+    nhis_(0), nunder_(0), nover_(0),
+    swunder_(0.), swover_(0.)
+{
+  Define(xMin, xMax, nxb);
+}
+
+/* --Methode-- */
+void QHis1D::Define(r_8 xMin,r_8 xMax,int_4 nxb)
+{
+  nx_=nxb; 
+  xmin_=xMin;  xmax_=xMax;  
+  dxb_=(xmax_-xmin_)/(double)nx_;
+  sa_size_t sz[5];  sz[0]=nx_; 
+  aw_.ReSize(1,sz);
+  nhis_=nunder_=nover_=0;
+  swunder_=swover_=0.;
+  return;
+}
+
+/* --Methode-- */
+sa_size_t QHis1D::Add(r_8 x, r_8 w)
+{
+  sa_size_t ix = (sa_size_t)((x-xmin_)/dxb_);
+  if (ix<0) {
+    nunder_++;  swunder_+=w;
+    return 0.;
+  }
+  if (ix>=nx_) {
+    nover_++;  swover_+=w;
+    return 0.;
+  }
+  nhis_++;
+  aw_(ix) += w;
+  return nhis_;
+}
+
+Histo QHis1D::Convert()
+{
+  Histo h(xmin_,xmax_,nx_);
+  for(sa_size_t i=0; i<nx_; i++)  h(i) = aw_(i);
+  cout << "QHis1D::Convert()/Info: Nx" << nx_ << " NUnder=" << nunder_ << " NHis=" << nhis_ 
+       << " NOver=" << nover_ << "\n" 
+       << " SumWUnder=" << swunder_ << " SumW=" << aw_.Sum() << " SumWOver=" << swover_ << endl;
+  return h;
+}
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 2D 
+// -----------------------------------
+QHis2D::QHis2D()
+  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
+{
+  ixb0 = jyb0 = 0;
+}
+QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
+  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
+{
+  Define(xMin, xMax, nxb, yMin, yMax, nyb);
+}
+void QHis2D::Define(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
+{
+  nx=nxb; ny=nyb;
+  xmin=xMin;  xmax=xMax;  
+  ymin=yMin;  ymax=yMax;  
+  dxb=(xmax-xmin)/(double)nx;
+  dyb=(ymax-ymin)/(double)ny;
+  sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
+  aw.ReSize(2,sz);
+  SetZeroBin();
+  sumw0=0.;
+  return;
+}
+double QHis2D::Add(r_8 x, r_8 y, r_8 w, bool fgfh)
+{
+  sa_size_t ix = (sa_size_t)((x-xmin)/dxb);
+  sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
+  if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
+  double rw = ((ix==ixb0)&&(jy==jyb0)) ? w : 0.;
+  sumw0 += rw; 
+  if (fgfh) aw(ix,jy) += w;
+  return rw;
+}
+void QHis2D::SetZeroBin(r_8 x, r_8 y)
+{
+  ixb0 = (sa_size_t)((x-xmin)/dxb);
+  jyb0 = (sa_size_t)((y-ymin)/dyb);
+}
+Histo2D QHis2D::Convert()
+{
+  int_4 imn,jmn,imx,jmx;
+  r_8 min = aw(0,0);
+  r_8 max = aw(0,0);
+  imn=jmn=imx=jmx=0;
+  Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
+  for(int_4 j=0; j<ny; j++) 
+    for(int_4 i=0; i<nx; i++) {
+      h2(i,j) = aw(i,j);
+      if (aw(i,j)>max) {
+	imx=i;  jmx=j;  max=aw(i,j);
+      }
+      if (aw(i,j)<min) {
+	imn=i;  jmn=j;  min=aw(i,j);
+      }
+    }
+  cout << "QHis2D::Convert()/Info: Nx,Ny=" << nx << "," << ny << " SumW=" << sumw0 
+       << "\n ... Max:" << imx << "," << jmx << " ->" << max 
+       << " @" << imx*dxb+xmin << "," << jmx*dyb+ymin
+       << "\n ...Min:" << imn << "," << jmn << " ->" << min 
+       << " @" << imn*dxb+xmin << "," << jmn*dyb+ymin << endl;
+  return h2;
+}
Index: /trunk/Cosmo/RadioBeam/qhist.h
===================================================================
--- /trunk/Cosmo/RadioBeam/qhist.h	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/qhist.h	(revision 3783)
@@ -0,0 +1,58 @@
+//  Simple specific histo-like classes 
+//  R. Ansari - Avril-Mai 2010 
+
+#ifndef QHIST_SEEN
+#define QHIST_SEEN
+
+#include "histats.h"       // SOPHYA .h 
+
+#include "machdefs.h"      // SOPHYA .h  
+#include "sopnamsp.h"      // SOPHYA .h  
+#include <iostream>
+
+#include "histats.h"       // SOPHYA .h 
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 1D 
+class QHis1D {
+  QHis1D();
+  QHis1D(QHis1D const& a);
+  QHis1D(r_8 xMin,r_8 xMax,int_4 nxBin);
+  void Define(r_8 xMin,r_8 xMax,int_4 nxBin); 
+  sa_size_t Add(r_8 x, r_8 w=1.);
+  inline double WBinX() { return dxb_; }
+  inline double NBinX() { return nx_; }
+
+  Histo Convert();
+  
+  r_8 xmin_,xmax_;
+  r_8 dxb_;
+  sa_size_t nx_;
+  TArray<r_8> aw_;
+  sa_size_t nhis_, nunder_, nover_;
+  r_8 swunder_, swover_ ; 
+};
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 2D 
+class QHis2D {
+public:
+  QHis2D();
+  QHis2D(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
+  void Define(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
+  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
+  void SetZeroBin(r_8 x=0, r_8 y=0);
+  inline double WBinX() { return dxb; }
+  inline double WBinY() { return dyb; }
+  Histo2D Convert();
+  
+  r_8 xmin,xmax,ymin,ymax;
+  r_8 dxb,dyb;
+  sa_size_t nx,ny;
+  sa_size_t ixb0, jyb0;
+  TArray<r_8> aw;
+  double sumw0;
+};
+
+
+#endif
Index: /trunk/Cosmo/RadioBeam/radutil.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/radutil.cc	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/radutil.cc	(revision 3783)
@@ -0,0 +1,29 @@
+#include "radutil.h"
+
+double H21Conversions::SpeedOfLight_Cst = 2.99792458e8;  // metres/sec
+double H21Conversions::Freq021cm_Cst = 1420.4;  // MHz 
+double H21Conversions::k_Boltzman_Cst = 1.3806503e-23;  // Boltzman constant SI (J K^-1)
+
+H21Conversions::H21Conversions(double freq, double opix)
+{
+  if (freq<1.e-6)  freq=Freq021cm_Cst;
+  setFrequency(freq);
+  setOmegaPix(opix);
+}
+
+void H21Conversions::setFrequency(double nu)
+{
+  freq_=nu;
+  z_=Freq021cm_Cst/freq_-1.;
+  lambda_=SpeedOfLight_Cst/freq_/1.e6;
+}
+
+double H21Conversions::toJansky(double temp)
+{
+  return 2.*k_Boltzman_Cst*temp/lambda_/lambda_ * omegapix_ * 1.e26;
+}
+
+double H21Conversions::toKelvin(double jy)
+{
+  return jy*1e-26 / omegapix_ * lambda_ * lambda_ / 2. / k_Boltzman_Cst ;
+}
Index: /trunk/Cosmo/RadioBeam/radutil.h
===================================================================
--- /trunk/Cosmo/RadioBeam/radutil.h	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/radutil.h	(revision 3783)
@@ -0,0 +1,54 @@
+//  Classes to compute simple quantities for radio / 21 cm 
+//  R. Ansari - June 2010 
+
+#ifndef RADUTIL_SEEN
+#define RADUTIL_SEEN
+
+#include "machdefs.h"      // SOPHYA .h  
+#include "sopnamsp.h"      // SOPHYA .h  
+
+#include "vector3d.h"      // conversion d'angle 
+
+class H21Conversions {
+public:
+  // redshift nul, angle solide du pixel en steradian 
+  H21Conversions(double freq=0., double opix=1.e-5); 
+
+  double toJansky(double temp);  // Conversion d'une temperature (Kelvin) en jansky 
+  double toKelvin(double jy);    // Conversion de Jansky en  temperature (Kelvin) 
+
+  void setFrequency(double nu);  // on definit la frequence en MHz 
+  inline void setRedshift(double z)    // on definit le redshift
+  { setFrequency(Freq021cm_Cst/(1+z)); }
+  inline void setLambda(double lam)    // on definit la longueur d'onde en m  
+  { setFrequency(SpeedOfLight_Cst/lam); }
+
+  inline void setOmegaPix(double opix)   // angle solide en steradian
+  { omegapix_ = opix; }
+  inline void setOmegaPixDeg2(double opix)  // angle solide en Deg^2
+  { double cf=Angle(1.,Angle::Degree).ToRadian(); omegapix_ = opix*cf*cf; }
+  inline void setOmegaPixArcmin2(double opix)  // angle solide en Arcmin^2
+  { double cf=Angle(1.,Angle::ArcMin).ToRadian(); omegapix_ = opix*cf*cf; }
+
+  inline double getRedshift() { return z_; }
+  inline double getFrequency() { return freq_; }
+  inline double getLambda() { return lambda_; } 
+
+  inline double getOmegaPix()   // angle solide en steradian
+  { return omegapix_; }
+  inline double getOmegaPixArcmin2()  // angle solide en Arcmin^2
+  { double cf=Angle(1.,Angle::ArcMin).ToRadian(); return omegapix_/cf/cf; }
+  inline double getOmegaPixDeg2()  // angle solide en Deg^2
+  { 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
+  static double k_Boltzman_Cst;         // Boltzmann constant  (SI Units)
+
+  double z_;
+  double freq_;
+  double lambda_;
+  double omegapix_;
+};
+
+#endif
Index: /trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/specpk.cc	(revision 3782)
+++ /trunk/Cosmo/RadioBeam/specpk.cc	(revision 3783)
@@ -225,9 +225,10 @@
   double kxx, kyy, kzz;
   // sa_size_t is large integer type  
-  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+  // We ignore 0th term in all frequency directions ...
+  for(sa_size_t kz=1; kz<fourAmp.SizeZ(); kz++) {
     kzz =  (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_; 
-    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+    for(sa_size_t ky=1; ky<fourAmp.SizeY(); ky++) {
       kyy =  (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
-      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+      for(sa_size_t kx=1; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
 	double kxx=(double)kx*dkx_;
 	complex<TF> za = fourAmp(kx, ky, kz);
Index: /trunk/Cosmo/RadioBeam/syncube.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3783)
@@ -0,0 +1,189 @@
+/*  ------------------------ 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"
+
+//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------   
+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("syncube");
+
+  if (narg < 3) {
+    cout << "Usage: syncube InFitsName Out3DPPFName [Out2DMapName] \n" << endl;
+    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 
+  string outname = arg[2];
+  string inname = arg[1];
+  int rc = 91;
+
+  cout << " ====== syncube :   Input map name= " << inname << " OutName=" << outname;
+  bool fginmap=true;
+  try {
+    SphereHEALPix<r_4> inmap0;
+    if (fginmap) {  // Lecture eventuelle de la carte en entree 
+      cout << " syncube[0]: Reading input map : " << inname << endl;
+      FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
+      fis >> inmap0;
+      cout << inmap0;
+    }
+    // On fait une carte en doublant la resolution 
+    int nside = inmap0.SizeIndex() ;  // Parametre de pixelisation HEALPix - definit la resolution
+    SphereHEALPix<r_4> inmap(2*nside);
+    for(int_4 kk=0; kk<inmap.NbPixels(); kk++) {
+      double theta, phi;  // Theta, Phi en radians
+      inmap.PixThetaPhi(kk, theta, phi);
+      inmap(kk) = inmap0(theta, phi);
+    }    
+    //    Sph2Sph(inmap0, inmap);
+    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.);
+    for(int_4 kk=0; kk<outmap.NbPixels(); kk++) {
+      double theta, phi;  // Theta, Phi en radians
+      outmap.PixThetaPhi(kk, theta, phi);
+      outmap(kk) = inmap(theta, phi);
+    }    
+
+    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;
+    for (sa_size_t j=0; j<omap.SizeY(); j++)  {
+      double theta = j*dtet+tet0;
+      for (sa_size_t i=0; i<omap.SizeX(); i++)  {
+	double phi = i*dphi+phi0;
+	omap(i,j) = inmap(theta, phi);
+      }
+    }
+    //    Sph2Sph(inmap,outmap);
+    
+    cout << "syncube[2]: Output rectangular map computed " << outmap << endl;
+    if (narg > 3) {
+      string ppfname = arg[3];
+      cout << " syncube[3]: Saving inmap/outmap tp PPF file-> " << ppfname << endl;
+      POutPersist po(ppfname);
+      po << PPFNameTag("inmap") << inmap;
+      po << PPFNameTag("outmap") << outmap;
+      po << PPFNameTag("omap") << omap;
+    }
+
+    sa_size_t NFreq = 128;
+    TArray<r_4> ocube(NPhi,NTet,NFreq);
+
+    double infreq = 400.; // Frequence carte input en MHz
+    double freq0 = 840.;  // Freq0 du cube de sortie 
+    double dfreq = 1.;  
+
+    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(sigPLidx1,PLidx1);
+	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 ajoute une autre composante avec un indice spectral different
+	freqexpo = rg.Gaussian(sigPLidx2,PLidx2); 
+	double famp = rg.Flat01()*AmpPL2;
+	for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
+	  double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
+	  ocube(i,j,k) += famp*omap(i,j)*rapfreq;
+	}
+
+      }
+    }
+    
+    // On sauve le cube de sortie
+    {
+      cout << " syncube[2]: Saving output cube to -> " << outname << endl;
+      POutPersist poc(outname);
+      poc << ocube;
+    }
+
+    rc = 0;
+  }
+  catch (PThrowable& exc) {
+    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n syncube.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+  }
+  catch (...) {
+    cerr << " syncube.cc catched unknown (...) exception  " << endl; 
+    rc = 78; 
+  } 
+
+  cout << ">>>> syncube[9] ------- FIN ----------- Rc=" << rc << endl;
+  return rc;
+}
+
+
Index: /trunk/Cosmo/RadioBeam/tjyk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3783)
+++ /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3783)
@@ -0,0 +1,18 @@
+#include "radutil.h"
+#include <iostream>
+
+int main(int narg, char* arg[])
+{
+
+  H21Conversions conv;
+  
+  conv.setRedshift(0.6);
+  conv.setOmegaPixArcmin2(10.*10.);
+  cout << " --------------------------------------------------------------- " << endl;
+  cout << " H21Conversions: z=" << conv.getRedshift() << " Freq=" << conv.getFrequency()
+       << " 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 << " --------------------------------------------------------------- " << endl;
+
+  return 0;
+}
