Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3788)
@@ -25,6 +25,7 @@
 
 A.3/ List of other common files :
+  - fgndsub.h fgndsub.h : radio source, synchrotron subtraction class 
+  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - 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 3788)
@@ -68,5 +68,5 @@
       pin >> incube;
     }
-    cout << incube;
+    incube.Show(); 
 
     double dxdeg = ThetaSizeDegre/(double)NTheta;
@@ -85,10 +85,8 @@
 
     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);
+    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
     BeamEffect beam(fresp);
 
@@ -96,8 +94,8 @@
     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 << "applobe[4]: calling ReSample(incube, 0.5, 0.5, 1.) ... " << endl;
+    TArray< r_4 > outcube = beam.ReSample(incube, 0.5, 0.5, 1.);
     
-    cout << outcube;
+    outcube.Show(); 
     MeanSigma(outcube, mean, sigma);
     cout << " OutCube 3D- : Mean=" << mean << " Sigma=" << sigma <<  endl;
Index: /trunk/Cosmo/RadioBeam/calcpk2.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3788)
@@ -7,5 +7,6 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk 
+               [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]
 ---------------------------------------------------------------  */
 
@@ -20,8 +21,11 @@
 #include "specpk.h"
 #include "histats.h"
+#include "vector3d.h"
 
 #include "qhist.h" 
 #include "lobe.h" 
 #include "cubedef.h" 
+#include "fgndsub.h" 
+#include "radutil.h" 
 
 #include "histinit.h"
@@ -32,7 +36,4 @@
 
 typedef ThSDR48RandGen RandomGenerator ;
-//--- Declaration des fonctions 
-TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq, 
-                            TArray<r_8>& synctemp, TArray<r_8>& specidx);
 
 //-------------------------------------------------------------------------
@@ -43,5 +44,6 @@
 {
   if (narg<6) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk [PixNoiseLevel] " << endl;
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk \n" 
+         << "        [PixNoiseLevel] [TargetBeamArcmin] [NSigSrcThr]" << endl;
     return 1;
   }
@@ -61,13 +63,25 @@
       fgaddnoise=true;
     }
+    
+    bool fgcorrbeam=true;
+    double tbeamarcmin=30.;
+    if (narg>7) {
+      tbeamarcmin=atof(arg[7]);
+      if (tbeamarcmin<1.e-6)  fgcorrbeam=false;
+    }
+    bool fgclnsrc=true;
+    double nsigsrc=5.;
+    if (narg>8) {
+      nsigsrc=atof(arg[8]);
+      if (nsigsrc<1.e-6)  fgclnsrc=false;
+    }
 
-
-    TArray<r_4> inmaplss, inmapsync;
-    const char * tits[2]={"LSS", "Sync"};
+    TArray<r_4> maplss, mapsync;
+    const char * tits[2]={"LSS", "Sync/RadioSrc"};
     for(int ks=0; ks<2; ks++)  {
       string& ppfname=inppflss;
       r_4 rfac=rfaclss;
-      TArray<r_4>* inmap=&inmaplss; 
-      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&inmapsync; }
+      TArray<r_4>* inmap=&maplss; 
+      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&mapsync; }
       cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname 
 	   << "  RenormFactor=" << rfac << endl;
@@ -82,27 +96,56 @@
     }
 
+    bool smo; 
+    if (!maplss.CompareSizes(mapsync,smo) ) {
+      cout << " calcpk2/ERROR sizes " << endl;
+      maplss.Show();  mapsync.Show(); 
+      return 99;
+    }
+    
+    TArray<r_4> skycube(mapsync);
+    skycube += maplss;
+
     if (fgaddnoise) { 
-      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
-      BeamEffect::AddNoise(inmaplss, pixsignoise);
+      cout << " calcpk2: adding noise to skycube cube ... " << endl;
+      BeamEffect::AddNoise(skycube, pixsignoise);
     }
 
-    tm.Split(" After read ");
-    TArray<r_8> synctemp, specidx; 
-    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
-    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);
-    cout << " After cleaning: Mean=" << mean << " Sigma=" << sigma << endl;      
+    MeanSigma(skycube, mean, sigma);
+    cout << " input sky cube : Mean=" << mean << " Sigma=" << sigma << endl;      
+    tm.Split(" After input ");
+
+    H21Conversions conv;
+    conv.setFrequency(Freq0MHz);
+    double lambda = conv.getLambda();
+    Four2DResponse arep(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
+    double DoL = 1.22/ArcminToRadian(tbeamarcmin);
+    Four2DResponse tbeam(2, DoL, DoL );
+    ForegroundCleaner  cleaner(arep, tbeam, skycube); 
+    if (fgcorrbeam) {
+      cout << "calcpk2[3.a] : calling cleaner.BeamCorrections() for target beam (" << tbeamarcmin << " arcmin)" 
+	   << " Diam/Lambda=" << DoL << endl;
+      cleaner.BeamCorrections();
+    }
+    if (fgclnsrc) {
+      cout << "calcpk2[3.b] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
+      cleaner.CleanPointSources(nsigsrc);
+    }
+
+    cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
+    TArray<r_4> synctemp, specidx; 
+    TArray<r_4> exlss = cleaner.extractLSSCube(synctemp, specidx);
+
+    MeanSigma(exlss, mean, sigma);
+    cout << " After cleaning/extractLSS: Mean=" << mean << " Sigma=" << sigma << endl;      
     tm.Split(" After CleanForeground");
 
-    cout << "calcpk2[4] : computing 3D Fourier coefficients ... " << endl;
+    cout << "calcpk2[5] : computing 3D Fourier coefficients ... " << endl;
     FFTWServer ffts;
     TArray< complex<r_4> > four3d;
-    ffts.FFTForward(inmap, four3d);
+    ffts.FFTForward(exlss, four3d);
     tm.Split(" After FFTForward ");
 
-    cout << "calcpk2[5] : computing power spectrum ... " << endl;
+    cout << "calcpk2[6] : computing power spectrum ... " << endl;
     RandomGenerator rg;
     Four3DPk pkc(four3d, rg);
@@ -112,9 +155,12 @@
     tm.Split(" Done ComputePk ");    
     
-    cout << "calcpk2[4] : writing profile P(k) and foreground maps to  " << outname << endl;
+    cout << "calcpk2[7.a] : writing profile P(k)  to  " << outname << endl;
     POutPersist po(outname);
-    po << PPFNameTag("Pk") << hp;
-    po << PPFNameTag("Tsync") << synctemp;
-    po << PPFNameTag("async") << specidx;
+    po << hp;
+    outname = "fgm_" + outname;
+    cout << "calcpk2[7.b] : writing foreground maps to  " << outname << endl;
+    POutPersist pom(outname);
+    pom << PPFNameTag("Tsync") << synctemp;
+    pom << PPFNameTag("async") << specidx;
 
   }  // End of try bloc 
@@ -132,56 +178,7 @@
     rc = 97;
   }
-  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
+  cout << " ==== End of calcpk2.cc program  Rc= " << rc << endl;
   return rc;    
 }
 
 
-/* --Fonction-- */
-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;
-  if (!maplss.CompareSizes(mapsync,smo) ) {
-    cout << " CleanForeground/ERROR sizes " << endl;
-    maplss.Show();  mapsync.Show(); 
-    throw ParmError("CleanForeground- maplss , mapsync not the same size !");
-  }
-  sa_size_t sz[5];   sz[0]=maplss.SizeX();  sz[1]=maplss.SizeY();
-  synctemp.SetSize(2, sz); 
-  specidx.SetSize(2, sz);
-  TArray<r_4>& omap=maplss;
-  Vector vlnf(maplss.SizeZ());
-  int nprt = 0;
-  // 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++)  {
-      r_8 s1, sx, sx2, sy, sxy;
-      s1=sx=sx2=sy=sxy=0.;
-      for(sa_size_t k=0; k<maplss.SizeZ(); k++)  {
-	double lnf=log((double)k*dfreq+freq0);
-	vlnf(k)=lnf;
-	double ttot=(r_8)(mapsync(i,j,k)+maplss(i,j,k));
-	if (ttot < 1.e-5) continue;
-	double lntt=log(ttot);
-	s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
-	sy+=lntt;    sxy+=(lnf*lntt);
-      }
-      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
-      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
-      double T0 = exp(beta+alpha*vlnf(0));
-      if ((i%6==0)&&(j%7==0)) 
-	cout << "CleanForeground[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
-	     << " (mapsync=" << mapsync(i,j,0) << " ... " << mapsync(i,j,125) << ")" << endl; 
-      synctemp(i,j) = T0;
-      specidx(i,j) = alpha;
-      for(sa_size_t k=0; k<maplss.SizeZ(); k++) {
-	r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
-	omap(i,j,k) = mapsync(i,j,k)+maplss(i,j,k)-fittedtemp;
-      }
-    }
-  return omap;
-}
Index: /trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/cubedef.h	(revision 3788)
@@ -32,4 +32,7 @@
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
+// Parametres pour la reponse en Fourier de l'instrument :
+static double InterfArrayDiametre = 64.;
+
 /*
 static sa_size_t NTheta=256;
Index: /trunk/Cosmo/RadioBeam/fgndsub.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/fgndsub.cc	(revision 3788)
+++ /trunk/Cosmo/RadioBeam/fgndsub.cc	(revision 3788)
@@ -0,0 +1,123 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+    Classe  ForegroundCleaner
+    R. Ansari , C. Magneville - Juin 2010 
+---------------------------------------------------------------  */
+
+#include "fgndsub.h"
+#include "lobe.h"
+
+#include "cubedef.h"
+#include "matharr.h"
+
+#include "ctimer.h"
+
+/* --Methode-- */
+ForegroundCleaner::ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube)
+  : arrep_(arrep) , tbeam_(tbeam), skycube_(skycube)
+{
+  double dxdeg = ThetaSizeDegre/(double)NTheta;
+  double dydeg = PhiSizeDegre/(double)NPhi;
+  dx_ = DegreeToRadian(dxdeg);
+  dy_ = DegreeToRadian(dydeg);
+  freq0_ = Freq0MHz;
+  dfreq_ = FreqSizeMHz/(double)NFreq;
+  cout << " ForegroundCleaner: " << " dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx_ << " dy_rad=" << dy_ << ")" 
+       << " Freq0=" << freq0_ << " deltaFreq=" << dfreq_ << " MHz" << endl;
+  skycube.Show();
+}
+
+/* --Methode-- */
+void ForegroundCleaner::BeamCorrections()
+{
+  BeamEffect beam(arrep_);
+  beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_);
+  cout << " ForegroundCleaner::BeamCorrections() done " << endl;
+}
+
+/* --Methode-- */
+int ForegroundCleaner::CleanPointSources(double nsigmas)
+{
+  Timer tm("ForegroundCleaner::CleanPointSources");
+  TArray< TF > sky2d(skycube_.SizeX(), skycube_.SizeY());
+  for(sa_size_t ky=0; ky<sky2d.SizeY(); ky++) 
+    for(sa_size_t kx=0; kx<sky2d.SizeX(); kx++) 
+      sky2d(kx, ky) = skycube_(Range(kx), Range(ky), Range::all()).Sum();
+
+
+  double mean, sigma;
+  
+
+  TArray< TF > amz(1,1,skycube_.SizeZ()), asz(1,1,skycube_.SizeZ());
+  for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++)  {
+    TArray< TF > slice = skycube_(Range::all() , Range::all(), Range(kz));
+    MeanSigma(slice, mean, sigma);
+    amz(0,0,kz)=mean;  asz(0,0,kz)=sigma; 
+  }
+
+  MeanSigma(sky2d, mean, sigma);
+  cout << " ForegroundCleaner::CleanPointSources 2D Sky projection, mean=" << mean << " sigma=" << sigma << endl;  
+  
+  TF seuil = (TF)(mean+nsigmas*sigma);
+  
+  sa_size_t srccnt=0;
+  for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++) 
+    for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) {
+      if (sky2d(kx,ky)>seuil) {
+	srccnt++; 
+	skycube_(Range(kx), Range(ky), Range::all()) = amz;
+      }
+    }
+  
+  cout << " Cleaned NSrc= " << srccnt << " 2D source/pixels (TotNPix=" << sky2d.Size() 
+       << ")-> " << 100.*srccnt/sky2d.Size() <<  "% with S>" << seuil << " NSigmas=" << nsigmas << endl;
+  return (int)srccnt;
+}
+
+
+/* --Methode-- */
+TArray< TF >  ForegroundCleaner::extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx)
+{
+  Timer tm("ForegroundCleaner::extractLSSCube");
+// Inputs : maplss, mapsyc, freq0, dfreq
+// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
+// Return_Array : foreground subtracted LSS signal 
+{
+  sa_size_t sz[5];   sz[0]=skycube_.SizeX();  sz[1]=skycube_.SizeY();
+  synctemp.SetSize(2, sz); 
+  specidx.SetSize(2, sz);
+  TArray<r_4> omap; 
+  omap.SetSize(skycube_, true);
+  Vector vlnf(skycube_.SizeZ());
+  int nprt = 0;
+  // double freq0 : Frequence premier index en k (MHz)
+  // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
+  for(sa_size_t i=0; i<skycube_.SizeX(); i++) 
+    for(sa_size_t j=0; j<skycube_.SizeY(); j++)  {
+      r_8 s1, sx, sx2, sy, sxy;
+      s1=sx=sx2=sy=sxy=0.;
+      for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
+	double lnf=log((double)k*dfreq_+freq0_);
+	vlnf(k)=lnf;
+	double ttot=(r_8)(skycube_(i,j,k));
+	if (ttot < 1.e-5) continue;
+	double lntt=log(ttot);
+	s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
+	sy+=lntt;    sxy+=(lnf*lntt);
+      }
+      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
+      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
+      double T0 = exp(beta+alpha*vlnf(0));
+      if ((i%16==0)&&(j%27==0)) 
+	cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
+	     << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,125) << ")" << endl; 
+      synctemp(i,j) = T0;
+      specidx(i,j) = alpha;
+      for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
+	r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
+	omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+      }
+    }
+  return omap;
+}
+
+}
Index: /trunk/Cosmo/RadioBeam/fgndsub.h
===================================================================
--- /trunk/Cosmo/RadioBeam/fgndsub.h	(revision 3788)
+++ /trunk/Cosmo/RadioBeam/fgndsub.h	(revision 3788)
@@ -0,0 +1,40 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+    Estimation et soustraction d'avant plan radio 
+    R. Ansari , C. Magneville - Juin 2010 
+---------------------------------------------------------------  */
+
+#ifndef FGNDSUB_SEEN
+#define FGNDSUB_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 "lobe.h"
+#include "mdish.h"
+
+#ifndef TF 
+#define TF r_4 
+#endif
+
+class ForegroundCleaner {
+public:
+  ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube);
+  void BeamCorrections();
+  int CleanPointSources(double nsigmas=5.);
+  TArray< TF > extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx);
+
+  Four2DResponse& arrep_;   // Array/Instrument beam response
+  Four2DResponse& tbeam_;    // Target beam to which would be corrected using BeamCorrections()
+  TArray< TF > skycube_;
+  double dx_, dy_;   // taille des pixels (radians) de skycube
+  double freq0_, dfreq_;    // 1ere frequence et bin en frequence de skycube_;
+  
+};
+
+#endif
Index: /trunk/Cosmo/RadioBeam/lobe.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/lobe.cc	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/lobe.cc	(revision 3788)
@@ -17,19 +17,4 @@
 }
 
-/* --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;
-}
 
 
@@ -52,7 +37,50 @@
     ffts.FFTForward(slice, fourAmp);
     conv.setFrequency(f0+kz*df);
-    ApplyLobeK2D(fourAmp, dkx, dky, conv.getLambda());
+    fresp_.setLambda(conv.getLambda());
+    ApplyLobeK2D(fresp_, fourAmp, dkx, dky);
     ffts.FFTBackward(fourAmp, slice, true);
-    if (kz%10==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+    if (kz%20==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  return;
+}
+
+/* --Methode-- */
+void BeamEffect::Correct2RefLobe(Four2DResponse& rep, TArray< TF >& a, double dx, double dy, double f0, double df)
+// dx, dy en radioans, f0, df en MHz
+{
+  Timer tm("BeamEffect::Correct2RefLobe");
+  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);
+    fresp_.setLambda(conv.getLambda());
+    Four2DRespRatio rratio(fresp_, rep);
+    ApplyLobeK2D(rratio, fourAmp, dkx, dky);
+    ffts.FFTBackward(fourAmp, slice, true);
+    if (kz%20==0)  cout << "BeamEffect::Correct2RefLobe() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  return;
+}
+
+/* --Methode-- */
+void BeamEffect::ApplyLobeK2D(Four2DResponse& rep, TArray< complex<TF> >& fourAmp, double dkx, double dky)
+//  dx, dy en radians, lambda en metres 
+{
+  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>(rep(kxx, kyy), 0.);   	
+    }
   }
   return;
Index: /trunk/Cosmo/RadioBeam/lobe.h
===================================================================
--- /trunk/Cosmo/RadioBeam/lobe.h	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/lobe.h	(revision 3788)
@@ -1,5 +1,5 @@
 /*  ------------------------ Projet BAORadio -------------------- 
     Calcul de l'effet de lobe sur Carte2D (angleX,angleY) et 
-    cube 3D (angleX,angleY)
+    cube 3D (angleX,angleY, frquences)
     R. Ansari , C. Magneville - Juin 2010 
 ---------------------------------------------------------------  */
@@ -25,9 +25,20 @@
 class BeamEffect {
 public:
+  // Definition de l'objet avec la reponse en frequence de l'instrument 
   BeamEffect(Four2DResponse& resp);
-  void ApplyLobeK2D(TArray< complex<TF> >& f2d, double dkx, double dky, double lambda=1.);
+  // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement  
   void ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df);  
 
+  // Corrige de l'effet de l'effet de lobe, pour chaque plan de frequence, pour tout ramener au lobe defini par "rep" 
+  void Correct2RefLobe(Four2DResponse& rep, TArray< TF >& a, double dx, double dy, double f0, double df);
+
+  // Applique l'effet de lobe "rep" dans le plan de Fourier  pour une frequence (longueur d'onde) fixee
+  static void ApplyLobeK2D(Four2DResponse& rep, TArray< complex<TF> >& f2d, double dkx, double dky);
+
+  // Re-echntillonnage du cube 3D en appliquant les facteurs xfac,yfac,zfac selon chaque direction 
+  // fac = 2 ---> on double le nombre d'echantillon , fac=0.5 : un echantillon sur deux 
   static TArray< TF > ReSample(TArray< TF >& a, double xfac, double yfac, double zfac);
+
+  // On ajoute du bruit gaussien au cube 3D - espace des positions 
   static void AddNoise(TArray< TF >& a, double pixsignoise, bool fgcmsig=true);
 
Index: /trunk/Cosmo/RadioBeam/makefile
===================================================================
--- /trunk/Cosmo/RadioBeam/makefile	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/makefile	(revision 3788)
@@ -11,6 +11,6 @@
 	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 
+PKGOLIST =  Objs/fgndsub.o Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o 
+PKGHLIST =  fgndsub.h lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h 
 
 ### les executables 
@@ -87,11 +87,14 @@
 
 ### les classes / fonctions 
-Objs/specpk.o : specpk.cc specpk.h  mdish.h qhist.h 
+Objs/fgndsub.o : fgndsub.cc $(PKGHLIST)
+	$(CXXCOMPILE) -o Objs/fgndsub.o fgndsub.cc
+
+Objs/specpk.o : specpk.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/specpk.o specpk.cc
 
-Objs/mdish.o : mdish.cc mdish.h qhist.h 
+Objs/mdish.o : mdish.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/mdish.o mdish.cc
 
-Objs/lobe.o : lobe.cc lobe.h mdish.h qhist.h 
+Objs/lobe.o : lobe.cc $(PKGHLIST)
 	$(CXXCOMPILE) -o Objs/lobe.o lobe.cc
 
Index: /trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/mdish.cc	(revision 3788)
@@ -75,4 +75,21 @@
 }
 
+//---------------------------------------------------------------
+// -- Four2DRespRatio : rapport de la reponse entre deux objets Four2DResponse
+//---------------------------------------------------------------
+Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b)
+  : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b) 
+{
+}
+
+double Four2DRespRatio::Value(double kx, double ky)
+{
+  double ra = a_.Value(kx,ky);
+  double rb = b_.Value(kx,ky);
+  if (rb<1.e-19) rb = 1.e-19;
+  return (ra/rb);
+}
+
+//---------------------------------------------------------------
 //--- Classe simple pour le calcul de rotation 
 class Rotation {
Index: /trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.h	(revision 3787)
+++ /trunk/Cosmo/RadioBeam/mdish.h	(revision 3788)
@@ -55,4 +55,5 @@
 };
 
+
 // -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
 class Four2DRespTable : public  Four2DResponse {
@@ -63,4 +64,15 @@
   virtual double Value(double kx, double ky); 
   Histo2D hrep_;
+};
+
+
+// -- Four2DRespRatio: Retourne le rapport de la reponse entre deux objets Four2DResponse
+class Four2DRespRatio : public  Four2DResponse {
+public:
+  Four2DRespRatio(Four2DResponse& a, Four2DResponse& b);
+  // Return the ratio a.Value(kx,ky) / b.Value(kx, ky) - with protection against divide by zero 
+  virtual double Value(double kx, double ky); 
+  Four2DResponse& a_;
+  Four2DResponse& b_;
 };
 
