Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3787)
@@ -16,6 +16,8 @@
   - plpkn.pic : spiapp script for plotting pknoise.cc program output 
 A.2/ power spectrum computation and foreground/radio source 
+  - subtractradsrc.cmd : Explanation and commands to compute various power spectra and foreground subtraction 
   - 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
+  - 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
@@ -24,4 +26,5 @@
 A.3/ List of other common files :
   - specpk.h specpk.cc : 3D power spectrum computation classes  
+  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
   - radutil.h radutil.cc : utilitaire de conversion a 21 cm
Index: /trunk/Cosmo/RadioBeam/applobe.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/applobe.cc	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/applobe.cc	(revision 3787)
@@ -0,0 +1,132 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de convolution avec le lobe d'un cube 3D (angles,freq) 
+    R. Ansari , C. Magneville - Juin 2010 
+
+  Usage: applobe In3DPPFName Out3DPPFName
+---------------------------------------------------------------  */
+
+#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 "vector3d.h"
+
+// #include "xastropack.h"     -- Pour faire les conversions de coordonnees celestes
+
+#include "radutil.h"
+#include "lobe.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("applobe");
+
+  if (narg < 3) {
+    cout << "Usage: applobe In3DPPFName Out3DPPFName \n" << endl;
+    return 1;
+  }
+
+  // decodage arguments 
+  string outname = arg[2];
+  string inname = arg[1];
+  int rc = 91;
+
+  cout << " ====== applobe :   Input NVSS catalog name= " << inname << " OutName=" << outname;
+  bool fginmap=true;
+  try {
+    TArray<r_4> incube;
+    cout << "applobe[1]: reading input 3D map (cube) from file " << inname << endl;
+    {
+      PInPersist pin(inname);
+      pin >> incube;
+    }
+    cout << incube;
+
+    double dxdeg = ThetaSizeDegre/(double)NTheta;
+    double dydeg = PhiSizeDegre/(double)NPhi;
+    double dx = DegreeToRadian(dxdeg);
+    double dy = DegreeToRadian(dydeg);
+    double dfreq = FreqSizeMHz/(double)NFreq;
+
+    cout << " X,Y map size in degrees , X/Phi=" << PhiSizeDegre << " Y/Theta=" << ThetaSizeDegre
+	 << " \n dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx << " dy_rad=" << dy << ")" 
+         << " FreqSize=" << FreqSizeMHz << " dfreq=dz= " << dfreq << " MHz" << endl;
+
+    double mean, sigma;
+    MeanSigma(incube, mean, sigma);
+    cout << " InCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;
+
+    cout << "applobe[2]: creating Four2DResponse and BeamEffect objects..." << endl;
+    // Dish de 64 m de diametre
+    double Diametre = 64.;
+    H21Conversions conv;
+    conv.setFrequency(Freq0MHz);
+    double lambda = conv.getLambda();
+    Four2DResponse fresp(2, Diametre/lambda, Diametre/lambda);
+    BeamEffect beam(fresp);
+
+    cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
+    beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
+
+    cout << "applobe[4]: calling ReSample(incube, 0.25, 0.25, 1.) ... " << endl;
+    TArray< r_4 > outcube = beam.ReSample(incube, 0.25, 0.25, 1.);
+    
+    cout << outcube;
+    MeanSigma(outcube, mean, sigma);
+    cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;
+
+    // On sauve le cube de sortie
+    {
+      cout << " applobe[5]: Saving output cube to -> " << outname << endl;
+      POutPersist poc(outname);
+      poc << outcube;
+    }
+
+    rc = 0;
+  }
+  catch (PThrowable& exc) {
+    cerr << " applobe.cc catched SOPHYA Exception " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n applobe.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+  }
+  catch (...) {
+    cerr << " applobe.cc catched unknown (...) exception  " << endl; 
+    rc = 78; 
+  } 
+
+  cout << ">>>> applobe[9] ------- FIN ----------- Rc=" << rc << endl;
+  return rc;
+}
+
+
Index: /trunk/Cosmo/RadioBeam/calcpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3787)
@@ -19,4 +19,6 @@
 
 #include "qhist.h" 
+#include "lobe.h" 
+
 #include "histinit.h"
 #include "fftwserver.h"
@@ -33,5 +35,5 @@
 {
   if (narg<3) {
-    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] " << endl;
+    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
     return 1;
   }
@@ -41,5 +43,5 @@
     string inppfname = arg[1];
     string outname = arg[2];
-    r_4 renfact=1.;
+    double renfact=1.;
     bool fgrenfact=false;
     if (narg>3) {
@@ -47,4 +49,11 @@
       fgrenfact=true;
     }
+    double pixsignoise = 0.;
+    bool fgaddnoise=false;
+    if (narg>4) {
+      pixsignoise=atof(arg[3]);
+      fgaddnoise=true;
+    }
+
     cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
     TArray<r_4> inmap;
@@ -53,5 +62,8 @@
       pin >> inmap;
     }
-    if (fgrenfact) inmap *= renfact;
+    if (fgrenfact) { 
+      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
+      inmap *= renfact;
+    }
     double mean, sigma;
     MeanSigma(inmap, mean, sigma);
@@ -59,4 +71,8 @@
     inmap.Show(); 
     cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
+    if (fgaddnoise) { 
+      BeamEffect::AddNoise(inmap, pixsignoise);
+    }
+
     tm.Split(" After read ");
     
Index: /trunk/Cosmo/RadioBeam/calcpk2.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3787)
@@ -7,5 +7,5 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk 
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
 ---------------------------------------------------------------  */
 
@@ -22,4 +22,7 @@
 
 #include "qhist.h" 
+#include "lobe.h" 
+#include "cubedef.h" 
+
 #include "histinit.h"
 #include "fftwserver.h"
@@ -30,5 +33,6 @@
 typedef ThSDR48RandGen RandomGenerator ;
 //--- Declaration des fonctions 
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx);
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
+                            TArray<r_8>& synctemp, TArray<r_8>& specidx);
 
 //-------------------------------------------------------------------------
@@ -39,5 +43,5 @@
 {
   if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk  " << endl;
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk [PixNoiseLevel] " << endl;
     return 1;
   }
@@ -50,4 +54,13 @@
     r_4 rfacsync = atof(arg[4]);
     string outname = arg[5];
+
+    double pixsignoise = 0.;
+    bool fgaddnoise=false;
+    if (narg>6) {
+      pixsignoise=atof(arg[6]);
+      fgaddnoise=true;
+    }
+
+
     TArray<r_4> inmaplss, inmapsync;
     const char * tits[2]={"LSS", "Sync"};
@@ -68,8 +81,16 @@
       cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
     }
+
+    if (fgaddnoise) { 
+      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
+      BeamEffect::AddNoise(inmaplss, pixsignoise);
+    }
+
     tm.Split(" After read ");
     TArray<r_8> synctemp, specidx; 
     cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
-    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, synctemp, specidx);
+    double freq0 = Freq0MHz;
+    double dfreq = FreqSizeMHz/(double)NFreq;
+    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, freq0, dfreq, synctemp, specidx);
     double mean, sigma;
     MeanSigma(inmap, mean, sigma);
@@ -117,5 +138,9 @@
 
 /* --Fonction-- */
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx)
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
+			    TArray<r_8>& synctemp, TArray<r_8>& specidx)
+// Inputs : maplss, mapsyc, freq0, dfreq
+// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
+// Return_Array : foreground subtracted LSS signal 
 {
   bool smo;
@@ -131,6 +156,6 @@
   Vector vlnf(maplss.SizeZ());
   int nprt = 0;
-  double freq0=840.;  // Frequence premier index en k (MHz)
-  double dfreq=1.;   // largeur en frequence de chaque plan (Mhz)  
+  // double freq0 : Frequence premier index en k (MHz)
+  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
   for(sa_size_t i=0; i<maplss.SizeX(); i++) 
     for(sa_size_t j=0; j<maplss.SizeY(); j++)  {
Index: /trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3787)
@@ -6,15 +6,16 @@
 
 // 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 sa_size_t NTheta=360;
+static sa_size_t NPhi=360;
+static sa_size_t NFreq = 200;
+
+static double Theta0Degre = 75.;  // -> Delta = +30 deg
+static double Phi0Degre = 135.;   // -> alpha = 9h00
+static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
+static double PhiSizeDegre = 30.;    // 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)
+static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
 
 
@@ -31,3 +32,16 @@
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
+/*
+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 = 100.;    // Taille de la carte en MHz (axe Z)
+*/
+
 #endif
Index: /trunk/Cosmo/RadioBeam/lobe.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/lobe.cc	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/lobe.cc	(revision 3787)
@@ -0,0 +1,103 @@
+
+#include "lobe.h"
+#include "radutil.h"
+#include "randfmt.h"      
+typedef FMTRandGen RandomGenerator ;
+
+
+#include "fftwserver.h"    
+#include "ctimer.h"    
+
+
+/* --Methode-- */
+BeamEffect::BeamEffect(Four2DResponse& resp)
+  : fresp_(resp)
+// resp doit avoir sa longueur d'onde de reference en metres
+{
+}
+
+/* --Methode-- */
+void BeamEffect::ApplyLobeK2D(TArray< complex<TF> >& fourAmp, double dkx, double dky, double lambda)
+//  dx, dy en radians, lambda en metres 
+{
+  fresp_.setLambda(lambda);
+  double kxx, kyy;
+  for(sa_size_t ky=0; 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++) {
+      kxx=(double)kx*dkx;
+      fourAmp(kx, ky) *= complex<TF>(fresp_(kxx, kyy), 0.);   	
+    }
+  }
+  return;
+}
+
+
+/* --Methode-- */
+void BeamEffect::ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df)
+// dx, dy en radioans, f0, df en MHz
+{
+  Timer tm("BeamEffect::ApplyLobe3D");
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  
+  H21Conversions conv;
+  
+  TArray< complex<TF> > fourAmp;
+  double dkx = DeuxPI/(double)a.SizeX()/dx;
+  double dky = DeuxPI/(double)a.SizeY()/dy;
+  
+  for(sa_size_t kz=0; kz<a.SizeZ(); kz++) {
+    TArray< TF > slice( a(Range::all(), Range::all(), kz) );
+    ffts.FFTForward(slice, fourAmp);
+    conv.setFrequency(f0+kz*df);
+    ApplyLobeK2D(fourAmp, dkx, dky, conv.getLambda());
+    ffts.FFTBackward(fourAmp, slice, true);
+    if (kz%10==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  return;
+}
+
+/* --Methode-- */
+TArray< TF > BeamEffect::ReSample(TArray< TF >& a, double xfac, double yfac, double zfac)
+{
+  Timer tm("BeamEffect::ReSample");
+
+  sa_size_t szx = a.SizeX()*xfac;
+  sa_size_t szy = a.SizeY()*yfac;
+  sa_size_t szz = a.SizeZ()*zfac;
+
+  TArray<TF> rsa(szx, szy, szz);
+  for(sa_size_t kz=0; kz<rsa.SizeZ(); kz++) {
+    sa_size_t kza=kz/zfac;
+    if ((kza<0)||(kza>=a.SizeZ()))  continue;
+    for(sa_size_t ky=0; ky<rsa.SizeY(); ky++) {
+      sa_size_t kya=ky/yfac;
+      if ((kya<0)||(kya>=a.SizeY()))  continue;
+      for(sa_size_t kx=0; kx<rsa.SizeX(); kx++) {
+	sa_size_t kxa=kx/xfac;
+	if ((kxa<0)||(kxa>=a.SizeX()))  continue;
+	rsa(kx,ky,kz)=a(kxa,kya,kza);
+      }
+    }
+  }
+  return rsa;
+}
+
+
+/* --Methode-- */
+void BeamEffect::AddNoise(TArray< TF >& a, double pixsignoise, bool fgcmsig)
+{
+  cout << "BeamEffect::AddNoise() PixelSigmaNoise=" << pixsignoise << endl;
+  RandomGenerator rg;
+  for(sa_size_t kz=0; kz<a.SizeZ(); kz++) 
+    for(sa_size_t ky=0; ky<a.SizeY(); ky++) 
+      for(sa_size_t kx=0; kx<a.SizeX(); kx++) 
+	a(kx,ky,kz) += rg.Gaussian(pixsignoise);
+  if (fgcmsig) {
+    double mean, sigma;
+    MeanSigma(a, mean, sigma);
+    cout << "BeamEffect::AddNoise()-done,  Mean=" << mean << " Sigma=" << sigma << endl;
+  }
+  return;
+}
Index: /trunk/Cosmo/RadioBeam/lobe.h
===================================================================
--- /trunk/Cosmo/RadioBeam/lobe.h	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/lobe.h	(revision 3787)
@@ -0,0 +1,37 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+    Calcul de l'effet de lobe sur Carte2D (angleX,angleY) et 
+    cube 3D (angleX,angleY)
+    R. Ansari , C. Magneville - Juin 2010 
+---------------------------------------------------------------  */
+
+#ifndef BEAMH_SEEN
+#define BEAMH_SEEN
+
+#include "machdefs.h"      // SOPHYA .h  
+#include "sopnamsp.h"      // SOPHYA .h  
+#include <math.h>
+#include <iostream>
+#include <vector>
+#include <string> 
+
+#include "tarray.h"
+#include "vector3d.h"
+#include "mdish.h"
+
+#ifndef TF 
+#define TF r_4 
+#endif
+
+class BeamEffect {
+public:
+  BeamEffect(Four2DResponse& resp);
+  void ApplyLobeK2D(TArray< complex<TF> >& f2d, double dkx, double dky, double lambda=1.);
+  void ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df);  
+
+  static TArray< TF > ReSample(TArray< TF >& a, double xfac, double yfac, double zfac);
+  static void AddNoise(TArray< TF >& a, double pixsignoise, bool fgcmsig=true);
+
+  Four2DResponse& fresp_;
+};
+
+#endif
Index: /trunk/Cosmo/RadioBeam/makefile
===================================================================
--- /trunk/Cosmo/RadioBeam/makefile	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/makefile	(revision 3787)
@@ -6,8 +6,11 @@
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : pknoise calcpk calcpk2 syncube srcat2cube tjyk
+all : pknoise calcpk calcpk2 syncube srcat2cube tjyk applobe
 
 clean : 
 	rm Objs/* 
+
+PKGOLIST =  Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o 
+PKGHLIST =  lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h 
 
 ### les executables 
@@ -30,28 +33,31 @@
 	echo 'makefile : tjyk made'
 
+applobe : Objs/applobe
+	echo 'makefile : applobe 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 : Objs/pknoise.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/pknoise Objs/pknoise.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
 
-Objs/pknoise.o : pknoise.cc specpk.h  mdish.h qhist.h 
+Objs/pknoise.o : pknoise.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/pknoise.o pknoise.cc
 
 ### 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 : Objs/calcpk.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/calcpk Objs/calcpk.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
 
-Objs/calcpk.o : calcpk.cc specpk.h  mdish.h qhist.h 
+Objs/calcpk.o : calcpk.cc $(PKGHLIST)
 	$(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 : Objs/calcpk2.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/calcpk2.o : calcpk2.cc specpk.h  mdish.h qhist.h 
+Objs/calcpk2.o : calcpk2.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/calcpk2.o calcpk2.cc
 
 ### programme syncube
-Objs/syncube : Objs/syncube.o 
-	$(CXXLINK) -o Objs/syncube Objs/syncube.o $(SOPHYAEXTSLBLIST)
+Objs/syncube : Objs/syncube.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/syncube Objs/syncube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
 Objs/syncube.o : syncube.cc 
@@ -59,16 +65,24 @@
 
 ### programme srcat2cube
-Objs/srcat2cube : Objs/srcat2cube.o 
-	$(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o  Objs/radutil.o $(SOPHYAEXTSLBLIST)
+Objs/srcat2cube : Objs/srcat2cube.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/srcat2cube.o : srcat2cube.cc radutil.h 
+Objs/srcat2cube.o : srcat2cube.cc $(PKGHLIST)
 	$(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)
+Objs/tjyk : Objs/tjyk.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/tjyk Objs/tjyk.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/tjyk.o : tjyk.cc radutil.h 
+Objs/tjyk.o : tjyk.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/tjyk.o tjyk.cc
+
+### programme applobe
+Objs/applobe : Objs/applobe.o $(PKGOLIST)
+	$(CXXLINK) -o Objs/applobe Objs/applobe.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
+
+Objs/applobe.o : applobe.cc $(PKGHLIST)
+	$(CXXCOMPILE) -o Objs/applobe.o applobe.cc
+
 
 ### les classes / fonctions 
@@ -79,4 +93,7 @@
 	$(CXXCOMPILE) -o Objs/mdish.o mdish.cc
 
+Objs/lobe.o : lobe.cc lobe.h mdish.h qhist.h 
+	$(CXXCOMPILE) -o Objs/lobe.o lobe.cc
+
 Objs/qhist.o : qhist.cc qhist.h 
 	$(CXXCOMPILE) -o Objs/qhist.o qhist.cc
Index: /trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3787)
@@ -12,4 +12,6 @@
   : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
 {
+  setLambdaRef();
+  setLambda();
 }
 
@@ -17,4 +19,6 @@
 double Four2DResponse::Value(double kx, double ky)
 {
+  kx *= lambda_ratio_;
+  ky *= lambda_ratio_;
   double wk,wkx,wky;
   switch (typ_) 
@@ -62,4 +66,6 @@
 double Four2DRespTable::Value(double kx, double ky)
 {
+  kx *= lambda_ratio_;
+  ky *= lambda_ratio_;
   int_4 i,j;
   if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) || 
Index: /trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.h	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/mdish.h	(revision 3787)
@@ -34,4 +34,9 @@
   { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }
 
+  inline void setLambdaRef(double lambda=1.)
+  { lambdaref_ = lambda; }
+  inline void setLambda(double lambda=1.)
+  { lambda_ = lambda;  lambda_ratio_ = lambdaref_/lambda_; }
+  
   // Return the 2D response for wave vector (kx,ky)
   virtual double Value(double kx, double ky); 
@@ -45,4 +50,7 @@
   int typ_;
   double dx_, dy_; 
+  double lambdaref_, lambda_;
+  double lambda_ratio_;     // lambdaref_/lambda_;
+
 };
 
Index: /trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/specpk.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/specpk.cc	(revision 3787)
@@ -151,5 +151,5 @@
   // fourAmp represent 3-D fourier transform of a real input array. 
   // The second half of the array along Y and Z contain negative frequencies
-  double kxx, kyy, kzz;
+  double kxx, kyy, kzz, rep, amp;
   // sa_size_t is large integer type  
   for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
@@ -158,10 +158,10 @@
       kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
       for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
-	double kxx=(double)kx*dkx_;
-	double rep = resp(kxx, kyy);
+	kxx=(double)kx*dkx_;
+	rep = resp(kxx, kyy);
 	if (crmask&&(kz==0))  mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
 	if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
 	else {
-	  double amp = 1./sqrt(rep)/sqrt(2.);
+	  amp = 1./sqrt(rep)/sqrt(2.);
 	  fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
 	}
Index: /trunk/Cosmo/RadioBeam/subtractradsrc.cmd
===================================================================
--- /trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3787)
@@ -0,0 +1,81 @@
+#####################################################################################
+####  Commands to run the different programs to produce foreground maps 
+####  and compute radio-source subtracted P(k)
+#####################################################################################
+
+### Cube definition in file cubedef.h 
+
+### Step 1/ Produce an LSS data cube  with appropriate size and redshift using SimLSS
+# 1.a/ Run SimLSS
+csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,4 -y 360,4 -z 128,5 -Z 0.6 -8 1. -n 10000 -O 0,2 -o map3dz06B -T 2
+# 1.b/ 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) 
+#  Execucte the following script in spiapp :
+
+csh> cat > racube.pic
+set f map3dz06B
+readfits ${f}_r.fits
+rename ${f}_r map
+print map 
+c++exec \
+  TArray<r_4> omap(map.SizeY(),map.SizeZ(),map.SizeX()-2 ); \
+  for(sa_size_t i=0;i<omap.SizeX();i++) \
+    for(sa_size_t j=0;j<omap.SizeY();j++) \
+      for(sa_size_t k=0;k<omap.SizeZ();k++) \
+        omap(i,j,k)=map(k+1,i,j);  \
+  KeepObj(omap); 
+
+print omap
+saveppf omap lsscube.ppf 
+
+csh> spiapp -term -exec lsscube.ppf
+
+## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
+# 2.a/ Synchrotron map from HASLAM 400 MHz map
+csh> ./Objs/syncube syncmap_eq.fits syncube.ppf
+# 2.b/ radio source cube from NVSS catalog
+csh> ./Objs/srcat2cube nvss.fits nvsscube.ppf
+# 2.c/ Add the two cubes using the following spiapp script 
+csh> cat > sumcubes.pic
+openppf syncube.ppf
+openppf nvsscube.ppf
+expmeansig syncube val 
+expmeansig nvsscube val
+c++exec TArray<r_4> fgndcube = syncube+nvsscube; KeepObj(fgndcube);
+print fgndcube
+expmeansig fgndcube val
+saveppf fgndcube fgndcube.ppf
+
+csh> spiapp -term -exec sumcubes.pic 
+
+## Step 3/ Apply lobe effect on foreground cube and LSS cube
+csh> ./Objs/applobe fgndcube.ppf fgndcube_lobe.ppf
+csh> ./Objs/applobe lsscube.ppf lsscube_lobe.ppf
+
+### Step 4/ Compute power spectra 
+## mass to temperature converion factor   CT21 ~= 0.2 mK 
+## Foreground maps are in temperature 
+## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
+## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
+## sigma_noise ~ 0.25 mK 
+# 4.a/ LSS power spectrum  without noise
+csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.2 
+# and with noise  
+csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.25 
+#  with the lobe effect 
+csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2  
+# 4.b/ Foreground power spectrum 
+csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
+csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
+
+# 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.25
+# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect 
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.25
+
+### Step 5 / Check the results using spiapp 
+openppf fgndpk.ppf
+openppf lsspk.ppf
+openppf lsspkwn.ppf
+openppf lsspklobe.ppf
Index: /trunk/Cosmo/RadioBeam/syncube.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/syncube.cc	(revision 3787)
@@ -115,5 +115,5 @@
     double mean, sigma;
     MeanSigma(omap, mean, sigma);
-    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
+    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << " Sizes:" << endl;
     omap.Show(); 
 
Index: /trunk/Cosmo/RadioBeam/tjyk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3786)
+++ /trunk/Cosmo/RadioBeam/tjyk.cc	(revision 3787)
@@ -6,4 +6,8 @@
 int main(int narg, char* arg[])
 {
+  if (narg<2) {
+    cout << " Usage: tjyk redshift [resol_arcmin] " << endl;
+    return 1;
+  }
 
   H21Conversions conv;
