source: Sophya/trunk/Cosmo/RadioBeam/multicyl.cc@ 3618

#include "multicyl.h"
#include "fftpserver.h"
#include "vector3d.h"
#include "matharr.h"
#include "srandgen.h"
#include "ctimer.h"
#include "resusage.h"
#include "datacards.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.; 


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

MultiCylinders::MultiCylinders(int nr, int ns)
{
  NR_ = nr;
  NS_ = ns;

  SetPrintLevel(0);

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

//-----------------------------------------------------------------------------
MultiCylinders::MultiCylinders(const char* fileName)
{
  adfg_ = false; src_ = NULL;
  SetSources(new BRSourceGen, true);
  
//      read telescope parameters and fill variable members
  DataCards dc;
  dc.ReadFile(fileName);
//      in old versions frequences were in units of 1/T = 0.5 GHz
//      and distances in units of cT =3E8 * 2E-9=0.60 m
//  double fUnit=0.5;   // 0.5 GHz <=> T = 2 ns
//  double dUnit=0.6;   // distance unit in m.
//      now f and d in real units 
  double fUnit=1.;      // 0.5 GHz <=> T = 2 ns
  double dUnit=1.;      // distance unit in m.
  
  NR_=dc.IParam("nAntenna");
  NS_=dc.IParam("nSample");
  PrtLev_=dc.IParam("printLevel");
  Da_=dc.DParam("dAntenna")/dUnit;
  freq0_=dc.DParam("freq0")/fUnit;
  timejitter_=dc.DParam("sigmaTimeJitt");
  toffsig_=dc.DParam("sigmaClockJitt");
  signoise_=dc.DParam("noiseSigma");
  gain_=dc.DParam("meanGain");
  siggain_=dc.DParam("sigmaGain");
  ngainzero_=dc.IParam("nDeadAntenna");
  
//  tClock=dc.DParam("tClock");
  int nCyl=dc.IParam("nCyl");
  for (int i=0; i<nCyl; i++){
    double xCyl=dc.DParam("xCyl",i)/dUnit;
    double yCyl=dc.DParam("yCyl",i)/dUnit;
    AddCylinder(xCyl,yCyl);
  }
//  maxangX=dc.DParam("angMaxX");
//  double cylDiam=dc.DParam("cylinderDiam")/dUnit;
//// thetaMax = lambda_M/d = c/freq_min/d;  freq_min = freq0 + 1/2T
//  maxangY=cLight/(freq0+1./2./tClock)/cylDiam;
//  halfNY=dc.IParam("halfNY");
//  NX=dc.IParam("NX");
}

//-----------------------------------------------------------------------------
MultiCylinders::~MultiCylinders()
{
  if (adfg_ && src_) delete src_;
  for(int k=0; k<(int)mCyl_.size(); k++) delete mCyl_[k];
}

//-----------------------------------------------------------------------------
MultiBeamCyl& MultiCylinders::GetCylinder(int k)
{
  if ((k < 0) || (k >= (int)mCyl_.size())) {
    cout <<"******************************************* k="<<k<<endl;
    throw RangeCheckError("MultiCylinders::GetCylinder(k) k<0 OR k>=NCyl");
  }
  return (*mCyl_[k]);
}

//-----------------------------------------------------------------------------
void MultiCylinders::SetSources(BRSourceGen* brs, bool ad)
{
  if (brs == NULL) return;
  if (adfg_ && src_) delete src_;
  src_ = brs;  adfg_=ad;
}


//-----------------------------------------------------------------------------
void MultiCylinders::ConfigureCylinders()
{
  cout << " MultiCylinders::ConfigureCylinders() with NCyl= " << mCyl_.size() << endl;
  for(int k=0; k<(int)mCyl_.size(); k++) {
    mCyl_[k]->SetPrintLevel(PrtLev_);
    mCyl_[k]->SetBaseFreqDa(freq0_, Da_);
    mCyl_[k]->SetNoiseSigma(signoise_);
    mCyl_[k]->SetTimeJitter(timejitter_);
    mCyl_[k]->SetTimeOffsetSigma(toffsig_);
    mCyl_[k]->SetGains(gain_, siggain_, ngainzero_);
    mCyl_[k]->SetSources(src_, false);
  }
}


//-----------------------------------------------------------------------------
void MultiCylinders::ReconstructCylinderPlaneS(bool fgzerocentre)
{
  Timer tm("RecCylPlaneS");
  ResourceUsage resu;
  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - NCyl= " << mCyl_.size() 
       << " MemSize=" << resu.getMemorySize() << " kb" << endl;       
  ConfigureCylinders();

  char buff[128];
  for(int k=0; k<(int)mCyl_.size(); k++) {
    cout << "---- Cyl[" << k << "] Calling MultiBeamCyl.ReconstructSourcePlane() ..." << endl;
    mCyl_[k]->ReconstructSourcePlane(fgzerocentre);
    sprintf(buff,"Cyl[%d].RecSrcPlane()",k);
    tm.Split(buff);
  }
  resu.Update();
  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - Done " 
       << " MemSize=" << resu.getMemorySize() << " kb" << endl;
}

//-----------------------------------------------------------------------------
void MultiCylinders::ReconstructSourceBox(int halfny, double stepangy, 
        int nx, double stepangx)
{
  nx=256;
  ReconstructCylinderPlaneS(true);
  TMatrix< complex<r_4> > & mtx = GetCylinder(0).getRecSrcPlane(); 
// boite 3D X:angX, Y:angY , Z: freq
//      if all cylinders at same x position NX is set to zero
//      => x-size = mtx.NCols() = numbers of receptors per cylinders
//              move that to readParam() ?
  sa_size_t sz[5] = {0,0,0,0,0};
  sz[0] = nx;
  sz[1] = halfny*2+1;
  sz[2] = mtx.NRows();
  TArray< complex<r_4> > & box = getRecSrcBox();
  box.ReSize(3, sz);            //  values initialized to zero ?
  cout << " MultiCylinders::ReconstructSourceBox(" << halfny << "," << stepangy
       << "=" << Angle(stepangy).ToArcMin() << " ) srcbox:" << endl;
  Timer tm("RecSrcBox");
  box.Show(cout);


//  int pmod = mtx.NRows()/10;
  
  double fstep = 1.0/(double)NS_/tClock;  // pas en frequence, 
                                // attention, on a vire la composante continu
//  bool first=true;
  for(int kf=0; kf<mtx.NRows(); kf++) { // Loop over frequencies
    double frq = (double)(kf+1.)*fstep + freq0_;  // + 1 car f=0 a ete vire
//            then up to frq =   (mtx.NRows() +1 ) * fstep            !!?
//    cout<<"************"<<mCyl.size()
    for(int lc=0; lc<(int)mCyl_.size(); lc++) {  // Loop over cylinders 
      MultiBeamCyl& cyl = GetCylinder(lc);

      TMatrix< complex<r_4> > & recp = cyl.getRecSrcPlane(); 

      double facl_y = - 2*M_PI*frq*cyl.getCylinderYPos()/cLight;  // attention signe -
      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da_; 
      complex< r_4 > phasefactor;
      int jyy = 0;
      
      for(int jy=-halfny; jy<=halfny; jy++) {  // Loop over Y angular steps
      for(int ix=0; ix<nx; ix++) {  // Loop over AngX directions
        double dphi = facl_y * sin( (double)jy*stepangy )
                      + facl_x*(double(ix)/double(nx)-1./2.);
        phasefactor = complex< r_4 > ((r_4) cos(dphi) , (r_4) sin(dphi));
        // sur recp : index ligne -> frequence  ,    index colonne -> angX ,
        int ixx=(int)(ix*(double)cyl.NR_/double(nx));
        box(ix, jyy, kf) += recp(kf, ixx)*phasefactor;
      }  //  
      jyy++;
      } // End of Loop over Y angles
    } // End of loop over cylinders

    tm.Nop();
//    if ( (PrtLev_>0) && (kf%pmod == 0) )   
    if ( (PrtLev_>0) && (kf%(mtx.NRows()/10) == 0) )   
      cout << " OK box(angx,angy, freq=kf) done for kf=" << kf 
         << " / Max_kf=" << mtx.NRows() << endl;
  } // End of loop over over frequencies

  cout << " MultiCylinders::ReconstructSourceBox() done " << endl; 
} 

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

//-----------------------------------------------------------------------------
TMatrix< r_4 >  MultiCylinders::getRecXYSlice(int kfmin, int kfmax)
{
  TArray< complex<r_4> > & box = getRecSrcBox();  
  if ((kfmin < 0) || (kfmin >=  box.SizeZ()) || (kfmax < kfmin) || (kfmax >=  box.SizeZ()) )
    throw RangeCheckError("MultiCylinders::getRecXYSlice(kfmin, kfmax)");
  TMatrix< r_4> rmtx(box.SizeY(), box.SizeX());
  for(int kf=kfmin; kf<=kfmax; kf++) {
    for(int jy=0; jy<box.SizeY(); jy++) 
      for(int ix=0; ix<box.SizeX(); ix++) 
        rmtx(jy, ix) += myZmodule(box(ix, jy, kf));
  }
  return(rmtx);
}
//-----------------------------------------------------------------------------
TMatrix< r_4 >  MultiCylinders::getRecYXSlice(int kfmin, int kfmax)
{
  TArray< complex<r_4> > & box = getRecSrcBox();  
  if ((kfmin < 0) || (kfmin >=  box.SizeZ()) || (kfmax < kfmin) || (kfmax >=  box.SizeZ()) )
    throw RangeCheckError("MultiCylinders::getRecXYSlice(kfmin, kfmax)");
  TMatrix< r_4> rmtx(box.SizeX(), box.SizeY());
  for(int kf=kfmin; kf<=kfmax; kf++) {
    for(int jy=0; jy<box.SizeY(); jy++) 
      for(int ix=0; ix<box.SizeX(); ix++) 
        rmtx(ix, jy) += myZmodule(box(ix, jy, kf));
  }
  return(rmtx);
}