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

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:

• trunk/Cosmo/RadioBeam/multicyl.cc (modified) (8 diffs)

trunk/Cosmo/RadioBeam/multicyl.cc

@@ -6 +6 @@
 #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.; 
 
 
@@ -12 +18 @@
 MultiCylinders::MultiCylinders(int nr, int ns)
 {
-  NR = nr;
-  NS = ns;
+  NR_ = nr;
+  NS_ = ns;
 
   SetPrintLevel(0);
@@ -22 +28 @@
   SetTimeOffsetSigma(0.);
   SetGains(1., 0., 0);
-  adfg = false; src = NULL;
+  adfg_ = false; src_ = NULL;
   SetSources(new BRSourceGen, true);
 }
 
+//-----------------------------------------------------------------------------
+MultiCylinders::MultiCylinders(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<mCyl.size(); k++) delete mCyl[k];
-}
-
+  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 >= mCyl.size()))
+  if ((k < 0) || (k >= (int)mCyl_.size())) {
+    cout <<"******************************************* k="<<k<<endl;
     throw RangeCheckError("MultiCylinders::GetCylinder(k) k<0 OR k>=NCyl");
-  return (*mCyl[k]);
-}
-
+  }
+  return (*mCyl_[k]);
+}
+
+//-----------------------------------------------------------------------------
 void MultiCylinders::SetSources(BRSourceGen* brs, bool ad)
 {
   if (brs == NULL) return;
-  if (adfg && src) delete src;
-  src = brs;  adfg=ad;
-}
-
-
+  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);
-  }
-}
-
-
+  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() 
+  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - NCyl= " << mCyl_.size() 
        << " MemSize=" << resu.getMemorySize() << " kb" << endl;       
   ConfigureCylinders();
 
   char buff[128];
-  for(int k=0; k<mCyl.size(); k++) {
+  for(int k=0; k<(int)mCyl_.size(); k++) {
     cout << "---- Cyl[" << k << "] Calling MultiBeamCyl.ReconstructSourcePlane() ..." << endl;
-    mCyl[k]->ReconstructSourcePlane(fgzerocentre);
+    mCyl_[k]->ReconstructSourcePlane(fgzerocentre);
     sprintf(buff,"Cyl[%d].RecSrcPlane()",k);
     tm.Split(buff);
@@ -82 +139 @@
 }
 
-void MultiCylinders::ReconstructSourceBox(int halfny, double stepangy)
-{
+//-----------------------------------------------------------------------------
+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
+// 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};
-//  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 );
+  box.ReSize(3, sz);            //  values initialized to zero ?
   cout << " MultiCylinders::ReconstructSourceBox(" << halfny << "," << stepangy
        << "=" << Angle(stepangy).ToArcMin() << " ) srcbox:" << endl;
@@ -102 +162 @@
 
 
-  int pmod = mtx.NRows()/10;
-  
-  double fstep = 1.0/(double)NS;  // pas en frequence, attention, on a vire la composante continu
+//  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
-    
-    for(int lc=0; lc<mCyl.size(); lc++) {  // Loop over cylinders 
+    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 = - 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;
+      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 * 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++;
+      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%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;
+         << " / Max_kf=" << mtx.NRows() << endl;
   } // End of loop over over frequencies
 
@@ -140 +204 @@
 } 
 
+//-----------------------------------------------------------------------------
 inline float myZmodule(complex<r_4>& z) 
 {
@@ -145 +210 @@
 }
 
+//-----------------------------------------------------------------------------
 TMatrix< r_4 >  MultiCylinders::getRecXYSlice(int kfmin, int kfmax)
 {
@@ -158 +224 @@
   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);
+}