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

#include "multicyl.h"
#include "fftpserver.h"
#include "vector3d.h"
#include "matharr.h"
#include "srandgen.h"
#include "ctimer.h"
#include "resusage.h"


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

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()
{
  if (adfg && src) delete src;
  for(int k=0; k<mCyl.size(); k++) delete mCyl[k];
}

MultiBeamCyl& MultiCylinders::GetCylinder(int k)
{
  if ((k < 0) || (k >= mCyl.size()))
    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<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<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)
{
  ReconstructCylinderPlaneS(true);
  TMatrix< complex<r_4> > & mtx = GetCylinder(0).getRecSrcPlane(); 
  // boite 3D X:angX, Y:angY , Z: freq
  sa_size_t sz[5] = {0,0,0,0,0};
//  int nx=mtx.NCols(); // uncomment to go back to nxbin = nantenna
  int nx=256;
  sz[0] = nx;
  sz[1] = halfny*2+1;
  sz[2] = mtx.NRows();
  TArray< complex<r_4> > & box = getRecSrcBox();
  box.ReSize(3, sz);
  //  box = complex< r_4 >( 0.f, 0.f );
  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;  // pas en frequence, attention, on a vire la composante continu
  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
    
    for(int lc=0; lc<mCyl.size(); lc++) {  // Loop over cylinders 
      MultiBeamCyl& cyl = GetCylinder(lc);

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

      double facl = - 2*M_PI*frq*cyl.getCylinderYPos();  // attention au signe -
//      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da/(double)cyl.NR; 
      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da; 
      double dphi;
      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 * 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) )   
      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);
}