source: Sophya/trunk/Cosmo/RadioBeam/mbeamcyl.cc@ 3583

#include "mbeamcyl.h"
#include "fftpserver.h"
#include "vector3d.h"
#include "matharr.h"
#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) ,
    gain_(nr), signal_(ns), sigjitt_(ns) 
{
  NR_ = nr;
  NS_ = ns;
  posY_ = posy;
  posX_ = posx;

  SetPrintLevel(0);

  SetBaseFreqDa(2., 0.25);
  SetNoiseSigma(0.);
  SetTimeJitter(0.);
  SetTimeOffsetSigma(0.);
  SetGains(1., 0., 0);
  adfg_ = false; src_ = NULL;
  SetSources(new BRSourceGen, true);
}

//-----------------------------------------------------------------------------
MultiBeamCyl::~MultiBeamCyl()
{
  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) 
    cout << " MultiBeamCyl::SetSources(brs=" << brs <<" ,bool) NbSrc=" 
         << 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);
  int k;
  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) 
    cout << " MultiBeamCyl::SetGains(g=" << g <<" ,sigg=" << sigg 
         << " ,nzg=" << nzerogain << " )" << endl;
}

//-----------------------------------------------------------------------------
void MultiBeamCyl::ComputeTimeVectors()
{
  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);
  
}

//-----------------------------------------------------------------------------
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);
    if ((fr < 0.) || (fr > 0.5)) continue;
    nok++;
    // Pour le dephasage entre recepteurs, on doit utiliser la frequence vraie,
    // 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_)/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);
      if (fgsignojit) 
        signal_(k) += amprep*sin(2.*M_PI*fr*texact_(k)+dephasage);
    }
  }
  // Application du gain du detecteur 
  r_4 ga = gain_(num);
  if (fabs(ga-1.) > 1.e-9) {
    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_);
    }
  }
//............      FFT in time
  FFTPackServer ffts;
  ffts.FFTForward(sigjitt_, f_sigjit_);
  if (fgsignojit) ffts.FFTForward(signal_, f_sig_);

  return nok;
}


//-----------------------------------------------------------------------------
/*  --- a supprimer ?
inline float myZmodule(complex<r_4>& z) 
{
  return (float)sqrt((double)(z.real()*z.real()+z.imag()*z.imag()));
}
----- */

//-----------------------------------------------------------------------------
void MultiBeamCyl::ReconstructSourcePlane(bool fgzerocentre)
{
  ComputeTimeVectors();
  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;
  }
  // On ne s'occupe pas de la composante continue
  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);
    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;
  }

  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_);
  FFTPackServer ffts;
  TVector<complex<r_4> > fcf;
  pmod = vvfc.size()/10;
  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;
      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;
      for (int ir=0; ir<=milieu; 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));
      
    }
    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)) 
      cout << " OK rec_srcplane(jf, ir) for jf=" << jf << endl;
  }

  
  cout << "ReconstructSourcePlane()/Info: rec_srcplane computed OK" << endl;
}