source: Sophya/trunk/SophyaPI/PIext/piacmd.cc@ 171

Last change on this file since 171 was 171, checked in by ercodmgr, 27 years ago

Adaptation a PDynLinkMgr Reza 21/12/98

File size: 37.0 KB
RevLine 
[165]1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4
5#include <iostream.h>
6#include <string>
7#include <vector>
8#include <map>
9#if defined(__KCC__)
10using std::string ;
11#include <vector.h>
12#include <map.h>
13#endif
14
15#include "ctimer.h"
16
17#include "nobjmgr.h"
18#include "piacmd.h"
19#include "pistdimgapp.h"
20
21#include "histos.h"
22#include "histos2.h"
23#include "hisprof.h"
24#include "ntuple.h"
25#include "generaldata.h"
26
27
28// Pour le decoupage des commandes en lignes
29typedef vector<string> cmdtok;
30
31// Fonction pour faire le Fit de PSF sur image
32static int FitPSFImg(cmdtok& tokens);
33static int Photom_Ouv(cmdtok& tokens);
34
35
36/* --Methode-- */
37PIACmd::PIACmd(NamedObjMgr* omg, PIStdImgApp* app)
38{
39if (omg) { mOmg= false; mObjMgr = omg; }
40else { mObjMgr = new NamedObjMgr; mOmg = true; }
41mImgApp = app;
42system("cp history.pic hisold.pic");
43hist.open("history.pic");
44trace = false; timing = false;
45gltimer = NULL;
[171]46dynlink = NULL;
[165]47}
48
49/* --Methode-- */
50PIACmd::~PIACmd()
51{
52if (mOmg) delete mObjMgr;
53hist.close();
54if (gltimer) { delete gltimer; gltimer = NULL; }
55}
56
57/* --Methode-- */
58int PIACmd::Exec(string& file)
59{
60char line_buff[512];
61FILE *fip;
62
63if ( (fip = fopen(file.c_str(),"r")) == NULL ) {
64 cerr << "PIACmd::Exec() Error opening file " << file << endl;
65 hist << "##! PIACmd::Exec() Error opening file " << file << endl;
66 return(0);
67 }
68
69hist << "### Executing commands from " << file << endl;
70if (trace) {
71 mImgApp->GetConsole()->AddStr("### Executing commands from ", PIVA_Magenta);
72 mImgApp->GetConsole()->AddStr(file.c_str(), PIVA_Magenta);
73 mImgApp->GetConsole()->AddStr("\n", PIVA_Magenta);
74 }
75
76while (fgets(line_buff,511,fip) != NULL)
77 {
78 if (trace) mImgApp->GetConsole()->AddStr(line_buff, PIVA_Magenta);
79 line_buff[strlen(line_buff)-1] = '\0'; /* LF/CR de la fin */
80 string line(line_buff);
81 Do(line);
82 }
83hist << "### End of Exec( " << file << " ) " << endl;
84if (trace) {
85 mImgApp->GetConsole()->AddStr("### End of Exec( ", PIVA_Magenta);
86 mImgApp->GetConsole()->AddStr(file.c_str(), PIVA_Magenta);
87 mImgApp->GetConsole()->AddStr(" ) \n", PIVA_Magenta);
88 }
89
90return(0);
91}
92
93/* --Methode-- */
94int PIACmd::Do(string& s)
95{
96
97cmdtok tokens;
98if (s.length() < 1) return(0);
99
100hist << s << endl; // On enregistre les commandes
101
102if (s[0] == '#') {
103 cout << "PIACmd::Do() Comment-Line:" << s << endl;
104 return(0);
105 }
106string toks,kw;
107size_t p = s.find_first_not_of(" ");
108s = s.substr(p);
109p = 0;
110size_t q = s.find_first_of(" ");
111size_t l = s.length();
112
113if (q < l)
114 { kw = s.substr(p,q-p); toks = s.substr(q, l-q); }
115else { kw = s.substr(p,l-p); toks = ""; }
116
117q = 0;
118while (q < l) {
119 p = toks.find_first_not_of(" ",q+1); // au debut d'un token
120 if (p>=l) break;
121 q = toks.find_first_of(" ",p); // la fin du token;
122 string token = toks.substr(p,q-p);
123 tokens.push_back(token);
124 }
125
126for(int k=0; k<tokens.size(); k++) { // On remplace les $varname par la valeur de la variable
127 if ((tokens[k])[0] != '$') continue;
128 CmdVarList::iterator it = mVars.find(tokens[k].substr(1));
129 if (it != mVars.end()) tokens[k] = (*it).second;
130}
131
132// cout << "PIACmd::Do() DBG KeyW= " << kw << " NbArgs= " << tokens.size() << endl;
133// for(int ii=0; ii<tokens.size(); ii++)
134// cout << "arg[ " << ii << " ] : " << tokens[ii] << endl;
135
136// >>>>>>>>>>> Commande d'interpreteur
137if (kw == "help") {
138 cout << "\n PIACmd::Do() ----- Help : Liste des commandes ------------ \n";
139 cout << "help shelp exec set listvars traceon traceoff \n";
140 cout << "timingon timingoff shell \n";
141 cout << "link call zone gratt stacknext \n" ;
142 cout << "openfits openppf saveall listobjs \n";
143 cout << "print rename del delobjs \n";
144 cout << "newh1d newh2d newprof newgfd func func2d \n" ;
145 cout << "disp surf nt2d nt3d gfd2d gfd3d \n" ;
146 cout << "plot2d plot3d projh1d projh2d projprof \n";
147 cout << "fillvec fillnt fillgd1 fillgd2 fft fit fitpsf aper\n";
148 cout << "*** Commande sans args -> liste d'arguments\n";
149 cout << "*** dopt: DisplayOption = next/same/win/stack,color,... \n";
150 cout << "*** $varname remplace par la valeur de la variable varname\n" << endl;
151 }
152
153else if (kw == "shelp" ) {
154 cout << "\n PIACmd::Do() ------------- SHelp ------------ " << endl;
155 cout << "exec filename -- set varname string -- clr varname -- listvars \n";
156 cout << "traceon/off -- timingon/off -- shell cmdline\n";
157 cout << "link fnameso f1 [f2 f3 f4 f5] -- call userf [arg1 arg2 ...]\n";
158 cout << "zone nx ny -- gratt attributes_list -- stacknext \n";
159 cout << "openfits file -- openppf file -- saveall file -- listobjs --\n";
160 cout << "print nomobj -- rename nomobj nomnew -- del nomobj -- delobjs pattern (*,?) \n";
161 cout << "newh1d nomh xmin xmax nbin -- newh2d nomh xmin xmax nbinx ymin ymax nbiny\n";
162 cout << "newprof nomh xmin xmax nbin [ymin ymax] \n";
163 cout << "newgfd nomgfd nvar nalloc [errx(0/1)]\n";
164 cout << "func expfunc xmin xmax npt opt\n";
165 cout << "func2d expfunc xmin xmax nptx ymin ymax npty opt\n";
166 cout << "disp nomobj opt -- surf nomobj opt\n";
167 cout << "nt2d nomobj varx vary errx erry opt\n";
168 cout << "nt3d nomobj varx vary varz errx erry errz dopt\n";
169 cout << "gfd2d nomobj numvarx erreur=(x y xy) opt\n";
170 cout << "gfd3d nomobj numvarx numvary erreur=(x y z xy xz yz xyz) opt\n";
171 cout << "plot2d nomobj expx expy [experrx experry] expcut dopt\n";
172 cout << "plot3d nomobj expx expy expz expcut dopt\n";
173 cout << "projh1d nomobj expx expwt expcut nomh1 dopt\n";
174 cout << "projh2d nomobj expx expy expwt expcut nomh2 dopt\n";
175 cout << "projprof nomobj expx expy expwt expcut nomprof dopt\n";
176 cout << "fillnt nomobj expx expy expz expt expcut nomnt dopt\n";
177 cout << "fillvec nomobj expx expcut nomvec dopt\n";
178 cout << "fillgd1 nomobj expx expy experry expcut nomgfd\n";
179 cout << "fillgd2 nomobj expx expy expz experrz expcut nomgfd\n";
180 cout << "fft nomobj dopt\n";
181 cout << "fit nomobj func [p:p1,.,pn s:s1,.,sn m:m1,.,mn M:M1,.,Mn I:i1,i2 o:... o:...]\n";
182 cout << "fitpsf psf_typ [options_fit]\n";
183 cout << "aper fond_typ [options_calcul]\n";
184}
185
186else if (kw == "set") {
187 if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: set varname string" << endl; return(0); }
188 mVars[tokens[0]] = tokens[1];
189 }
190else if (kw == "clr") {
191 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: clr varname" << endl; return(0); }
192 CmdVarList::iterator it = mVars.find(tokens[0]);
193 if (it != mVars.end()) mVars.erase(it);
194 else cerr << "PIACmd::Do() Pas de variable de nom " << tokens[0] << endl;
195 }
196else if (kw == "listvars") {
197 cout << "PIACmd::Do() Variable List , VarName = Value \n";
198 CmdVarList::iterator it;
199 for(it = mVars.begin(); it != mVars.end(); it++)
200 cout << (*it).first << " = " << (*it).second << "\n";
201 cout << endl;
202 }
203else if (kw == "traceon") { cout << "PIACmd::Do() -> Trace ON mode " << endl; trace = true; }
204else if (kw == "traceoff") { cout << "PIACmd::Do() -> Trace OFF mode " << endl; trace = false; }
205else if (kw == "timingon") {
206 cout << "PIACmd::Do() -> Timing ON mode " << endl;
207 if (gltimer) delete gltimer; gltimer = new Timer("PIA-CmdInterpreter "); timing = true; }
208else if (kw == "timingoff") {
209 cout << "PIACmd::Do() -> Timing OFF mode " << endl;
210 if (gltimer) delete gltimer; gltimer = NULL; timing = false; }
211else if (kw == "exec") {
212 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: exec filename" << endl; return(0); }
213 Exec(tokens[0]);
214 }
215else if (kw == "shell") {
216 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: shell cmdline" << endl; return(0); }
217 system(toks.c_str());
218 }
219
220// >>>>>>>>>>> Fenetre graphique , changement d'attributs graphiques
221else if (kw == "zone") {
222 if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: zone nx ny" << endl; return(0); }
223 int nx, ny;
224 nx = ny = 1;
225 nx = atoi(tokens[0].c_str()); ny = atoi(tokens[1].c_str());
226 mObjMgr->SetGraphicWinZone(nx, ny);
227 }
228else if (kw == "gratt") {
229 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: gratt attributes_list (att=def->defaut)" << endl; return(0); }
230 mObjMgr->SetGraphicAttributes(tokens[0]);
231 }
232else if (kw == "stacknext") mImgApp->StackWinNext();
233
234// >>>>>>>>>>> Link dynamique de fonctions C++
235else if (kw == "link" ) {
236 if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: link fnameso f1 [f2 f3 f4 f5]" << endl; return(0); }
237 string sph = "";
238 for(int gg=0; gg<4; gg++) tokens.push_back(sph);
239 int rc = LinkUserFuncs(tokens[0], tokens[1], tokens[2], tokens[3], tokens[4], tokens[5]);
240 if (rc == 0) cout << "PIAcmd::Do(): Link from " << tokens[0] << " OK " << endl;
241}
242else if (kw == "call" ) {
243 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: call userf [arg1 arg2 ...]" << endl; return(0); }
244 UsFmap::iterator it = usfmap.find(tokens[0]);
245 if (it == usfmap.end()) {
246 cerr << "PIAcmd::Do() No User Function " << tokens[0] << endl;
247 return(0);
248 }
249 cout << "PIAcmd::Do() Call " << tokens[0] << "( ... )" << endl;
250// on est oblige de faire un cast etant donne qu'on
251// utilise donc des DlFunction (Reza 20/08/98) voir commentaire ds .h (pb g++)
252 DlUserProcFunction fuf = (DlUserProcFunction)(*it).second;
253 TRY {
254 tokens.erase(tokens.begin());
255 fuf(tokens);
256 } CATCH(merr) {
257 fflush(stdout);
258 cout << endl;
259 cerr << endl;
260 string es = PeidaExc(merr);
261 cerr << "PIACmd::Do() Call UserFunc Exception :" << merr << es;
262 }
263}
264
265// >>>>>>>>>>> lecture/ecriture des objets, gestion des objets
266else if (kw == "openfits" ) {
267 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: openfits file " << endl; return(0); }
268 else mObjMgr->ReadFits(tokens[0]);
269}
270else if (kw == "openppf" ) {
271 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: openppf file " << endl; return(0); }
272 mObjMgr->ReadAll(tokens[0]);
273}
274else if (kw == "saveall" ) {
275 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: saveall file " << endl; return(0); }
276 mObjMgr->SaveAll(tokens[0]);
277}
278else if (kw == "print" ) {
279 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: print nomobj " << endl; return(0); }
280 mObjMgr->PrintObj(tokens[0]);
281}
282else if (kw == "rename" ) {
283 if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: rename nomobj nomnew" << endl; return(0); }
284 mObjMgr->RenameObj(tokens[0], tokens[1]);
285}
286else if (kw == "del" ) {
287 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: del nomobj " << endl; return(0); }
288 mObjMgr->DelObj(tokens[0]);
289}
290else if (kw == "delobjs" ) {
291 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: delobjs nomobjpattern (*,?) " << endl; return(0); }
292 mObjMgr->DelObjects(tokens[0]);
293}
294else if (kw == "listobjs") mObjMgr->ListObjs();
295
296// >>>>>>>>>>> Creation d'histos 1D-2D
297else if (kw == "newh1d") {
298 if (tokens.size() < 4) { cout << "PIACmd::Do() Usage: newh1d nom xmin xmax nbin" << endl; return(0); }
299 int nbx;
300 float xmin, xmax;
301 nbx = 100;
302 xmin = 0.; xmax = 1.;
303 nbx = atoi(tokens[3].c_str());
304 xmin = atof(tokens[1].c_str()); xmax = atof(tokens[2].c_str());
305 Histo* h = new Histo(xmin, xmax, nbx);
306 mObjMgr->AddObj(h, tokens[0]);
307 }
308else if (kw == "newh2d") {
309 if (tokens.size() < 7) {
310 cout << "PIACmd::Do() Usage:newh2d nomh xmin xmax nbinx ymin ymax nbiny" << endl;
311 return(0);
312 }
313 int nbx, nby;
314 float xmin, xmax;
315 float ymin, ymax;
316 nbx = nby = 50;
317 xmin = 0.; xmax = 1.;
318 ymin = 0.; ymax = 1.;
319 nbx = atoi(tokens[3].c_str());
320 nby = atoi(tokens[6].c_str());
321 xmin = atof(tokens[1].c_str()); xmax = atof(tokens[2].c_str());
322 ymin = atof(tokens[4].c_str()); ymax = atof(tokens[5].c_str());
323 Histo2D* h = new Histo2D(xmin, xmax, nbx, ymin, ymax, nby);
324 mObjMgr->AddObj(h, tokens[0]);
325 }
326else if (kw == "newprof") {
327 if (tokens.size() < 4)
328 { cout << "PIACmd::Do() Usage: newprof nom xmin xmax nbin [ymin ymax]" << endl; return(0); }
329 int nbx;
330 float xmin, xmax, ymin = 1., ymax = -1.;
331 if(tokens.size() > 5)
332 {ymin = atof(tokens[4].c_str()); ymax = atof(tokens[5].c_str());}
333 nbx = 100;
334 xmin = 0.; xmax = 1.;
335 nbx = atoi(tokens[3].c_str());
336 xmin = atof(tokens[1].c_str()); xmax = atof(tokens[2].c_str());
337 HProf* h = new HProf(xmin, xmax, nbx, ymin, ymax);
338 mObjMgr->AddObj(h, tokens[0]);
339 }
340else if (kw == "newgfd") {
341 if (tokens.size() < 3)
342 { cout << "PIACmd::Do() Usage: newgfd nvar nalloc [errx(0/1)]" << endl; return(0); }
343 int nvar, nalloc, errx=0;
344 if (tokens.size() > 3)
345 { errx = atoi(tokens[3].c_str()); if(errx>0) errx=1; else errx = 0;}
346 nvar = atoi(tokens[1].c_str()); nalloc = atoi(tokens[2].c_str());
347 if(nvar>0 && nalloc>0) {
348 GeneralFitData* gfd = new GeneralFitData(nvar,nalloc,errx);
349 mObjMgr->AddObj(gfd, tokens[0]);
350 }
351 }
352
353// >>>>>>>>>>> Affichage des objets
354else if ( (kw == "disp") || (kw == "surf") ) {
355 if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: disp nomobj opt" << endl; return(0); }
356 string opt = "n";
357 if (tokens.size() > 1) opt = tokens[1];
358 if (kw == "disp") mObjMgr->DisplayObj(tokens[0], opt);
359 else if (kw == "surf") mObjMgr->DisplaySurf3D(tokens[0], opt);
360 }
361
362else if (kw == "nt2d") {
363 if (tokens.size() < 5) { cout << "PIACmd::Do() Usage: nt2d nomobj varx vary errx erry opt" << endl; return(0); }
364 string opt = "n";
365 if (tokens.size() > 5) opt = tokens[5];
366 string ph = "";
367 mObjMgr->DisplayNT(tokens[0],tokens[1],tokens[2], ph, tokens[3], tokens[4], ph, opt);
368 }
369else if (kw == "nt3d") {
370 if (tokens.size() < 7) { cout << "PIACmd::Do() Usage: nt3d nomobj varx vary varz errx erry errz opt" << endl; return(0); }
371 string opt = "n";
372 if (tokens.size() > 7) opt = tokens[7];
373 mObjMgr->DisplayNT(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], tokens[6], opt);
374 }
375
376else if (kw == "gfd2d") {
377 if(tokens.size()<2)
378 {cout<<"PIACmd::Do() Usage: gfd2d nomobj numvarx erreur=(x y xy) opt"<<endl;
379 return(0);}
380 string numvary = "";
381 string err = "";
382 string opt = "n";
383 if(tokens.size()>2) err = tokens[2];
384 if(tokens.size()>3) opt = tokens[3];
385 mObjMgr->DisplayGFD(tokens[0],tokens[1],numvary,err,opt);
386 }
387else if (kw == "gfd3d") {
388 if(tokens.size()<3)
389 {cout<<"PIACmd::Do() Usage: gfd3d nomobj numvarx numvary erreur=(x y z xy xz yz xyz) opt"<<endl;
390 return(0);}
391 string err = "";
392 string opt = "n";
393 if(tokens.size()>3) err = tokens[3];
394 if(tokens.size()>4) opt = tokens[4];
395 mObjMgr->DisplayGFD(tokens[0],tokens[1],tokens[2],err,opt);
396 }
397
398// >>>>>>>>>>> Trace de fonctions
399else if ( (kw == "func") ) {
400 if (tokens.size() < 4) { cout << "PIACmd::Do() Usage: func expr xmin xmax npt opt" << endl; return(0); }
401 string opt = "n";
402 if (tokens.size() > 4) opt = tokens[4];
403 int np;
404 float xmin, xmax;
405 np = 100;
406 xmin = 0.; xmax = 1.;
407 np = atoi(tokens[3].c_str());
408 xmin = atof(tokens[1].c_str()); xmax = atof(tokens[2].c_str());
409 mObjMgr->PlotFunc(tokens[0], xmin, xmax, np, opt);
410 }
411else if ( (kw == "func2d") ) {
412 if (tokens.size() < 7) {
413 cout << "PIACmd::Do() Usage: func2d expr xmin xmax nptx ymin ymax npty opt" << endl;
414 return(0);
415 }
416 int npx, npy;
417 float xmin, xmax;
418 float ymin, ymax;
419 npx = npy = 50;
420 xmin = 0.; xmax = 1.;
421 ymin = 0.; ymax = 1.;
422 npx = atoi(tokens[3].c_str());
423 npy = atoi(tokens[6].c_str());
424 xmin = atof(tokens[1].c_str()); xmax = atof(tokens[2].c_str());
425 ymin = atof(tokens[4].c_str()); ymax = atof(tokens[5].c_str());
426 string opt = "n";
427 if (tokens.size() > 7) opt = tokens[7];
428 mObjMgr->PlotFunc2D(tokens[0], xmin, xmax, ymin, ymax, npx, npy, opt);
429 }
430
431// >>>>>>>>>>> Trace d'expressions de N_Tuple, StarList, etc ...
432else if (kw == "plot2d" ) {
433 if (tokens.size() < 4) {
434 cout << "PIACmd::Do() Usage: plot2d nomobj expx expy [experrx experry] expcut opt" << endl;
435 return(0);
436 }
437 string errx = ""; string erry = ""; string ecut = "1";
438 string opt = "n";
439 if (tokens.size() < 6) { // Plot sans les erreurs
440 ecut = tokens[3];
441 if (tokens.size() > 4) opt = tokens[4];
442 }
443 else { // Plot avec les erreurs
444 errx = tokens[3]; erry = tokens[4]; ecut = tokens[5];
445 if (tokens.size() > 6) opt = tokens[6];
446 }
447 mObjMgr->DisplayPoints2D(tokens[0],tokens[1],tokens[2],errx,erry,ecut,opt);
448 }
449
450else if (kw == "plot3d" ) {
451 if (tokens.size() < 5) {
452 cout << "PIACmd::Do() Usage: plot3d nomobj expx expy expz expcut opt" << endl;
453 return(0);
454 }
455 string opt = "n";
456 if (tokens.size() > 5) opt = tokens[5];
457 mObjMgr->DisplayPoints3D(tokens[0],tokens[1],tokens[2],tokens[3], tokens[4], opt);
458 }
459
460else if (kw == "projh1d" ) {
461 if (tokens.size() < 5) {
462 cout << "PIACmd::Do() Usage: projh1d nomobj expx expwt expcut nomh1 opt" << endl;
463 return(0);
464 }
465 string opt = "n";
466 if (tokens.size() > 5) opt = tokens[5];
467 mObjMgr->ProjectH1(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], opt);
468 }
469
470else if (kw == "projh2d" ) {
471 if (tokens.size() < 6) {
472 cout << "PIACmd::Do() Usage: projh2 nomobj expx expy expwt expcut nomh2 opt" << endl;
473 return(0);
474 }
475 string opt = "n";
476 if (tokens.size() > 6) opt = tokens[6];
477
478 mObjMgr->ProjectH2(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], opt);
479 }
480
481else if (kw == "projprof" ) {
482 if (tokens.size() < 6) {
483 cout << "PIACmd::Do() Usage: projprof nomobj expx expy expwt expcut nomprof opt" << endl;
484 return(0);
485 }
486 string opt = "n";
487 if (tokens.size() > 6) opt = tokens[6];
488
489 mObjMgr->ProjectHProf(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], opt);
490 }
491
492else if (kw == "fillnt" ) {
493 if (tokens.size() < 7) {
494 cout << "PIACmd::Do() Usage: fillnt nomobj expx expy expz expt expcut nomnt" << endl;
495 return(0);
496 }
497 mObjMgr->FillNT(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], tokens[6] );
498 }
499
500else if (kw == "fillvec" ) {
501 if (tokens.size() < 4) {
502 cout << "PIACmd::Do() Usage: fillvec nomobj expx expcut nomvec opt" << endl;
503 return(0);
504 }
505 string opt = "n";
506 if (tokens.size() > 4) opt = tokens[4];
507 mObjMgr->FillVect(tokens[0],tokens[1],tokens[2], tokens[3], opt);
508 }
509
510else if (kw == "fillgd1" ) {
511 if (tokens.size() < 5) {
512 cout << "PIACmd::Do() Usage: fillgd1 nomobj expx expy experry expcut nomgfd" << endl;
513 return(0);
514 }
515 string nomgfd = "";
516 if (tokens.size() > 5) nomgfd = tokens[5];
517 string expy = "";
518 mObjMgr->FillGFD(tokens[0],tokens[1], expy, tokens[2], tokens[3], tokens[4],nomgfd);
519 }
520
521else if (kw == "fillgd2" ) {
522 if (tokens.size() < 6) {
523 cout << "PIACmd::Do() Usage: fillgd2 nomobj expx expy expz experrz expcut nomgfd" << endl;
524 return(0);
525 }
526 string nomgfd = "";
527 if (tokens.size() > 6) nomgfd = tokens[6];
528 mObjMgr->FillGFD(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5],nomgfd);
529 }
530
531// Calcul TF
532else if (kw == "fft" ) {
533 if (tokens.size() < 1) {
534 cout << "PIACmd::Do() Usage: fft nomobj opt" << endl;
535 return(0);
536 }
537 string opt = "n";
538 if (tokens.size() > 1) opt = tokens[1];
539 mObjMgr->FFT(tokens[0], opt);
540 }
541
542// Fit 1D sur objets 1D. Egalement Fit 2D sur objets 2D.
543else if (kw == "fit") {
544 if (tokens.size() < 2) {
545 cout <<"PIACmd::Do() Usage:fit nomobj func \n"
546 <<" [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]\n";
547 return(0);
548 }
549 string p=""; string s=""; string m=""; string M=""; string O="";
550 if (tokens.size()>2)
551 for(int ip=2;ip<tokens.size();ip++) {
552 if(tokens[ip].length()<=2) continue;
553 const char *c = tokens[ip].c_str();
554 if(c[1]!=':') continue;
555 if(c[0]=='p') p=c+2;
556 else if(c[0]=='s') s=c+2;
557 else if(c[0]=='m') m=c+2;
558 else if(c[0]=='M') M=c+2;
559 else if(c[0]=='o') {O += ","; O += c+2;}
560 }
561 mObjMgr->Fit12D(tokens[0],tokens[1],p,s,m,M,O);
562}
563
564#ifndef SANS_StarReco
565// FITPSF sur image
566else if (kw == "fitpsf" ) {
567 if (tokens.size() < 1) {
568 cout << "PIACmd::Do() Usage: fitpsf psf_typ [options_fit]" << endl;
569 return(0);
570 }
571 FitPSFImg(tokens);
572 }
573
574// Photometrie d'ouvertue sur image
575else if (kw == "aper" ) {
576 if (tokens.size() < 1) {
577 cout << "PIACmd::Do() Usage: aper fond_typ [options_calcul]" << endl;
578 return(0);
579 }
580 Photom_Ouv(tokens);
581 }
582#endif
583
584else {
585 cerr << "PIACmd::Do() Erreur - Commande " << kw << " inconuue ! " << endl;
586 return(-1);
587 }
588if (timing) gltimer->Split();
589
590return(0);
591}
592
593
594
595/* --Methode-- */
596int PIACmd::LinkUserFuncs(string& fnameso, string& func1, string& func2, string& func3,
597 string& func4, string& func5)
598{
599string cmd;
600int rc;
601
[171]602if (dynlink) delete dynlink; dynlink = NULL;
[165]603usfmap.clear();
604
[171]605dynlink = new PDynLinkMgr(fnameso, true);
606if (dynlink == NULL) {
607 string sn = fnameso;
[165]608 cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur ouverture SO " << sn << endl;
609 sn = dlerror();
610 cerr << "Loader Error (dlerror()) :" << sn;
611 return(2);
612 }
613
614int nok=0;
615// on utilise donc des DlFunction (Reza 20/08/98) voir commentaire ds .h (pb g++)
616// DlUserProcFunction f = NULL;
617DlFunction f = NULL;
618if ((func1.length() < 1) || (func1 == "-") || (func1 == ".") ) goto fin;
619// f = (DlUserProcFunction) dlsym(dlhandle, func1.c_str());
[171]620f = dynlink->GetFunction(func1);
[165]621if (f) { nok++; usfmap[func1] = f; }
622else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func1 << endl;
623
624if ((func2.length() < 1) || (func2 == "-") || (func2 == ".") ) goto fin;
625// f = (DlUserProcFunction) dlsym(dlhandle, func2.c_str());
[171]626f = dynlink->GetFunction(func2);
[165]627if (f) { nok++; usfmap[func2] = f; }
628else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func2 << endl;
629
630if ((func3.length() < 1) || (func3 == "-") || (func3 == ".") ) goto fin;
631// f = (DlUserProcFunction) dlsym(dlhandle, func3.c_str());
[171]632f = dynlink->GetFunction(func3);
[165]633if (f) { nok++; usfmap[func3] = f; }
634else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func3 << endl;
635
636if ((func4.length() < 1) || (func4 == "-") || (func4 == ".") ) goto fin;
637// f = (DlUserProcFunction) dlsym(dlhandle, func4.c_str());
[171]638f = dynlink->GetFunction(func4);
[165]639if (f) { nok++; usfmap[func4] = f; }
640else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func4 << endl;
641
642if ((func5.length() < 1) || (func5 == "-") || (func5 == ".") ) goto fin;
643// f = (DlUserProcFunction) dlsym(dlhandle, func5.c_str());
[171]644f = dynlink->GetFunction(func5);
[165]645if (f) { nok++; usfmap[func5] = f; }
646else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func5 << endl;
647
648fin:
[171]649if (nok < 1) { if (dynlink) delete dynlink; dynlink = NULL; return(3); }
[165]650else return(0);
651}
652
653/////////////////////////////////////////////////////////////////////////
654/////////////////////////////////////////////////////////////////////////
655/////////////////////////////// Fonctions ///////////////////////////////
656/////////////////////////////////////////////////////////////////////////
657/////////////////////////////////////////////////////////////////////////
658#ifndef SANS_StarReco
659# include "psf.h"
660#include "dynccd.h"
661#include "pimgadapter.h"
662#include "nbmath.h"
663#include "nbtri.h"
664
665/* Nouvelle-Fonction rz + cmv 15/10/98 */
666int FitPSFImg(cmdtok& tokens)
667// Fit de PSF sur une Image.
668//| tokens[0] = Type de PSF ()
669//| : G = fit gaussienne+fond (volume)
670//| : Gi = fit gaussienne+fond integree (volume)
671//| : d = fit DL de gaussienne+fond (volume)
672//| : di = fit DL de gaussienne+fond integree (volume)
673//| : D = fit DL de gaussienne+fond avec coeff variable p6 (volume)
674//| : Di = fit DL de gaussienne+fond integree avec coeff variable p6 (volume)
675//| : M = fit Moffat+fond (expos=p6) (volume)
676//| : Mi = fit Moffat+fond integree (expos=p6) (volume)
677//| tokens[1] = options pour le fit ex: l1.0,G2.,R7.,m10.,M65000.,P3.,f,n,r,i,w10/450
678//| : r = ecriture de l'image des residus (def=non)
679//| : f = ecriture de l'image de de la PSF fitee (def=non)
680//| : i = ecriture du pave initial selectionne pour le fit (def=non)
681//| : n = ecriture de l'image du bruit (def=non)
682//| : Gaa.bbb = gain en e-/ADU (def=1)
683//| : Raa.bbb = readout noise en e- (def=0)
684//| : maa.bbb = limite basse dynamique en ADU (def=-9.e19)
685//| : Maa.bbb = limite haute dynamique en ADU (def=9.e19)
686//| : Paa.bbb = valeur du parametre supplementaire pour
687//| : la.a = print a=FitPSFImg, b=PSF.Fit PSF.IniFit
688//| : wa/b = "a" est le facteur de zoom de la fenetre de representation
689//| des resultats (def=6) et b sa taille maximale (def=500).
690//| Return: 0 si Ok, <0 sinon.
691{
692// Decodage des arguments
693if(tokens.size()<1) return(-1);
694int typ_psf=PSF::psfGaussienne;
695bool wres=false, wpsf=false, wpav=false, wnois=false;
696float gain=1., readout=0., pixmin=-9.e19, pixmax=9.e19;
697int zoom=6, wmax=500;
698float par6=-1.;
699int lp=1, lpg=0;
700string dummy,dum;
701
702if(tokens.size()>=2) {
703 dummy = "," + tokens[1]+ ",";
704 if(strstr(dummy.c_str(),",r,")) wres=true;
705 if(strstr(dummy.c_str(),",f,")) wpsf=true;
706 if(strstr(dummy.c_str(),",i")) wpav=true;
707 if(strstr(dummy.c_str(),",n")) wnois=true;
708 if(strstr(dummy.c_str(),",G")) {
709 size_t p = dummy.find(",G"); size_t q = dummy.find_first_of(',',p+1);
710 dum = dummy.substr(p,q-p);
711 if(dum.length()>2) sscanf(dum.c_str(),",G%f",&gain);
712 }
713 if(strstr(dummy.c_str(),",P")) {
714 size_t p = dummy.find(",P"); size_t q = dummy.find_first_of(',',p+1);
715 dum = dummy.substr(p,q-p);
716 if(dum.length()>2) sscanf(dum.c_str(),",P%f",&par6);
717 }
718 if(strstr(dummy.c_str(),",R")) {
719 size_t p = dummy.find(",R"); size_t q = dummy.find_first_of(',',p+1);
720 dum = dummy.substr(p,q-p);
721 if(dum.length()>2) sscanf(dum.c_str(),",R%f",&readout);
722 }
723 if(strstr(dummy.c_str(),",m")) {
724 size_t p = dummy.find(",m"); size_t q = dummy.find_first_of(',',p+1);
725 dum = dummy.substr(p,q-p);
726 if(dum.length()>2) sscanf(dum.c_str(),",m%f",&pixmin);
727 }
728 if(strstr(dummy.c_str(),",M")) {
729 size_t p = dummy.find(",M"); size_t q = dummy.find_first_of(',',p+1);
730 dum = dummy.substr(p,q-p);
731 if(dum.length()>2) sscanf(dum.c_str(),",M%f",&pixmax);
732 }
733 if(strstr(dummy.c_str(),",l")) { // niveau de print
734 size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
735 dum = dummy.substr(p,q-p);
736 float ab;
737 if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
738 if(ab<0) ab = 0.;
739 lp = (int) ab; lpg = (int) 10.*(ab-lp);
740 }
741 if(strstr(dummy.c_str(),",w")) { // niveau de print
742 size_t p = dummy.find(",w"); size_t q = dummy.find_first_of(',',p+1);
743 dum = dummy.substr(p,q-p);
744 if(dum.length()>2) sscanf(dum.c_str(),",w%d/%d",&zoom,&wmax);
745 }
746}
747
748dummy = "," + tokens[0]+ ",";
749if(strstr(dummy.c_str(),",G,")) typ_psf=PSF::psfGaussienne;
750else if(strstr(dummy.c_str(),",Gi,")) typ_psf=PSF::psfGaussInteg;
751else if(strstr(dummy.c_str(),",d,") ) typ_psf=PSF::psfDelGaussienne;
752else if(strstr(dummy.c_str(),",di,")) typ_psf=PSF::psfDelGaussinteg;
753else if(strstr(dummy.c_str(),",D,") ) {typ_psf=PSF::psfDel1Gaussienne; if(par6<=0.) par6=1.;}
754else if(strstr(dummy.c_str(),",Di,")) {typ_psf=PSF::psfDel1GaussInteg; if(par6<=0.) par6=1.;}
755else if(strstr(dummy.c_str(),",M,") ) {typ_psf=PSF::psfMoffat; if(par6<=1.) par6=2.5;}
756else if(strstr(dummy.c_str(),",Mi,")) {typ_psf=PSF::psfMoffInteg; if(par6<=1.) par6=2.5;}
757
758if(gain<=0.) gain=1.; if(readout<0.) readout=0.;
759if(zoom<1) zoom=1; if(zoom>10) zoom=10; // Il semble que le zoom soit limite a 10 !!!
760if(wmax>500) wmax=500;
761if(lp<0) lp=1; if(lpg<0) lpg=0;
762
763// Wigdet,piimage et image courants
764NamedObjMgr omg;
765PIStdImgApp* app = omg.GetImgApp();
766if(app == NULL)
767 {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
768PIBaseWdg* curw= app->CurrentBaseWdg() ;
769if(curw == NULL)
770 {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
771if(curw->kind() != PIImage::ClassId)
772 {cout<<"Current window not an Image"<<endl; return(-4);}
773PIImage* curpimg = (PIImage*)curw;
774if(curpimg == NULL)
775 {cout<<"Pas de PIImage*"<<endl; return(-5);}
776P2DArrayAdapter* img = curpimg->Image();
777if(img == NULL)
778 {cout<<"Pas de Image*"<<endl; return(-6);}
779string nom = "fitpsf_";
780
781// Definition du pave courant sur la piimage
782int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();
783int x0 = curpimg->XPave()-dx/2, y0 = curpimg->YPave()-dy/2;
784if(x0<0) x0 = 0; if(y0<0) y0 = 0;
785int x1 = x0+dx, y1 = y0+dy;
786if(x1>img->XSize()) x1=img->XSize(); if(y1>img->YSize()) y1=img->YSize();
787dx = x1-x0; dy = y1-y0;
788int dxy = (dx>dy)?dx:dy;
789ImageR4* pav = new ImageR4(dx, dy);
790{for(int j=0;j<dy;j++) for(int i=0;i<dx;i++) (*pav)(i,j)=(*img)(i+x0,j+y0);}
791
792// Bruit sur le pave
793DynCCD mdynccd(kPhotonNoise,pixmin,pixmax,gain,readout);
794ImageR4 *nois = NoiseImage(pav, &mdynccd);
795
796// Print
797int w=2*zoom*dxy+25; if(w>wmax){w=wmax;zoom=(wmax-25)/(2*dxy);if(zoom<1)zoom=1;}
798if(lp>1)
799 {cout<<"FitPSFImg: fun="<<tokens[0]<<","<<typ_psf
800 <<" Write res="<<wres<<" psf="<<wpsf<<" pav="<<wpav<<" nois="<<wnois
801 <<" Par6="<<par6<<" lp="<<lp<<","<<lpg<<endl;
802 mdynccd.Print();
803 cout<<"... dx,y="<<dxy<<" zoom="<<zoom<<" w="<<w<<endl;}
804
805// PSF et Fit
806PSF mypsf(typ_psf,100.0);
807int npar = mypsf.NParm();
808PSFStar mypsfstar(&mypsf);
809if(npar>7) mypsfstar.Parm(6)=(double)par6; // Cas moffat et Dl gauss avec coeff var.
810int rc = mypsfstar.IniFit(*pav,*nois,lpg);
811if(rc != 0)
812 {cout<<"PIDemoApp::FitPSFImg(), Echec IniFit (rc="<<rc<<")"<<endl;
813 delete nois; delete pav; return(-10); }
814rc = mypsfstar.Fit(*pav,*nois,lpg);
815if(rc != 0)
816 {cout<<"PIDemoApp::FitPSFImg(), Echec Fit (rc="<<rc<<")"<<endl;
817 delete nois; delete pav; return(-11); }
818if(lp) mypsfstar.Print(2);
819float smax,axisrat,tiltdeg;
820rc = paramga((float) mypsfstar.SigX(),(float) mypsfstar.SigY()
821 ,(float) mypsfstar.Rho(),&smax,&axisrat,&tiltdeg);
822if(rc==0 && lp)
823 printf("... sa=%g sb=%g (%g) teta=%g\n",smax,smax*axisrat,axisrat,tiltdeg);
824else printf("Erreur paramga \n");
825
826// Creation des pave resultats
827ImageR4* res = new ImageR4(dx, dy);
828ImageR4* psf = new ImageR4(dx, dy);
829CopieImage(*res,*pav); *res -= mypsfstar;
830*psf = 0.; *psf += mypsfstar;
831
832// Si on demande d'ecrire au moins un des resultats (pav,nois,res,psf) => menage!
833bool delobj = (wpav||wnois||wres||wpsf)?true:false;
834
835// Display des resultats, fenetre 2x2 {{pav,nois},{res,psf}}
836app->CreateGraphWin(2,2,w,w); app->SetZoomAtt(zoom);
837dummy = nom+"pav"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
838if(wpav) {
839 // Si on veut garder l'objet
840 omg.AddObj(pav,dummy);
841 omg.DisplayObj(omg.LastObjName());
842} else {
843 // "true" pour tuer (delete) de l'objet ImageR4 lors de la destruction de PIImage
844 app->DispImage(new ImageAdapter<r_4>(pav,true),dummy,Disp_Next);
845}
846dummy = nom+"nois"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
847if(wnois) {
848 omg.AddObj(nois,dummy);
849 omg.DisplayObj(omg.LastObjName());
850} else {
851 app->DispImage(new ImageAdapter<r_4>(nois,true),dummy,Disp_Next);
852}
853dummy = nom+"res"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
854if(wres) {
855 omg.AddObj(res,dummy);
856 omg.DisplayObj(omg.LastObjName());
857} else {
858 app->DispImage(new ImageAdapter<r_4>(res,true),dummy,Disp_Next);
859}
860dummy = nom+"psf"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
861if(wpsf) {
862 omg.AddObj(psf,dummy);
863 omg.DisplayObj(omg.LastObjName());
864} else {
865 app->DispImage(new ImageAdapter<r_4>(psf,true),dummy,Disp_Next);
866}
867
868return(0);
869}
870
871/* Nouvelle-Fonction cmv 16/10/98 */
872int Photom_Ouv(cmdtok& tokens)
873// Photometrie d'ouverture sur une Image.
874//| tokens[0] = Type de calcul du fond et du flux.
875//| : b = fond calcule sur les pixels du bord du pave.
876//| : a = fond sur les pixels autres que ceux servant au calcul du flux.
877//| : faa.bbb = fond fixe a aa.bbb
878//| defaut : "b"
879//| tokens[1] = options pour le fit ex: l1.0
880//| : ra = "a" rayon minimal de scan (def=0)
881//| : ga = "a" taille de la fenetre glissante pour les extrapolations (def=4)
882//| : la = "a" niveau de print
883//| tokens[2] = options graphiques
884//| Return: 0 si Ok, <0 sinon.
885{
886// Decodage des arguments
887if(tokens.size()<1) return(-1);
888int typ_fond=1, rmin=0, rmax, npt_gliss=4, lp=0;
889float fondimpos=0., fondbord=0., fondautre;
890string opt="";
891string dummy; string dum;
892
893if(tokens.size()>=3) opt = tokens[2];
894if(tokens.size()>=2) {
895 dummy = "," + tokens[1]+ ",";
896 if(strstr(dummy.c_str(),",r")) {
897 size_t p = dummy.find(",r"); size_t q = dummy.find_first_of(',',p+1);
898 dum = dummy.substr(p,q-p);
899 if(dum.length()>2) sscanf(dum.c_str(),",r%d",&lp);
900 }
901 if(strstr(dummy.c_str(),",g")) {
902 size_t p = dummy.find(",g"); size_t q = dummy.find_first_of(',',p+1);
903 dum = dummy.substr(p,q-p);
904 if(dum.length()>2) sscanf(dum.c_str(),",g%d",&npt_gliss);
905 }
906 if(strstr(dummy.c_str(),",l")) {
907 size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
908 dum = dummy.substr(p,q-p);
909 if(dum.length()>2) sscanf(dum.c_str(),",l%d",&lp);
910 }
911}
912
913dummy = "," + tokens[0]+ ",";
914if(strstr(dummy.c_str(),",b,")) typ_fond=1;
915else if(strstr(dummy.c_str(),",a,")) typ_fond=2;
916else if(strstr(dummy.c_str(),",f")) {
917 typ_fond=3;
918 size_t p = dummy.find(",f"); size_t q = dummy.find_first_of(',',p+1);
919 dum = dummy.substr(p,q-p);
920 if(dum.length()>2) sscanf(dum.c_str(),",f%f",&fondimpos);
921}
922
923if(rmin<0) rmin=0; if(npt_gliss<=0) npt_gliss=4; if(lp<0) lp=1;
924
925// Wigdet,piimage et image courants
926NamedObjMgr omg;
927PIStdImgApp* app = omg.GetImgApp();
928if(app == NULL)
929 {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
930PIBaseWdg* curw= app->CurrentBaseWdg() ;
931if(curw == NULL)
932 {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
933if(curw->kind() != PIImage::ClassId)
934 {cout<<"Current window not an Image"<<endl; return(-4);}
935PIImage* curpimg = (PIImage*)curw;
936if(curpimg == NULL)
937 {cout<<"Pas de PIImage*"<<endl; return(-5);}
938P2DArrayAdapter* img = curpimg->Image();
939if(img == NULL)
940 {cout<<"Pas de Image*"<<endl; return(-6);}
941string nom = "ouver_";
942
943// Definition du pave courant sur la piimage
944int xc = curpimg->XPave(), yc = curpimg->YPave();
945int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();
946int dxym = (dx<dy)?dx:dy; if(dxym%2==0) dxym++; rmax = dxym/2;
947if(dxym<3||dxym>img->XSize()||dxym>img->YSize())
948 {cout<<"Pave trop petit/grand: "<<dxym<<endl; return(-10);}
949if(rmin>=rmax)
950 {cout<<"rayon mini="<<rmin<<">= rmax="<<rmax<<endl; return(-11);}
951int x0=xc-rmax, y0=yc-rmax, x1=xc+rmax, y1=yc+rmax;
952if(x0<0||y0<0||x1>=img->XSize()||y1>=img->YSize())
953 {cout<<"Pave hors image"<<endl; return(-12);}
954int ic = rmax, jc = rmax;
955ImageR4* pav = new ImageR4(dxym, dxym);
956{for(int j=-rmax;j<=rmax;j++) for(int i=-rmax;i<=rmax;i++)
957 (*pav)(ic+i,jc+j)=(*img)(xc+i,yc+j);}
958
959// Print
960if(lp>1)
961 cout<<"photom_ouv[typ="<<typ_fond<<"]: pav["<<dx<<","<<dy<<" ("<<dxym
962 <<"] xc,yc="<<xc<<","<<yc<<"] ic,jc="<<ic<<","<<jc<<endl
963 <<"] (fond="<<fondimpos<<") pxmil="<<(*pav)(ic,ic)
964 <<" -> Rph de "<<rmin<<" a "<<rmax<<endl;
965
966// Ouverture des buffers et objets necessaires
967float *p = new float[dxym*dxym];
968r_4 xnt[4];
969NTuple *ntfl=NULL, *ntex=NULL;
970{char* name[4]={"npx","flux","fond","mag"}; ntfl = new NTuple(4,name);}
971{char* name[3]={"k","a","f0"}; ntex = new NTuple(3,name);}
972
973// Calcul du fond sur les pixels du bord du pave
974{
975int nf = 0;
976for(int i=0;i<dxym;i++) {p[nf]=(*pav)(i,0); nf++; p[nf]=(*pav)(i,dxym-1); nf++;}
977for(int j=1;j<dxym-1;j++) {p[nf]=(*pav)(0,j); nf++; p[nf]=(*pav)(dxym-1,j); nf++;}
978qsort(p,(size_t) nf,sizeof(float),qSort_Float);
979fondbord = (p[nf/2]+p[(nf-1)/2])/2.;
980if(lp>0) cout<<"Fond_Bord: "<<fondbord<<" (avec 1 px au bord="<<nf<<")"<<endl;
981}
982
983// Photometrie d'ouverture variable.
984int rfin = -1;
985for(int r=rmin;r<=rmax;r++) {
986 double sum=0.; int ns=0, nf=0;
987 for(int i=0;i<dxym;i++) for(int j=0;j<dxym;j++) {
988 if(ic-r<=i && i<=ic+r && jc-r<=j && j<=jc+r ) {sum += (*pav)(i,j); ns++;}
989 else {p[nf] = (*pav)(i,j); nf++;}
990 }
991 if(typ_fond==2) {
992 if(nf<2) continue;
993 qsort(p,(size_t) nf,sizeof(float),qSort_Float);
994 fondautre = (p[nf/2]+p[(nf-1)/2])/2.;
995 }
996 double fond=-1.e19;
997 if(typ_fond==3) fond=fondimpos; else if(typ_fond==2) fond=fondautre; else fond=fondbord;
998 rfin = r;
999 xnt[0] = ns;
1000 xnt[1] = sum-fond*ns;
1001 xnt[2] = fond;
1002 xnt[3] = (xnt[1]>0.)? -2.5*log10((double)xnt[1]): -99.;
1003 if(lp>0)
1004 cout<<"Ouv[r="<<r<<",ns="<<xnt[0]<<",nf="<<nf<<"]: fd="<<xnt[2]
1005 <<" fl="<<xnt[1]<<" mag="<<xnt[3]<<endl;
1006 ntfl->Fill(xnt);
1007}
1008delete [] p;
1009int nr = ntfl->NEntry();
1010if(lp>1) cout<<"Nombre de calculs de flux: "<<nr<<endl;
1011if(nr<=0)
1012 {cout<<"Pas de flux calcules"<<endl; delete pav; delete ntfl; delete ntex; return(-13);}
1013
1014cout<<"Ouv[rouv="<<rfin<<",npix="<<xnt[0]<<"]: fond="<<xnt[2]
1015 <<" flux="<<xnt[1]<<" mag="<<xnt[3]<<endl;
1016
1017// Extrapolation: F_mes = F_vrai + npxflux*(Fd_vrai-Fd_mes)
1018// On fit F_mes=f(npxflux) sur une fenetre glissante de npt_gliss pts.
1019bool ex = false;
1020if(nr>=npt_gliss) {
1021 ex = true;
1022 double *npix = new double[nr], *flux = new double[nr], *err = new double[nr];
1023 for(int i=0;i<nr;i++) { ntfl->GetVec(i,xnt); npix[i]=xnt[0]; flux[i]=xnt[1]; err[i]=1.;}
1024 double a,b; int n;
1025 for(int k=0;k<=nr-npt_gliss;k++) {
1026 n = npt_gliss; FitLin(&npix[k],&flux[k],&err[k],&n,&b,&a);
1027 xnt[0]=k; xnt[1]=a; xnt[2]=b; ntex->Fill(xnt);
1028 if(lp>1)
1029 cout<<"Fit["<<k<<","<<k+npt_gliss-1<<"] f0="<<xnt[2]<<" pente="<<xnt[1]<<endl;
1030 }
1031 delete [] npix; delete [] flux; delete [] err;
1032}
1033
1034// Display des resultats, fenetre 2x2 {{flux%ns,fond%ns},{f0,pente}}
1035if(ex) app->CreateGraphWin(2,2); else app->CreateGraphWin(1,2);
1036dummy = nom+"*"; omg.DelObjects(dummy);
1037
1038dummy=nom+"pav"; omg.AddObj(pav,dummy);
1039
1040string xd=""; string yd=""; string zd,exd,eyd,ezd;
1041dummy=nom+"flux"; omg.AddObj(ntfl,dummy);
1042xd="npx"; yd="flux";
1043omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1044if(typ_fond==2)
1045 {xd="npx"; yd="fond"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
1046if(ex) {
1047 dummy=nom+"fit"; omg.AddObj(ntex,dummy);
1048 if(typ_fond!=2)
1049 {xd="a"; yd="f0"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
1050 xd = "k"; yd = "f0";
1051 omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1052 xd = "k"; yd = "a";
1053 omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1054}
1055
1056return(0);
1057}
1058#endif
Note: See TracBrowser for help on using the repository browser.