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

#include "mbeamcyl.h"
#include "fftpserver.h"
#include "vector3d.h"
#include "matharr.h"
#include "srandgen.h"

//=================================================

MultiBeamCyl::MultiBeamCyl(int nr, int ns, double posy)
  : texact(ns) , tjitt(ns) , toffset(nr) ,
    signal(ns), sigjitt(ns) , gain(nr)
{
  NR = nr;
  NS = ns;
  posY = posy;

  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., 1.);
  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 = num*Da*sin(src->angX(is)) + posY*sin(src->angY(is)) ;
    dephasage *= (2*M_PI*(fr+freq0));
    // 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);
    }
  }

  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<vvfc.size(); jf++) {
    ffts.FFTForward(vvfc[jf], fcf);
    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;
}