source: Sophya/trunk/Cosmo/RadioBeam/treccyl.cc@ 3167

#include "sopnamsp.h"
#include "machdefs.h"
#include <math.h>
#include <iostream>
#include <typeinfo>

#include "tvector.h"
#include "srandgen.h"
#include "fioarr.h"
#include "sopemtx.h"
#include "pexceptions.h"
#include "matharr.h"

#include "sambainit.h"

// #include "tarrinit.h"

#include "timing.h"

#include "multicyl.h"
#include "mbeamcyl.h"

/* 
   Projet BAORadio / HSHS
   Programme de simulation pour reconstruction de lobe radio.
   programme principal de test

   R. Ansari - LAL      Jan 2007

*/

// Declaration des fonctions de ce fichier
static int test1cyl(string& ppfname);
static int testmulticyl(string& ppfname, int ncyl=5);

//-----------------------------------------------------------
// -------------- Parametres de simulation  -----------------
//-----------------------------------------------------------
static int MR = 256;  // Nombre de recepteur
static int NE = 1024;  // Nombre d'echantillon en temps;
static double freq0 = 2.;  // frequence de base
static double da = 0.25;     // pas des antennes le long du cylindre
// ATTENTION : les parametres suivants sont relies a MR/da
static double maxangX = M_PI/3.; // angle max en X ( +/- )
static double maxangY = M_PI/60.; // angle max en Y ( +/- )
static int nsrcmax = 50;  // Nb total de sources - en un plan

static double snoise = 1.0;  // sigma du bruit 
static double tjit = 0.05;   // sigma du jitter en temps
static double tos = 0.02;    // sigma des offsets en temps
static double gmean = 1.;    // gain moyen
static double gsig = 0.;     // sigma des gains
static int nantgz = 0;       // nb d'antennes morts (-> gain=0)
static int prtlevel = 0;     // niveau de print 
//-----------------------------------------------------------


/* --------------------------------------------------------
  Le main programme de test des classes de reconstruction
  multilobe radio - R. Ansari , Sep06 -- 2007
  --------------------------------------------------------- */
   
int main(int narg, char* arg[])
{

SophyaInit();
InitTim();   // Initializing the CPU timer

string ppfname = "treccyl.ppf";
int act = 1;
int ncyl = 5; 
if (narg < 3) {
  cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
       << " -act= X ou XY5 ou XY12  (5 ou 12 cylindres) \n"
       << " -ppfname=  treccyl.ppf par defaut" << endl;
  return 1;
}
if (strcmp(arg[1],"XY5") == 0) { act = 2;  ncyl = 5; }
else if (strcmp(arg[1],"XY12") == 0) { act = 2;  ncyl = 12; }
if (narg > 2)  ppfname = arg[2];
if (narg > 3)  prtlevel = atoi(arg[3]);

int rc = 0;
cout << ">>>> treccyl : " << arg[1] << " PPFName=" << ppfname << endl;
try {
  if (act == 2) rc = testmulticyl(ppfname, ncyl);
  else rc = test1cyl(ppfname);
}
  catch (PThrowable& exc) {
    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n treccyl.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " treccyl.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

cout << ">>>> treccyl ------- FIN ----------- Rc=" << rc << endl;
return rc;
}


//--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
int test1cyl(string& ppfname)
{

  // BRSourceGen sg;
  int nsrc = 60; 
  BRSourceGen sg(nsrcmax, maxangX, 0.);
  //  sg.WritePPF(string("brsrc1.ppf"));

  cout << "=== test1cyl: BRSourceGen NbSrc= " << sg.NbSources()
       << " NbRecep=" << MR << " NSamples=" << NE << endl;

  // BRSourceGen  sg(string("brsrc1.ppf"));
  if (prtlevel > 1)  sg.Print(cout);
  
  
  MultiBeamCyl  mb(MR, NE);
  mb.SetPrintLevel(prtlevel);
  mb.SetBaseFreqDa(freq0, da);
  mb.SetNoiseSigma(snoise);
  mb.SetTimeJitter(tjit);
  mb.SetTimeOffsetSigma(tos);
  mb.SetGains(gmean, gsig, nantgz);

  mb.SetSources(sg);
  
  mb.ComputeTimeVectors();
  mb.ComputeSignalVector(0, true);
  cout << "treccy/test1cyl:  signal vectors OK " << endl;
  PrtTim("test1cyl:[1] ");
  
  cout << "--- treccy/test1cyl: Saving to PPF file " << ppfname << endl;

  POutPersist po(ppfname);
  po << PPFNameTag("signal") << mb.signal;
  po << PPFNameTag("sigjitt") << mb.sigjitt;
  po << PPFNameTag("f_sig") << mb.f_sig;
  po << PPFNameTag("f_sigjit") << mb.f_sigjit;

  NTuple ntsrc = sg.Convert2Table(freq0);
  po << PPFNameTag("ntsrc") << ntsrc;

  cout << "treccy/test1cyl: - sig/f_sig,ntsrc to OutPPF OK  " << endl;
  PrtTim("test1cyl[2] ");
  
  mb.ReconstructSourcePlane(true);
  {
  TMatrix<r_4> srcplane = module(mb.getRecSrcPlane() );
  po << PPFNameTag("recsrcplane") << srcplane;
  }
  PrtTim("test1cyl[3] ");

  return 0;

}


//--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
int testmulticyl(string& ppfname, int ncyl)
{

  if ((ncyl != 5) && (ncyl != 12)) ncyl = 5;

  // BRSourceGen sg;
  int nsf = 6;
  vector<double> frq;
  frq.push_back(0.1);
  frq.push_back(0.27);
  frq.push_back(0.38);
  
  
  cout << "testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf 
       << "," << maxangX << "," << maxangY << ")" << endl;
  BRSourceGen sg(frq, nsf, maxangX, maxangY);
  
  int is;
  double fay[6] = {-0.7,-0.5,0.,0.,0.5,0.7};
  for(is=0; is<3*nsf; is++) {
    int ism = is%nsf;
    sg.angX(is) = maxangX*(ism-2.5)/3.;
    sg.angY(is) = maxangY*fay[ism];
  }
  // sg.WritePPF(string("brsrcm.ppf"));
  // BRSourceGen  sg(string("brsrcm.ppf"));
  cout << "=== testmulticyl: NbSrc= " << sg.NbSources() 
       << " NbRecep=" << MR << " NSamples=" << NE << " NCyl=" << ncyl << endl;
  if (prtlevel > 1)  sg.Print(cout);
  
  
  MultiCylinders  mcyl (MR, NE);
  mcyl.SetPrintLevel(prtlevel);
  mcyl.SetBaseFreqDa(freq0, da);
  mcyl.SetNoiseSigma(snoise);
  mcyl.SetTimeJitter(tjit);
  mcyl.SetTimeOffsetSigma(tos);
  mcyl.SetGains(gmean, gsig, nantgz);

  if (ncyl == 12) {
    cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... " << (ncyl-1)*5. 
         << " (EqualSpacing)" << endl;
    for(int kkc=0; kkc<12; kkc++) {
      cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << 5.*kkc << " )" << endl;
      mcyl.AddCylinder(5.*kkc);
    }
  }
  else {
    cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... 55 (IrregularSpacing)" << endl;
    double posyc[5] = {0.,5.,15.,35.,55.};
    for(int kkc=0; kkc<5; kkc++) {
      cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << posyc[kkc] << " )" << endl;
      mcyl.AddCylinder(posyc[kkc]);
    }
  }

  mcyl.SetSources(sg);
  
  PrtTim("testmulticyl[1] ");

  //  mcyl.ReconstructCylinderPlaneS(true);
  mcyl.ReconstructSourceBox(10, maxangY/10.);

  cout << "--- treccy/testmulticyl: Saving to PPF file " << ppfname << endl;
  POutPersist po(ppfname);

  NTuple ntsrc = sg.Convert2Table(freq0);
  po << PPFNameTag("ntsrc") << ntsrc;

  {
    //  TMatrix<r_4> srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane());
  TMatrix< complex<r_4> > srcplane0 = mcyl.GetCylinder(0).getRecSrcPlane();
  po << PPFNameTag("recsrcplane0") << srcplane0;
  }

  {
    //  TMatrix<r_4> srcplane2 = module(mcyl.GetCylinder(3).getRecSrcPlane());
  TMatrix< complex<r_4> > srcplane2 = mcyl.GetCylinder(2).getRecSrcPlane();
  po << PPFNameTag("recsrcplane2") << srcplane2;
  }

  {
    //  TMatrix<r_4> srcplane3 = module(mcyl.GetCylinder(3).getRecSrcPlane());
  TMatrix< complex<r_4> > srcplane3 = mcyl.GetCylinder(0).getRecSrcPlane();
  po << PPFNameTag("recsrcplane3") << srcplane3;
  }

  PrtTim("testmulticyl[2] ");

  int kfmin, kfmax;
  po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();
  kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[0] - 2;  kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice0 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  {
  TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf0") << slice0;
  }
  kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[1] - 2;  kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice1 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  {
  TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf1") << slice1;
  }
  kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[2] - 2;  kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice2 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  {
  TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf2") << slice2;
  }

  PrtTim("testmulticyl[3] ");

  return 0;

}