| 1 | #include "sopnamsp.h"
 | 
|---|
| 2 | #include "machdefs.h"
 | 
|---|
| 3 | #include <math.h>
 | 
|---|
| 4 | #include <iostream>
 | 
|---|
| 5 | #include <typeinfo>
 | 
|---|
| 6 | 
 | 
|---|
| 7 | #include "tvector.h"
 | 
|---|
| 8 | #include "srandgen.h"
 | 
|---|
| 9 | #include "fioarr.h"
 | 
|---|
| 10 | #include "sopemtx.h"
 | 
|---|
| 11 | #include "pexceptions.h"
 | 
|---|
| 12 | #include "matharr.h"
 | 
|---|
| 13 | 
 | 
|---|
| 14 | #include "sambainit.h"
 | 
|---|
| 15 | 
 | 
|---|
| 16 | // #include "tarrinit.h"
 | 
|---|
| 17 | 
 | 
|---|
| 18 | #include "timing.h"
 | 
|---|
| 19 | 
 | 
|---|
| 20 | #include "multicyl.h"
 | 
|---|
| 21 | #include "mbeamcyl.h"
 | 
|---|
| 22 | 
 | 
|---|
| 23 | /* 
 | 
|---|
| 24 |    Projet BAORadio / HSHS
 | 
|---|
| 25 |    Programme de simulation pour reconstruction de lobe radio.
 | 
|---|
| 26 |    programme principal de test
 | 
|---|
| 27 | 
 | 
|---|
| 28 |    R. Ansari - LAL      Jan 2007
 | 
|---|
| 29 | 
 | 
|---|
| 30 | */
 | 
|---|
| 31 | 
 | 
|---|
| 32 | // Declaration des fonctions de ce fichier
 | 
|---|
| 33 | static int test1cyl(string& ppfname);
 | 
|---|
| 34 | static int testmulticyl(string& ppfname, int ncyl=5);
 | 
|---|
| 35 | 
 | 
|---|
| 36 | //-----------------------------------------------------------
 | 
|---|
| 37 | // -------------- Parametres de simulation  -----------------
 | 
|---|
| 38 | //-----------------------------------------------------------
 | 
|---|
| 39 | static int MR = 256;  // Nombre de recepteur
 | 
|---|
| 40 | static int NE = 1024;  // Nombre d'echantillon en temps;
 | 
|---|
| 41 | static double freq0 = 2.;  // frequence de base
 | 
|---|
| 42 | static double da = 0.25;     // pas des antennes le long du cylindre
 | 
|---|
| 43 | // ATTENTION : les parametres suivants sont relies a MR/da
 | 
|---|
| 44 | static double maxangX = M_PI/3.; // angle max en X ( +/- )
 | 
|---|
| 45 | static double maxangY = M_PI/60.; // angle max en Y ( +/- )
 | 
|---|
| 46 | static int nsrcmax = 50;  // Nb total de sources - en un plan
 | 
|---|
| 47 | 
 | 
|---|
| 48 | static double snoise = 1.0;  // sigma du bruit 
 | 
|---|
| 49 | static double tjit = 0.05;   // sigma du jitter en temps
 | 
|---|
| 50 | static double tos = 0.02;    // sigma des offsets en temps
 | 
|---|
| 51 | static double gmean = 1.;    // gain moyen
 | 
|---|
| 52 | static double gsig = 0.;     // sigma des gains
 | 
|---|
| 53 | static int nantgz = 0;       // nb d'antennes morts (-> gain=0)
 | 
|---|
| 54 | static int prtlevel = 0;     // niveau de print 
 | 
|---|
| 55 | //-----------------------------------------------------------
 | 
|---|
| 56 | 
 | 
|---|
| 57 | 
 | 
|---|
| 58 | /* --------------------------------------------------------
 | 
|---|
| 59 |   Le main programme de test des classes de reconstruction
 | 
|---|
| 60 |   multilobe radio - R. Ansari , Sep06 -- 2007
 | 
|---|
| 61 |   --------------------------------------------------------- */
 | 
|---|
| 62 |    
 | 
|---|
| 63 | int main(int narg, char* arg[])
 | 
|---|
| 64 | {
 | 
|---|
| 65 | 
 | 
|---|
| 66 | SophyaInit();
 | 
|---|
| 67 | InitTim();   // Initializing the CPU timer
 | 
|---|
| 68 | 
 | 
|---|
| 69 | string ppfname = "treccyl.ppf";
 | 
|---|
| 70 | int act = 1;
 | 
|---|
| 71 | int ncyl = 5; 
 | 
|---|
| 72 | if (narg < 3) {
 | 
|---|
| 73 |   cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
 | 
|---|
| 74 |        << " -act= X ou XY5 ou XY12  (5 ou 12 cylindres) \n"
 | 
|---|
| 75 |        << " -ppfname=  treccyl.ppf par defaut" << endl;
 | 
|---|
| 76 |   return 1;
 | 
|---|
| 77 | }
 | 
|---|
| 78 | if (strcmp(arg[1],"XY5") == 0) { act = 2;  ncyl = 5; }
 | 
|---|
| 79 | else if (strcmp(arg[1],"XY12") == 0) { act = 2;  ncyl = 12; }
 | 
|---|
| 80 | if (narg > 2)  ppfname = arg[2];
 | 
|---|
| 81 | if (narg > 3)  prtlevel = atoi(arg[3]);
 | 
|---|
| 82 | 
 | 
|---|
| 83 | int rc = 0;
 | 
|---|
| 84 | cout << ">>>> treccyl : " << arg[1] << " PPFName=" << ppfname << endl;
 | 
|---|
| 85 | try {
 | 
|---|
| 86 |   if (act == 2) rc = testmulticyl(ppfname, ncyl);
 | 
|---|
| 87 |   else rc = test1cyl(ppfname);
 | 
|---|
| 88 | }
 | 
|---|
| 89 |   catch (PThrowable& exc) {
 | 
|---|
| 90 |     cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
 | 
|---|
| 91 |     rc = 77;
 | 
|---|
| 92 |   }  
 | 
|---|
| 93 |   catch (std::exception& sex) {
 | 
|---|
| 94 |     cerr << "\n treccyl.cc std::exception :" 
 | 
|---|
| 95 |          << (string)typeid(sex).name() << "\n msg= " 
 | 
|---|
| 96 |          << sex.what() << endl;
 | 
|---|
| 97 |   }
 | 
|---|
| 98 |   catch (...) {
 | 
|---|
| 99 |     cerr << " treccyl.cc catched unknown (...) exception  " << endl; 
 | 
|---|
| 100 |     rc = 78; 
 | 
|---|
| 101 |   } 
 | 
|---|
| 102 | 
 | 
|---|
| 103 | cout << ">>>> treccyl ------- FIN ----------- Rc=" << rc << endl;
 | 
|---|
| 104 | return rc;
 | 
|---|
| 105 | }
 | 
|---|
| 106 | 
 | 
|---|
| 107 | 
 | 
|---|
| 108 | //--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
 | 
|---|
| 109 | int test1cyl(string& ppfname)
 | 
|---|
| 110 | {
 | 
|---|
| 111 | 
 | 
|---|
| 112 |   // BRSourceGen sg;
 | 
|---|
| 113 |   int nsrc = 60; 
 | 
|---|
| 114 |   BRSourceGen sg(nsrcmax, maxangX, 0.);
 | 
|---|
| 115 |   //  sg.WritePPF(string("brsrc1.ppf"));
 | 
|---|
| 116 | 
 | 
|---|
| 117 |   cout << "=== test1cyl: BRSourceGen NbSrc= " << sg.NbSources()
 | 
|---|
| 118 |        << " NbRecep=" << MR << " NSamples=" << NE << endl;
 | 
|---|
| 119 | 
 | 
|---|
| 120 |   // BRSourceGen  sg(string("brsrc1.ppf"));
 | 
|---|
| 121 |   if (prtlevel > 1)  sg.Print(cout);
 | 
|---|
| 122 |   
 | 
|---|
| 123 |   
 | 
|---|
| 124 |   MultiBeamCyl  mb(MR, NE);
 | 
|---|
| 125 |   mb.SetPrintLevel(prtlevel);
 | 
|---|
| 126 |   mb.SetBaseFreqDa(freq0, da);
 | 
|---|
| 127 |   mb.SetNoiseSigma(snoise);
 | 
|---|
| 128 |   mb.SetTimeJitter(tjit);
 | 
|---|
| 129 |   mb.SetTimeOffsetSigma(tos);
 | 
|---|
| 130 |   mb.SetGains(gmean, gsig, nantgz);
 | 
|---|
| 131 | 
 | 
|---|
| 132 |   mb.SetSources(sg);
 | 
|---|
| 133 |   
 | 
|---|
| 134 |   mb.ComputeTimeVectors();
 | 
|---|
| 135 |   mb.ComputeSignalVector(0, true);
 | 
|---|
| 136 |   cout << "treccy/test1cyl:  signal vectors OK " << endl;
 | 
|---|
| 137 |   PrtTim("test1cyl:[1] ");
 | 
|---|
| 138 |   
 | 
|---|
| 139 |   cout << "--- treccy/test1cyl: Saving to PPF file " << ppfname << endl;
 | 
|---|
| 140 | 
 | 
|---|
| 141 |   POutPersist po(ppfname);
 | 
|---|
| 142 |   po << PPFNameTag("signal") << mb.signal;
 | 
|---|
| 143 |   po << PPFNameTag("sigjitt") << mb.sigjitt;
 | 
|---|
| 144 |   po << PPFNameTag("f_sig") << mb.f_sig;
 | 
|---|
| 145 |   po << PPFNameTag("f_sigjit") << mb.f_sigjit;
 | 
|---|
| 146 | 
 | 
|---|
| 147 |   NTuple ntsrc = sg.Convert2Table(freq0);
 | 
|---|
| 148 |   po << PPFNameTag("ntsrc") << ntsrc;
 | 
|---|
| 149 | 
 | 
|---|
| 150 |   cout << "treccy/test1cyl: - sig/f_sig,ntsrc to OutPPF OK  " << endl;
 | 
|---|
| 151 |   PrtTim("test1cyl[2] ");
 | 
|---|
| 152 |   
 | 
|---|
| 153 |   mb.ReconstructSourcePlane(true);
 | 
|---|
| 154 |   {
 | 
|---|
| 155 |   TMatrix<r_4> srcplane = module(mb.getRecSrcPlane() );
 | 
|---|
| 156 |   po << PPFNameTag("recsrcplane") << srcplane;
 | 
|---|
| 157 |   }
 | 
|---|
| 158 |   PrtTim("test1cyl[3] ");
 | 
|---|
| 159 | 
 | 
|---|
| 160 |   return 0;
 | 
|---|
| 161 | 
 | 
|---|
| 162 | }
 | 
|---|
| 163 | 
 | 
|---|
| 164 | 
 | 
|---|
| 165 | //--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
 | 
|---|
| 166 | int testmulticyl(string& ppfname, int ncyl)
 | 
|---|
| 167 | {
 | 
|---|
| 168 | 
 | 
|---|
| 169 |   if ((ncyl != 5) && (ncyl != 12)) ncyl = 5;
 | 
|---|
| 170 | 
 | 
|---|
| 171 |   // BRSourceGen sg;
 | 
|---|
| 172 |   int nsf = 6;
 | 
|---|
| 173 |   vector<double> frq;
 | 
|---|
| 174 |   frq.push_back(0.1);
 | 
|---|
| 175 |   frq.push_back(0.27);
 | 
|---|
| 176 |   frq.push_back(0.38);
 | 
|---|
| 177 |   
 | 
|---|
| 178 |   
 | 
|---|
| 179 |   cout << "testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf 
 | 
|---|
| 180 |        << "," << maxangX << "," << maxangY << ")" << endl;
 | 
|---|
| 181 |   BRSourceGen sg(frq, nsf, maxangX, maxangY);
 | 
|---|
| 182 |   
 | 
|---|
| 183 |   int is;
 | 
|---|
| 184 |   double fay[6] = {-0.7,-0.5,0.,0.,0.5,0.7};
 | 
|---|
| 185 |   for(is=0; is<3*nsf; is++) {
 | 
|---|
| 186 |     int ism = is%nsf;
 | 
|---|
| 187 |     sg.angX(is) = maxangX*(ism-2.5)/3.;
 | 
|---|
| 188 |     sg.angY(is) = maxangY*fay[ism];
 | 
|---|
| 189 |   }
 | 
|---|
| 190 |   // sg.WritePPF(string("brsrcm.ppf"));
 | 
|---|
| 191 |   // BRSourceGen  sg(string("brsrcm.ppf"));
 | 
|---|
| 192 |   cout << "=== testmulticyl: NbSrc= " << sg.NbSources() 
 | 
|---|
| 193 |        << " NbRecep=" << MR << " NSamples=" << NE << " NCyl=" << ncyl << endl;
 | 
|---|
| 194 |   if (prtlevel > 1)  sg.Print(cout);
 | 
|---|
| 195 |   
 | 
|---|
| 196 |   
 | 
|---|
| 197 |   MultiCylinders  mcyl (MR, NE);
 | 
|---|
| 198 |   mcyl.SetPrintLevel(prtlevel);
 | 
|---|
| 199 |   mcyl.SetBaseFreqDa(freq0, da);
 | 
|---|
| 200 |   mcyl.SetNoiseSigma(snoise);
 | 
|---|
| 201 |   mcyl.SetTimeJitter(tjit);
 | 
|---|
| 202 |   mcyl.SetTimeOffsetSigma(tos);
 | 
|---|
| 203 |   mcyl.SetGains(gmean, gsig, nantgz);
 | 
|---|
| 204 | 
 | 
|---|
| 205 |   if (ncyl == 12) {
 | 
|---|
| 206 |     cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... " << (ncyl-1)*5. 
 | 
|---|
| 207 |          << " (EqualSpacing)" << endl;
 | 
|---|
| 208 |     for(int kkc=0; kkc<12; kkc++) {
 | 
|---|
| 209 |       cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << 5.*kkc << " )" << endl;
 | 
|---|
| 210 |       mcyl.AddCylinder(5.*kkc);
 | 
|---|
| 211 |     }
 | 
|---|
| 212 |   }
 | 
|---|
| 213 |   else {
 | 
|---|
| 214 |     cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... 55 (IrregularSpacing)" << endl;
 | 
|---|
| 215 |     double posyc[5] = {0.,5.,15.,35.,55.};
 | 
|---|
| 216 |     for(int kkc=0; kkc<5; kkc++) {
 | 
|---|
| 217 |       cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << posyc[kkc] << " )" << endl;
 | 
|---|
| 218 |       mcyl.AddCylinder(posyc[kkc]);
 | 
|---|
| 219 |     }
 | 
|---|
| 220 |   }
 | 
|---|
| 221 | 
 | 
|---|
| 222 |   mcyl.SetSources(sg);
 | 
|---|
| 223 |   
 | 
|---|
| 224 |   PrtTim("testmulticyl[1] ");
 | 
|---|
| 225 | 
 | 
|---|
| 226 |   //  mcyl.ReconstructCylinderPlaneS(true);
 | 
|---|
| 227 |   mcyl.ReconstructSourceBox(10, maxangY/10.);
 | 
|---|
| 228 | 
 | 
|---|
| 229 |   cout << "--- treccy/testmulticyl: Saving to PPF file " << ppfname << endl;
 | 
|---|
| 230 |   POutPersist po(ppfname);
 | 
|---|
| 231 | 
 | 
|---|
| 232 |   NTuple ntsrc = sg.Convert2Table(freq0);
 | 
|---|
| 233 |   po << PPFNameTag("ntsrc") << ntsrc;
 | 
|---|
| 234 | 
 | 
|---|
| 235 |   {
 | 
|---|
| 236 |     //  TMatrix<r_4> srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane());
 | 
|---|
| 237 |   TMatrix< complex<r_4> > srcplane0 = mcyl.GetCylinder(0).getRecSrcPlane();
 | 
|---|
| 238 |   po << PPFNameTag("recsrcplane0") << srcplane0;
 | 
|---|
| 239 |   }
 | 
|---|
| 240 | 
 | 
|---|
| 241 |   {
 | 
|---|
| 242 |     //  TMatrix<r_4> srcplane2 = module(mcyl.GetCylinder(3).getRecSrcPlane());
 | 
|---|
| 243 |   TMatrix< complex<r_4> > srcplane2 = mcyl.GetCylinder(2).getRecSrcPlane();
 | 
|---|
| 244 |   po << PPFNameTag("recsrcplane2") << srcplane2;
 | 
|---|
| 245 |   }
 | 
|---|
| 246 | 
 | 
|---|
| 247 |   {
 | 
|---|
| 248 |     //  TMatrix<r_4> srcplane3 = module(mcyl.GetCylinder(3).getRecSrcPlane());
 | 
|---|
| 249 |   TMatrix< complex<r_4> > srcplane3 = mcyl.GetCylinder(0).getRecSrcPlane();
 | 
|---|
| 250 |   po << PPFNameTag("recsrcplane3") << srcplane3;
 | 
|---|
| 251 |   }
 | 
|---|
| 252 | 
 | 
|---|
| 253 |   PrtTim("testmulticyl[2] ");
 | 
|---|
| 254 | 
 | 
|---|
| 255 |   int kfmin, kfmax;
 | 
|---|
| 256 |   po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();
 | 
|---|
| 257 |   kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[0] - 2;  kfmax = kfmin+2;
 | 
|---|
| 258 |   cout << "testmulticyl/Info: slice0 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
 | 
|---|
| 259 |   {
 | 
|---|
| 260 |   TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
 | 
|---|
| 261 |   po << PPFNameTag("recXYf0") << slice0;
 | 
|---|
| 262 |   }
 | 
|---|
| 263 |   kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[1] - 2;  kfmax = kfmin+2;
 | 
|---|
| 264 |   cout << "testmulticyl/Info: slice1 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
 | 
|---|
| 265 |   {
 | 
|---|
| 266 |   TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
 | 
|---|
| 267 |   po << PPFNameTag("recXYf1") << slice1;
 | 
|---|
| 268 |   }
 | 
|---|
| 269 |   kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[2] - 2;  kfmax = kfmin+2;
 | 
|---|
| 270 |   cout << "testmulticyl/Info: slice2 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
 | 
|---|
| 271 |   {
 | 
|---|
| 272 |   TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
 | 
|---|
| 273 |   po << PPFNameTag("recXYf2") << slice2;
 | 
|---|
| 274 |   }
 | 
|---|
| 275 | 
 | 
|---|
| 276 |   PrtTim("testmulticyl[3] ");
 | 
|---|
| 277 | 
 | 
|---|
| 278 |   return 0;
 | 
|---|
| 279 | 
 | 
|---|
| 280 | }
 | 
|---|