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

#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"


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

//-----------------------------------------------------------
// -------------- 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
static double snoise = 0.; //1.0;   sigma du bruit 
static double tjit = 0.; //0.05;    sigma du jitter en temps
static double tos = 0.; //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;
if (narg < 3) {
  cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
       << "  act= X ou XY  , ppfname=  treccyl.ppf par defaut" << endl;
  return 1;
}
if (strcmp(arg[1],"XY") == 0) act = 2;
if (narg > 2)  ppfname = arg[2];
if (narg > 3)  prtlevel = atoi(arg[3]);

int rc = 0;
try {
  if (act == 2) rc = testmulticyl(ppfname);
  else rc = test1cyl(ppfname);
  cout << "treccy/Info - FIN  " << endl;
}
  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; 
  } 

 return rc;
}


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

  // BRSourceGen sg;
  int nsrc = 60; 
  BRSourceGen sg(nsrc);
  //  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] ");

  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)
{

  // BRSourceGen sg;
  int nsf = 7;
  vector<double> frq;
  frq.push_back(0.1);
  frq.push_back(0.27);
  frq.push_back(0.38);
  
  
  cout << "treccy/testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf << ")" << endl;
  BRSourceGen sg(frq, nsf, M_PI/3, M_PI/150);
  int is;
  for(is=0; is<nsf; is++)     sg.angY(is) = 0.;
  for(is=nsf; is<2*nsf; is++)     sg.angY(is) = M_PI/150.;
  for(is=2*nsf; is<3*nsf; is++)     sg.angY(is) = -M_PI/150.;
  // sg.WritePPF(string("brsrcm.ppf"));
  // BRSourceGen  sg(string("brsrcm.ppf"));
  cout << "=== testmulticyl: BRSourceGen NbSrc= " << sg.NbSources() 
       << " NbRecep=" << MR << " NSamples=" << NE << 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);

  mcyl.AddCylinder(0.);
  mcyl.AddCylinder(5.);
  mcyl.AddCylinder(10.);
  mcyl.AddCylinder(18.);
  mcyl.AddCylinder(30.);
  mcyl.AddCylinder(46.);
  mcyl.AddCylinder(64.);

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

  //  mcyl.ReconstructCylinderPlaneS(true);
  mcyl.ReconstructSourceBox(15, M_PI/500.);

  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()*frq[0] - 1;  kfmax = kfmin+2;
  {
  TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf0") << slice0;
  }
  kfmin = mcyl.getRecSrcBox().SizeZ()*frq[1] - 1;  kfmax = kfmin+2;
  {
  TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf1") << slice1;
  }
  kfmin = mcyl.getRecSrcBox().SizeZ()*frq[2] - 1;  kfmax = kfmin+2;
  {
  TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf2") << slice2;
  }

  PrtTim("testmulticyl[3] ");

  return 0;

}