Changeset 3192 in Sophya for trunk/Cosmo/RadioBeam/mbeamcyl.cc


Ignore:
Timestamp:
Mar 20, 2007, 11:48:15 AM (19 years ago)
Author:
legoff
Message:

input file read with DVcards by treccyl.cc

and by multiCylinders(char* filename)

using real units (meter and GHz)
using _ for member variables (e.g. x_)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/mbeamcyl.cc

    r3191 r3192  
    55#include "srandgen.h"
    66
     7static double cLight=0.3;       // in 1E9 m/s
     8static double tClock = 2.; // should come from param file !!!!
     9//static double cLight=1;       
     10//static double tClock = 1.;
     11
    712//=================================================
    813
    914MultiBeamCyl::MultiBeamCyl(int nr, int ns, double posx, double posy)
    10   : texact(ns) , tjitt(ns) , toffset(nr) ,
    11     signal(ns), sigjitt(ns) , gain(nr)
    12 {
    13   NR = nr;
    14   NS = ns;
    15   posY = posy;
    16   posX = posx;
     15  : texact_(ns) , tjitt_(ns) , toffset_(nr) ,
     16    gain_(nr), signal_(ns), sigjitt_(ns)
     17{
     18  NR_ = nr;
     19  NS_ = ns;
     20  posY_ = posy;
     21  posX_ = posx;
    1722
    1823  SetPrintLevel(0);
    … …  
    2328  SetTimeOffsetSigma(0.);
    2429  SetGains(1., 0., 0);
    25   adfg = false; src = NULL;
     30  adfg_ = false; src_ = NULL;
    2631  SetSources(new BRSourceGen, true);
    2732}
    2833
     34//-----------------------------------------------------------------------------
    2935MultiBeamCyl::~MultiBeamCyl()
    3036{
    31   if (adfg && src) delete src;
    32 }
    33 
     37  if (adfg_ && src_) delete src_;
     38}
     39
     40//-----------------------------------------------------------------------------
    3441void MultiBeamCyl::SetSources(BRSourceGen* brs, bool ad)
    3542{
    3643  if (brs == NULL) return;
    37   if (adfg && src) delete src;
    38   src = brs;  adfg=ad;
    39   if (PrtLev > 1)
     44  if (adfg_ && src_) delete src_;
     45  src_ = brs;  adfg_=ad;
     46  if (PrtLev_ > 1)
    4047    cout << " MultiBeamCyl::SetSources(brs=" << brs <<" ,bool) NbSrc="
    41          << src->NbSources() << endl;
    42 
    43 }
    44 
     48         << src_->NbSources() << endl;
     49
     50}
     51
     52//-----------------------------------------------------------------------------
    4553void MultiBeamCyl::SetGains(double g, double sigg, int nzerogain)
    4654{
    47   if (sigg < 1.e-6)  gain = g;
    48   else gain = RandomSequence(RandomSequence::Gaussian, g, sigg);
     55  if (sigg < 1.e-6)  gain_ = g;
     56  else gain_ = RandomSequence(RandomSequence::Gaussian, g, sigg);
    4957  int k;
    50   for (k=0; k<NR; k++) if (gain(k) < 0) gain(k) = 0.;
     58  for (k=0; k<NR_; k++) if (gain_(k) < 0) gain_(k) = 0.;
    5159  for(k=0; k<nzerogain; k++) {
    52     int zg = random()%NR;
    53     if ((zg >=0) && (zg < NR))  gain(zg) = 0.;
    54   }
    55   if (PrtLev > 1)
     60    int zg = random()%NR_;
     61    if ((zg >=0) && (zg < NR_))  gain_(zg) = 0.;
     62  }
     63  if (PrtLev_ > 1)
    5664    cout << " MultiBeamCyl::SetGains(g=" << g <<" ,sigg=" << sigg
    5765         << " ,nzg=" << nzerogain << " )" << endl;
    5866}
    5967
     68//-----------------------------------------------------------------------------
    6069void MultiBeamCyl::ComputeTimeVectors()
    6170{
    62   texact = RegularSequence(0., 1.);
    63   toffset = RandomSequence(RandomSequence::Gaussian, 0., toffsig);
     71  texact_ = RegularSequence(0., tClock);        // 0, tClock, 2*tClock, ...
     72//  for (int i=0;i<1024;i++) {cout << texact(i)<< endl;};
     73  toffset_ = RandomSequence(RandomSequence::Gaussian, 0., toffsig_);
    6474  NewTJitVector();
    6575}
    6676
     77//-----------------------------------------------------------------------------
    6778void MultiBeamCyl::NewTJitVector(int num)
    6879{
    69   if (timejitter > 1.e-19) {
    70     tjitt = RandomSequence(RandomSequence::Gaussian, 0., timejitter);
    71     tjitt += texact; 
    72   }
    73   else tjitt = texact; 
    74   if (num >= 0) tjitt += toffset(num);
     80  if (timejitter_ > 1.e-19) {
     81    tjitt_ = RandomSequence(RandomSequence::Gaussian, 0., timejitter_);
     82    tjitt_ += texact_; 
     83  }
     84  else tjitt_ = texact_; 
     85  if (num >= 0) tjitt_ += toffset_(num);
    7586 
    7687}
    7788
     89//-----------------------------------------------------------------------------
    7890int MultiBeamCyl::ComputeSignalVector(int num, bool fgsignojit)
    7991{
    8092  int nok = 0;
    81   signal = 0.;
    82   sigjitt = 0.;
    83   for(int is=0; is<src->freq.Size(); is++) {
    84     double fr = src->freq(is);
     93  signal_ = 0.;
     94  sigjitt_ = 0.;
     95  for(int is=0; is<src_->freq.Size(); is++) {
     96    double fr = src_->freq(is);
    8597    if ((fr < 0.) || (fr > 0.5)) continue;
    8698    nok++;
    … …  
    88100    // pas celle apres shift (freq-reduite)
    89101    // lambda = c T = c/freq avec c = 1, dephasage = 2*pi*num*Da*sin(ang)/lambda
    90     double dephasage = (posX+num*Da)*sin(src->angX(is)) +
    91                                 posY*sin(src->angY(is)) ;
    92     dephasage *= (2*M_PI*(fr+freq0));
     102    double dephasage = (posX_+num*Da_)*sin(src_->angX(is))
     103                                + posY_*sin(src_->angY(is)) ;
     104    dephasage *= 2*M_PI*(fr+freq0_)/cLight;
    93105    // On ajoute alors la phase propre de chaque source
    94     dephasage += src->phase(is);
    95     double amprep = src->amp(is)*AngResponse(src->angX(is), src->angY(is));
    96     for(int k=0; k<NS; k++) {
    97       sigjitt(k) += amprep*sin(2.*M_PI*fr*tjitt(k)+dephasage);
     106    dephasage += src_->phase(is);
     107    double amprep = src_->amp(is)*AngResponse(src_->angX(is), src_->angY(is));
     108    for(int k=0; k<NS_; k++) {
     109      sigjitt_(k) += amprep*sin(2.*M_PI*fr*tjitt_(k)+dephasage);
    98110      if (fgsignojit)
    99         signal(k) += amprep*sin(2.*M_PI*fr*texact(k)+dephasage);
     111        signal_(k) += amprep*sin(2.*M_PI*fr*texact_(k)+dephasage);
    100112    }
    101113  }
    102114  // Application du gain du detecteur
    103   r_4 ga = gain(num);
     115  r_4 ga = gain_(num);
    104116  if (fabs(ga-1.) > 1.e-9) {
    105     signal *= ga;
    106     sigjitt *= ga;
     117    signal_ *= ga;
     118    sigjitt_ *= ga;
    107119  }
    108120  // Ajout de bruit (ampli ...)
    109   if (signoise > 1.e-18) {
    110     for(int k=0; k<NS; k++) {
    111       sigjitt(k) += GauRnd(0., signoise);
    112       if (fgsignojit) signal(k) += GauRnd(0., signoise);
    113     }
    114   }
    115 
     121  if (signoise_ > 1.e-18) {
     122    for(int k=0; k<NS_; k++) {
     123      sigjitt_(k) += GauRnd(0., signoise_);
     124      if (fgsignojit) signal_(k) += GauRnd(0., signoise_);
     125    }
     126  }
     127//............      FFT in time
    116128  FFTPackServer ffts;
    117 
    118   ffts.FFTForward(sigjitt, f_sigjit);
    119   if (fgsignojit) ffts.FFTForward(signal, f_sig);
     129  ffts.FFTForward(sigjitt_, f_sigjit_);
     130  if (fgsignojit) ffts.FFTForward(signal_, f_sig_);
    120131
    121132  return nok;
    … …  
    123134
    124135
     136//-----------------------------------------------------------------------------
    125137/*  --- a supprimer ?
    126138inline float myZmodule(complex<r_4>& z)
    … …  
    130142----- */
    131143
     144//-----------------------------------------------------------------------------
    132145void MultiBeamCyl::ReconstructSourcePlane(bool fgzerocentre)
    133146{
    … …  
    135148  int noksrc = ComputeSignalVector(0, false);
    136149  vector<TVector< complex<r_4> > >  vvfc;
    137   cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR
    138        << " PosY=" << posY << " NFreq=" << f_sigjit.Size()-1 << " NOkSrc=" << noksrc << endl;
    139   if (PrtLev > 0) {
    140     cout << " ... posY= " << posY << " MeanGain=" << Mean(gain) << " MeanAmpSrc="
    141          << Mean(src->amp) << endl;
    142     cout << " ...SigNoise=" << signoise << " TimeJitter=" << timejitter
    143          << " TOffSig=" << toffsig << " Da=" << Da << " Freq0=" << freq0 << endl;
     150  cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR_
     151       << " PosY=" << posY_ << " NFreq=" << f_sigjit_.Size()-1 << " NOkSrc=" << noksrc << endl;
     152  if (PrtLev_ > 0) {
     153    cout << " ... posY= " << posY_ << " MeanGain=" << Mean(gain_) << " MeanAmpSrc="
     154         << Mean(src_->amp) << endl;
     155    cout << " ...SigNoise=" << signoise_ << " TimeJitter=" << timejitter_
     156         << " TOffSig=" << toffsig_ << " Da=" << Da_ << " Freq0=" << freq0_ << endl;
    144157  }
    145158  // On ne s'occupe pas de la composante continue
    146   for(int jf=1; jf<f_sigjit.Size(); jf++) {
    147     TVector<complex<r_4> > cf(NR);
     159  for(int jf=1; jf<f_sigjit_.Size(); jf++) {
     160    TVector<complex<r_4> > cf(NR_);
    148161    vvfc.push_back(cf);
    149162  }
    150163  cout << "ReconstructSourcePlane()/Info: computing s(t) for each receptor ..." << endl;
    151   int pmod = NR/10;
    152   for(int ir=0; ir<NR; ir++) {
    153     if (timejitter > 1.e-19) NewTJitVector(ir);
     164  int pmod = NR_/10;
     165  for(int ir=0; ir<NR_; ir++) {
     166    if (timejitter_ > 1.e-19) NewTJitVector(ir);
    154167    noksrc = ComputeSignalVector(ir, false);
    155     for(int jf=1; jf<f_sigjit.Size(); jf++)
    156       vvfc[jf-1](ir) = f_sigjit(jf);
    157     if ( (PrtLev>0) && (ir%pmod == 0) )
    158       cout << " OK s(t) for ir=" << ir << " / NR=" << NR << " NOkSrc=" << noksrc << endl;
     168    for(int jf=1; jf<f_sigjit_.Size(); jf++)
     169      vvfc[jf-1](ir) = f_sigjit_(jf);
     170    if ( (PrtLev_>0) && (ir%pmod == 0) )
     171      cout << " OK s(t) for ir=" << ir << " / NR=" << NR_ << " NOkSrc=" << noksrc << endl;
    159172  }
    160173
    161174  cout << "ReconstructSourcePlane()/Info: computing s(ang) for each freq by FFT" << endl;
    162   cmplx_srcplane.SetSize(f_sigjit.Size()-1, NR);
    163   // rec_srcplane.SetSize(f_sigjit.Size()-1, NR);
     175  cmplx_srcplane_.SetSize(f_sigjit_.Size()-1, NR_);
     176  // rec_srcplane.SetSize(f_sigjit_.Size()-1, NR_);
    164177  FFTPackServer ffts;
    165178  TVector<complex<r_4> > fcf;
    166179  pmod = vvfc.size()/10;
    167   for(int jf=0; jf<vvfc.size(); jf++) {
    168     ffts.FFTForward(vvfc[jf], fcf);
    169     if (fcf.Size() != NR) {
     180  for(int jf=0; jf<(int)vvfc.size(); jf++) {            // loop over frequencies
     181    ffts.FFTForward(vvfc[jf], fcf);     // FFT alomg cylinder
     182    if (fcf.Size() != NR_) {
    170183      cout << "ReconstructSourcePlane()/BUG jf=" << jf << " fcf.Size() != NR "
    171            << fcf.Size() << " != " << NR << endl;
     184           << fcf.Size() << " != " << NR_ << endl;
    172185      continue;
    173186    }
    174187    if (fgzerocentre) {  // On veut avoir la direction angle=0 au milieu de la matrice
    175       int milieu = (NR-1)/2;
    176       if (NR%2 == 0) milieu++;
    177       int decal = NR - milieu - 1;
     188      int milieu = (NR_-1)/2;
     189      if (NR_%2 == 0) milieu++;
     190      int decal = NR_ - milieu - 1;
    178191      for (int ir=0; ir<=milieu; ir++)
    179         cmplx_srcplane(jf, ir+decal) = fcf(ir);
     192        cmplx_srcplane_(jf, ir+decal) = fcf(ir);
    180193      //        rec_srcplane(jf, ir+decal) = myZmodule(fcf(ir));
    181194
    182       for (int ir=milieu+1; ir<NR; ir++)
    183         cmplx_srcplane(jf, decal-(NR-ir)) = fcf(ir);
    184       //        rec_srcplane(jf, decal-(NR-ir)) = myZmodule(fcf(ir));
     195      for (int ir=milieu+1; ir<NR_; ir++)
     196        cmplx_srcplane_(jf, decal-(NR_-ir)) = fcf(ir);
     197      //        rec_srcplane_(jf, decal-(NR_-ir)) = myZmodule(fcf(ir));
    185198     
    186199    }
    187     else for (int ir=0; ir<NR; ir++)
    188       cmplx_srcplane(jf, ir) = fcf(ir);
     200    else for (int ir=0; ir<NR_; ir++)
     201      cmplx_srcplane_(jf, ir) = fcf(ir);
    189202    //      rec_srcplane(jf, ir) = myZmodule(fcf(ir));
    190203
    191     if ( (PrtLev > 0) && (jf%pmod == 0))
     204    if ( (PrtLev_ > 0) && (jf%pmod == 0))
    192205      cout << " OK rec_srcplane(jf, ir) for jf=" << jf << endl;
    193206  }
Note: See TracChangeset for help on using the changeset viewer.