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

Last change on this file since 3164 was 3164, checked in by ansari, 19 years ago

modifs et corrections diverses ds treccyl.cc et recbeam.pic, Reza 31/01/2007

File size: 6.8 KB
Line 
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// Declaration des fonctions de ce fichier
25static int test1cyl(string& ppfname);
26static int testmulticyl(string& ppfname);
27
28//-----------------------------------------------------------
29// -------------- Parametres de simulation -----------------
30//-----------------------------------------------------------
31static int MR = 256; // Nombre de recepteur
32static int NE = 1024; // Nombre d'echantillon en temps;
33static double freq0 = 2.; // frequence de base
34static double da = 0.25; // pas des antennes le long du cylindre
35static double snoise = 1.0; // sigma du bruit
36static double tjit = 0.05; // sigma du jitter en temps
37static double tos = 0.02; // sigma des offsets en temps
38static double gmean = 1.; // gain moyen
39static double gsig = 0.; // sigma des gains
40static int nantgz = 0; // nb d'antennes morts (-> gain=0)
41static int prtlevel = 0; // niveau de print
42//-----------------------------------------------------------
43
44
45/* --------------------------------------------------------
46 Le main programme de test des classes de reconstruction
47 multilobe radio - R. Ansari , Sep06 -- 2007
48 --------------------------------------------------------- */
49
50int main(int narg, char* arg[])
51{
52
53SophyaInit();
54InitTim(); // Initializing the CPU timer
55
56string ppfname = "treccyl.ppf";
57int act = 1;
58if (narg < 3) {
59 cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
60 << " act= X ou XY , ppfname= treccyl.ppf par defaut" << endl;
61 return 1;
62}
63if (strcmp(arg[1],"XY") == 0) act = 2;
64if (narg > 2) ppfname = arg[2];
65if (narg > 3) prtlevel = atoi(arg[3]);
66
67int rc = 0;
68try {
69 if (act == 2) rc = testmulticyl(ppfname);
70 else rc = test1cyl(ppfname);
71 cout << "treccy/Info - FIN " << endl;
72}
73 catch (PThrowable& exc) {
74 cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
75 rc = 77;
76 }
77 catch (std::exception& sex) {
78 cerr << "\n treccyl.cc std::exception :"
79 << (string)typeid(sex).name() << "\n msg= "
80 << sex.what() << endl;
81 }
82 catch (...) {
83 cerr << " treccyl.cc catched unknown (...) exception " << endl;
84 rc = 78;
85 }
86
87 return rc;
88}
89
90
91//--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
92int test1cyl(string& ppfname)
93{
94
95 // BRSourceGen sg;
96 int nsrc = 60;
97 BRSourceGen sg(nsrc);
98 // sg.WritePPF(string("brsrc1.ppf"));
99
100 cout << "=== test1cyl: BRSourceGen NbSrc= " << sg.NbSources()
101 << " NbRecep=" << MR << " NSamples=" << NE << endl;
102
103 // BRSourceGen sg(string("brsrc1.ppf"));
104 if (prtlevel > 1) sg.Print(cout);
105
106
107 MultiBeamCyl mb(MR, NE);
108 mb.SetPrintLevel(prtlevel);
109 mb.SetBaseFreqDa(freq0, da);
110 mb.SetNoiseSigma(snoise);
111 mb.SetTimeJitter(tjit);
112 mb.SetTimeOffsetSigma(tos);
113 mb.SetGains(gmean, gsig, nantgz);
114
115 mb.SetSources(sg);
116
117 mb.ComputeTimeVectors();
118 mb.ComputeSignalVector(0, true);
119 cout << "treccy/test1cyl: signal vectors OK " << endl;
120 PrtTim("test1cyl:[1] ");
121
122 POutPersist po(ppfname);
123 po << PPFNameTag("signal") << mb.signal;
124 po << PPFNameTag("sigjitt") << mb.sigjitt;
125 po << PPFNameTag("f_sig") << mb.f_sig;
126 po << PPFNameTag("f_sigjit") << mb.f_sigjit;
127
128 NTuple ntsrc = sg.Convert2Table(freq0);
129 po << PPFNameTag("ntsrc") << ntsrc;
130
131 cout << "treccy/test1cyl: - sig/f_sig,ntsrc to OutPPF OK " << endl;
132 PrtTim("test1cyl[2] ");
133
134 mb.ReconstructSourcePlane(true);
135 {
136 TMatrix<r_4> srcplane = module(mb.getRecSrcPlane() );
137 po << PPFNameTag("recsrcplane") << srcplane;
138 }
139 PrtTim("test1cyl[3] ");
140
141 return 0;
142
143}
144
145
146//--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
147int testmulticyl(string& ppfname)
148{
149
150 // BRSourceGen sg;
151 int nsf = 6;
152 vector<double> frq;
153 frq.push_back(0.1);
154 frq.push_back(0.27);
155 frq.push_back(0.38);
156
157
158 cout << "treccy/testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf << ")" << endl;
159 BRSourceGen sg(frq, nsf, M_PI/3, M_PI/60);
160 /*
161 int is;
162 for(is=0; is<nsf; is++) sg.angY(is) = 0.;
163 for(is=nsf; is<2*nsf; is++) sg.angY(is) = M_PI/30.;
164 for(is=2*nsf; is<3*nsf; is++) sg.angY(is) = -M_PI/30.;
165 */
166 // sg.WritePPF(string("brsrcm.ppf"));
167 // BRSourceGen sg(string("brsrcm.ppf"));
168 cout << "=== testmulticyl: BRSourceGen NbSrc= " << sg.NbSources()
169 << " NbRecep=" << MR << " NSamples=" << NE << endl;
170 if (prtlevel > 1) sg.Print(cout);
171
172
173 MultiCylinders mcyl (MR, NE);
174 mcyl.SetPrintLevel(prtlevel);
175 mcyl.SetBaseFreqDa(freq0, da);
176 mcyl.SetNoiseSigma(snoise);
177 mcyl.SetTimeJitter(tjit);
178 mcyl.SetTimeOffsetSigma(tos);
179 mcyl.SetGains(gmean, gsig, nantgz);
180
181 /*
182 for(int kkc=0; kkc<12; kkc++)
183 mcyl.AddCylinder(5.*kkc);
184 */
185
186 mcyl.AddCylinder(0.);
187 mcyl.AddCylinder(5.);
188 mcyl.AddCylinder(15.);
189 mcyl.AddCylinder(35.);
190 mcyl.AddCylinder(65.);
191
192 mcyl.SetSources(sg);
193
194 PrtTim("testmulticyl[1] ");
195
196 // mcyl.ReconstructCylinderPlaneS(true);
197 mcyl.ReconstructSourceBox(10, M_PI/600.);
198
199 POutPersist po(ppfname);
200
201 NTuple ntsrc = sg.Convert2Table(freq0);
202 po << PPFNameTag("ntsrc") << ntsrc;
203
204 {
205 // TMatrix<r_4> srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane());
206 TMatrix< complex<r_4> > srcplane0 = mcyl.GetCylinder(0).getRecSrcPlane();
207 po << PPFNameTag("recsrcplane0") << srcplane0;
208 }
209
210 {
211 // TMatrix<r_4> srcplane2 = module(mcyl.GetCylinder(3).getRecSrcPlane());
212 TMatrix< complex<r_4> > srcplane2 = mcyl.GetCylinder(2).getRecSrcPlane();
213 po << PPFNameTag("recsrcplane2") << srcplane2;
214 }
215
216 {
217 // TMatrix<r_4> srcplane3 = module(mcyl.GetCylinder(3).getRecSrcPlane());
218 TMatrix< complex<r_4> > srcplane3 = mcyl.GetCylinder(0).getRecSrcPlane();
219 po << PPFNameTag("recsrcplane3") << srcplane3;
220 }
221
222 PrtTim("testmulticyl[2] ");
223
224 int kfmin, kfmax;
225 po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();
226 kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[0] - 2; kfmax = kfmin+2;
227 cout << "testmulticyl/Info: slice0 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
228 {
229 TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
230 po << PPFNameTag("recXYf0") << slice0;
231 }
232 kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[1] - 2; kfmax = kfmin+2;
233 cout << "testmulticyl/Info: slice1 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
234 {
235 TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
236 po << PPFNameTag("recXYf1") << slice1;
237 }
238 kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[2] - 2; kfmax = kfmin+2;
239 cout << "testmulticyl/Info: slice2 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
240 {
241 TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
242 po << PPFNameTag("recXYf2") << slice2;
243 }
244
245 PrtTim("testmulticyl[3] ");
246
247 return 0;
248
249}
Note: See TracBrowser for help on using the repository browser.