Index: trunk/Cosmo/RadioBeam/mbeamcyl.cc
===================================================================
--- trunk/Cosmo/RadioBeam/mbeamcyl.cc	(revision 3191)
+++ trunk/Cosmo/RadioBeam/mbeamcyl.cc	(revision 3192)
@@ -5,14 +5,19 @@
 #include "srandgen.h"
 
+static double cLight=0.3;	// in 1E9 m/s
+static double tClock = 2.; // should come from param file !!!!
+//static double cLight=1;	
+//static double tClock = 1.; 
+
 //=================================================
 
 MultiBeamCyl::MultiBeamCyl(int nr, int ns, double posx, double posy)
-  : texact(ns) , tjitt(ns) , toffset(nr) ,
-    signal(ns), sigjitt(ns) , gain(nr)
-{
-  NR = nr;
-  NS = ns;
-  posY = posy;
-  posX = posx;
+  : texact_(ns) , tjitt_(ns) , toffset_(nr) ,
+    gain_(nr), signal_(ns), sigjitt_(ns) 
+{
+  NR_ = nr;
+  NS_ = ns;
+  posY_ = posy;
+  posX_ = posx;
 
   SetPrintLevel(0);
@@ -23,64 +28,71 @@
   SetTimeOffsetSigma(0.);
   SetGains(1., 0., 0);
-  adfg = false; src = NULL;
+  adfg_ = false; src_ = NULL;
   SetSources(new BRSourceGen, true);
 }
 
+//-----------------------------------------------------------------------------
 MultiBeamCyl::~MultiBeamCyl()
 {
-  if (adfg && src) delete src;
-}
-
+  if (adfg_ && src_) delete src_;
+}
+
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::SetSources(BRSourceGen* brs, bool ad)
 {
   if (brs == NULL) return;
-  if (adfg && src) delete src;
-  src = brs;  adfg=ad;
-  if (PrtLev > 1) 
+  if (adfg_ && src_) delete src_;
+  src_ = brs;  adfg_=ad;
+  if (PrtLev_ > 1) 
     cout << " MultiBeamCyl::SetSources(brs=" << brs <<" ,bool) NbSrc=" 
-	 << src->NbSources() << endl;
-
-}
-
+	 << src_->NbSources() << endl;
+
+}
+
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::SetGains(double g, double sigg, int nzerogain)
 {
-  if (sigg < 1.e-6)  gain = g;
-  else gain = RandomSequence(RandomSequence::Gaussian, g, sigg);
+  if (sigg < 1.e-6)  gain_ = g;
+  else gain_ = RandomSequence(RandomSequence::Gaussian, g, sigg);
   int k;
-  for (k=0; k<NR; k++) if (gain(k) < 0) gain(k) = 0.;
+  for (k=0; k<NR_; k++) if (gain_(k) < 0) gain_(k) = 0.;
   for(k=0; k<nzerogain; k++) {
-    int zg = random()%NR;
-    if ((zg >=0) && (zg < NR))  gain(zg) = 0.; 
-  }
-  if (PrtLev > 1) 
+    int zg = random()%NR_;
+    if ((zg >=0) && (zg < NR_))  gain_(zg) = 0.; 
+  }
+  if (PrtLev_ > 1) 
     cout << " MultiBeamCyl::SetGains(g=" << g <<" ,sigg=" << sigg 
 	 << " ,nzg=" << nzerogain << " )" << endl;
 }
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::ComputeTimeVectors()
 {
-  texact = RegularSequence(0., 1.);
-  toffset = RandomSequence(RandomSequence::Gaussian, 0., toffsig);
+  texact_ = RegularSequence(0., tClock);	// 0, tClock, 2*tClock, ...
+//  for (int i=0;i<1024;i++) {cout << texact(i)<< endl;};
+  toffset_ = RandomSequence(RandomSequence::Gaussian, 0., toffsig_);
   NewTJitVector();
 }
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::NewTJitVector(int num)
 {
-  if (timejitter > 1.e-19) {
-    tjitt = RandomSequence(RandomSequence::Gaussian, 0., timejitter);
-    tjitt += texact;  
-  }
-  else tjitt = texact;  
-  if (num >= 0) tjitt += toffset(num);
+  if (timejitter_ > 1.e-19) {
+    tjitt_ = RandomSequence(RandomSequence::Gaussian, 0., timejitter_);
+    tjitt_ += texact_;  
+  }
+  else tjitt_ = texact_;  
+  if (num >= 0) tjitt_ += toffset_(num);
   
 }
 
+//-----------------------------------------------------------------------------
 int MultiBeamCyl::ComputeSignalVector(int num, bool fgsignojit)
 {
   int nok = 0;
-  signal = 0.;
-  sigjitt = 0.;
-  for(int is=0; is<src->freq.Size(); is++) {
-    double fr = src->freq(is);
+  signal_ = 0.;
+  sigjitt_ = 0.;
+  for(int is=0; is<src_->freq.Size(); is++) {
+    double fr = src_->freq(is);
     if ((fr < 0.) || (fr > 0.5)) continue;
     nok++;
@@ -88,34 +100,33 @@
     // pas celle apres shift (freq-reduite)
     // lambda = c T = c/freq avec c = 1, dephasage = 2*pi*num*Da*sin(ang)/lambda
-    double dephasage = (posX+num*Da)*sin(src->angX(is)) + 
-    				posY*sin(src->angY(is)) ;
-    dephasage *= (2*M_PI*(fr+freq0));
+    double dephasage = (posX_+num*Da_)*sin(src_->angX(is)) 
+    				+ posY_*sin(src_->angY(is)) ;
+    dephasage *= 2*M_PI*(fr+freq0_)/cLight;
     // On ajoute alors la phase propre de chaque source
-    dephasage += src->phase(is);
-    double amprep = src->amp(is)*AngResponse(src->angX(is), src->angY(is));
-    for(int k=0; k<NS; k++) {
-      sigjitt(k) += amprep*sin(2.*M_PI*fr*tjitt(k)+dephasage);
+    dephasage += src_->phase(is);
+    double amprep = src_->amp(is)*AngResponse(src_->angX(is), src_->angY(is));
+    for(int k=0; k<NS_; k++) {
+      sigjitt_(k) += amprep*sin(2.*M_PI*fr*tjitt_(k)+dephasage);
       if (fgsignojit) 
-	signal(k) += amprep*sin(2.*M_PI*fr*texact(k)+dephasage);
+	signal_(k) += amprep*sin(2.*M_PI*fr*texact_(k)+dephasage);
     }
   }
   // Application du gain du detecteur 
-  r_4 ga = gain(num);
+  r_4 ga = gain_(num);
   if (fabs(ga-1.) > 1.e-9) {
-    signal *= ga;
-    sigjitt *= ga;
+    signal_ *= ga;
+    sigjitt_ *= ga;
   }
   // Ajout de bruit (ampli ...) 
-  if (signoise > 1.e-18) {
-    for(int k=0; k<NS; k++) {
-      sigjitt(k) += GauRnd(0., signoise);
-      if (fgsignojit) signal(k) += GauRnd(0., signoise);
-    }
-  }
-
+  if (signoise_ > 1.e-18) {
+    for(int k=0; k<NS_; k++) {
+      sigjitt_(k) += GauRnd(0., signoise_);
+      if (fgsignojit) signal_(k) += GauRnd(0., signoise_);
+    }
+  }
+//............      FFT in time
   FFTPackServer ffts;
-
-  ffts.FFTForward(sigjitt, f_sigjit);
-  if (fgsignojit) ffts.FFTForward(signal, f_sig);
+  ffts.FFTForward(sigjitt_, f_sigjit_);
+  if (fgsignojit) ffts.FFTForward(signal_, f_sig_);
 
   return nok;
@@ -123,4 +134,5 @@
 
 
+//-----------------------------------------------------------------------------
 /*  --- a supprimer ?
 inline float myZmodule(complex<r_4>& z) 
@@ -130,4 +142,5 @@
 ----- */
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::ReconstructSourcePlane(bool fgzerocentre)
 {
@@ -135,59 +148,59 @@
   int noksrc = ComputeSignalVector(0, false);
   vector<TVector< complex<r_4> > >  vvfc;
-  cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR 
-       << " PosY=" << posY << " NFreq=" << f_sigjit.Size()-1 << " NOkSrc=" << noksrc << endl;
-  if (PrtLev > 0) {
-    cout << " ... posY= " << posY << " MeanGain=" << Mean(gain) << " MeanAmpSrc=" 
-	 << Mean(src->amp) << endl;
-    cout << " ...SigNoise=" << signoise << " TimeJitter=" << timejitter 
-	 << " TOffSig=" << toffsig << " Da=" << Da << " Freq0=" << freq0 << endl;
+  cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR_ 
+       << " PosY=" << posY_ << " NFreq=" << f_sigjit_.Size()-1 << " NOkSrc=" << noksrc << endl;
+  if (PrtLev_ > 0) {
+    cout << " ... posY= " << posY_ << " MeanGain=" << Mean(gain_) << " MeanAmpSrc=" 
+	 << Mean(src_->amp) << endl;
+    cout << " ...SigNoise=" << signoise_ << " TimeJitter=" << timejitter_ 
+	 << " TOffSig=" << toffsig_ << " Da=" << Da_ << " Freq0=" << freq0_ << endl;
   }
   // On ne s'occupe pas de la composante continue
-  for(int jf=1; jf<f_sigjit.Size(); jf++) {
-    TVector<complex<r_4> > cf(NR);
+  for(int jf=1; jf<f_sigjit_.Size(); jf++) {
+    TVector<complex<r_4> > cf(NR_);
     vvfc.push_back(cf);
   }
   cout << "ReconstructSourcePlane()/Info: computing s(t) for each receptor ..." << endl;
-  int pmod = NR/10;
-  for(int ir=0; ir<NR; ir++) {
-    if (timejitter > 1.e-19) NewTJitVector(ir);
+  int pmod = NR_/10;
+  for(int ir=0; ir<NR_; ir++) {
+    if (timejitter_ > 1.e-19) NewTJitVector(ir);
     noksrc = ComputeSignalVector(ir, false);
-    for(int jf=1; jf<f_sigjit.Size(); jf++)
-      vvfc[jf-1](ir) = f_sigjit(jf);
-    if ( (PrtLev>0) && (ir%pmod == 0) )
-      cout << " OK s(t) for ir=" << ir << " / NR=" << NR << " NOkSrc=" << noksrc << endl;
+    for(int jf=1; jf<f_sigjit_.Size(); jf++)
+      vvfc[jf-1](ir) = f_sigjit_(jf);
+    if ( (PrtLev_>0) && (ir%pmod == 0) )
+      cout << " OK s(t) for ir=" << ir << " / NR=" << NR_ << " NOkSrc=" << noksrc << endl;
   }
 
   cout << "ReconstructSourcePlane()/Info: computing s(ang) for each freq by FFT" << endl;
-  cmplx_srcplane.SetSize(f_sigjit.Size()-1, NR);
-  // rec_srcplane.SetSize(f_sigjit.Size()-1, NR);
+  cmplx_srcplane_.SetSize(f_sigjit_.Size()-1, NR_);
+  // rec_srcplane.SetSize(f_sigjit_.Size()-1, NR_);
   FFTPackServer ffts;
   TVector<complex<r_4> > fcf;
   pmod = vvfc.size()/10;
-  for(int jf=0; jf<vvfc.size(); jf++) {
-    ffts.FFTForward(vvfc[jf], fcf);
-    if (fcf.Size() != NR) {
+  for(int jf=0; jf<(int)vvfc.size(); jf++) {		// loop over frequencies
+    ffts.FFTForward(vvfc[jf], fcf);	// FFT alomg cylinder
+    if (fcf.Size() != NR_) {
       cout << "ReconstructSourcePlane()/BUG jf=" << jf << " fcf.Size() != NR "
-	   << fcf.Size() << " != " << NR << endl;
+	   << fcf.Size() << " != " << NR_ << endl;
       continue;
     }
     if (fgzerocentre) {  // On veut avoir la direction angle=0 au milieu de la matrice
-      int milieu = (NR-1)/2;
-      if (NR%2 == 0) milieu++;
-      int decal = NR - milieu - 1;
+      int milieu = (NR_-1)/2;
+      if (NR_%2 == 0) milieu++;
+      int decal = NR_ - milieu - 1;
       for (int ir=0; ir<=milieu; ir++) 
-	cmplx_srcplane(jf, ir+decal) = fcf(ir);
+	cmplx_srcplane_(jf, ir+decal) = fcf(ir);
       //	rec_srcplane(jf, ir+decal) = myZmodule(fcf(ir));
 
-      for (int ir=milieu+1; ir<NR; ir++) 
-	cmplx_srcplane(jf, decal-(NR-ir)) = fcf(ir);
-      //	rec_srcplane(jf, decal-(NR-ir)) = myZmodule(fcf(ir));
+      for (int ir=milieu+1; ir<NR_; ir++) 
+	cmplx_srcplane_(jf, decal-(NR_-ir)) = fcf(ir);
+      //	rec_srcplane_(jf, decal-(NR_-ir)) = myZmodule(fcf(ir));
       
     }
-    else for (int ir=0; ir<NR; ir++) 
-      cmplx_srcplane(jf, ir) = fcf(ir);
+    else for (int ir=0; ir<NR_; ir++) 
+      cmplx_srcplane_(jf, ir) = fcf(ir);
     //      rec_srcplane(jf, ir) = myZmodule(fcf(ir));
 
-    if ( (PrtLev > 0) && (jf%pmod == 0)) 
+    if ( (PrtLev_ > 0) && (jf%pmod == 0)) 
       cout << " OK rec_srcplane(jf, ir) for jf=" << jf << endl;
   }
