#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= (int)mCyl_.size())) { cout <<"******************************************* 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 > & 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 > & 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 > & 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 ((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& z) { return (float)sqrt((double)(z.real()*z.real()+z.imag()*z.imag())); } //----------------------------------------------------------------------------- TMatrix< r_4 > MultiCylinders::getRecXYSlice(int kfmin, int kfmax) { TArray< complex > & 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 MultiCylinders::getRecYXSlice(int kfmin, int kfmax) { TArray< complex > & 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