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

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

Petits pb de compil aCC Reza 23/12/98

File size: 37.0 KB
Line 
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;
46dynlink = NULL;
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] -- 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]" << endl; return(0); }
237 string sph = "";
238 for(int gg=0; gg<5; gg++) tokens.push_back(sph);
239 int rc = LinkUserFuncs(tokens[0], tokens[1], tokens[2], tokens[3]);
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
602if (dynlink) delete dynlink; dynlink = NULL;
603usfmap.clear();
604
605dynlink = new PDynLinkMgr(fnameso, true);
606if (dynlink == NULL) {
607 string sn = fnameso;
608 cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur ouverture SO " << sn << endl;
609 return(2);
610 }
611
612int nok=0;
613// on utilise donc des DlFunction (Reza 20/08/98) voir commentaire ds .h (pb g++)
614// DlUserProcFunction f = NULL;
615DlFunction f = NULL;
616if ((func1.length() < 1) || (func1 == "-") || (func1 == ".") ) goto fin;
617// f = (DlUserProcFunction) dlsym(dlhandle, func1.c_str());
618f = dynlink->GetFunction(func1);
619if (f) { nok++; usfmap[func1] = f; }
620else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func1 << endl;
621
622if ((func2.length() < 1) || (func2 == "-") || (func2 == ".") ) goto fin;
623// f = (DlUserProcFunction) dlsym(dlhandle, func2.c_str());
624f = dynlink->GetFunction(func2);
625if (f) { nok++; usfmap[func2] = f; }
626else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func2 << endl;
627
628if ((func3.length() < 1) || (func3 == "-") || (func3 == ".") ) goto fin;
629// f = (DlUserProcFunction) dlsym(dlhandle, func3.c_str());
630f = dynlink->GetFunction(func3);
631if (f) { nok++; usfmap[func3] = f; }
632else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func3 << endl;
633
634/* Pb compile g++ 2.7.2
635if ((func4.length() < 1) || (func4 == "-") || (func4 == ".") ) goto fin;
636// f = (DlUserProcFunction) dlsym(dlhandle, func4.c_str());
637f = dynlink->GetFunction(func4);
638if (f) { nok++; usfmap[func4] = f; }
639else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func4 << endl;
640
641if ((func5.length() < 1) || (func5 == "-") || (func5 == ".") ) goto fin;
642// f = (DlUserProcFunction) dlsym(dlhandle, func5.c_str());
643f = dynlink->GetFunction(func5);
644if (f) { nok++; usfmap[func5] = f; }
645else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func5 << endl;
646*/
647fin:
648if (nok < 1) { if (dynlink) delete dynlink; dynlink = NULL; return(3); }
649else return(0);
650}
651
652/////////////////////////////////////////////////////////////////////////
653/////////////////////////////////////////////////////////////////////////
654/////////////////////////////// Fonctions ///////////////////////////////
655/////////////////////////////////////////////////////////////////////////
656/////////////////////////////////////////////////////////////////////////
657#ifndef SANS_StarReco
658# include "psf.h"
659#include "dynccd.h"
660#include "pimgadapter.h"
661#include "nbmath.h"
662#include "nbtri.h"
663
664/* Nouvelle-Fonction rz + cmv 15/10/98 */
665int FitPSFImg(cmdtok& tokens)
666// Fit de PSF sur une Image.
667//| tokens[0] = Type de PSF ()
668//| : G = fit gaussienne+fond (volume)
669//| : Gi = fit gaussienne+fond integree (volume)
670//| : d = fit DL de gaussienne+fond (volume)
671//| : di = fit DL de gaussienne+fond integree (volume)
672//| : D = fit DL de gaussienne+fond avec coeff variable p6 (volume)
673//| : Di = fit DL de gaussienne+fond integree avec coeff variable p6 (volume)
674//| : M = fit Moffat+fond (expos=p6) (volume)
675//| : Mi = fit Moffat+fond integree (expos=p6) (volume)
676//| tokens[1] = options pour le fit ex: l1.0,G2.,R7.,m10.,M65000.,P3.,f,n,r,i,w10/450
677//| : r = ecriture de l'image des residus (def=non)
678//| : f = ecriture de l'image de de la PSF fitee (def=non)
679//| : i = ecriture du pave initial selectionne pour le fit (def=non)
680//| : n = ecriture de l'image du bruit (def=non)
681//| : Gaa.bbb = gain en e-/ADU (def=1)
682//| : Raa.bbb = readout noise en e- (def=0)
683//| : maa.bbb = limite basse dynamique en ADU (def=-9.e19)
684//| : Maa.bbb = limite haute dynamique en ADU (def=9.e19)
685//| : Paa.bbb = valeur du parametre supplementaire pour
686//| : la.a = print a=FitPSFImg, b=PSF.Fit PSF.IniFit
687//| : wa/b = "a" est le facteur de zoom de la fenetre de representation
688//| des resultats (def=6) et b sa taille maximale (def=500).
689//| Return: 0 si Ok, <0 sinon.
690{
691// Decodage des arguments
692if(tokens.size()<1) return(-1);
693int typ_psf=PSF::psfGaussienne;
694bool wres=false, wpsf=false, wpav=false, wnois=false;
695float gain=1., readout=0., pixmin=-9.e19, pixmax=9.e19;
696int zoom=6, wmax=500;
697float par6=-1.;
698int lp=1, lpg=0;
699string dummy,dum;
700
701if(tokens.size()>=2) {
702 dummy = "," + tokens[1]+ ",";
703 if(strstr(dummy.c_str(),",r,")) wres=true;
704 if(strstr(dummy.c_str(),",f,")) wpsf=true;
705 if(strstr(dummy.c_str(),",i")) wpav=true;
706 if(strstr(dummy.c_str(),",n")) wnois=true;
707 if(strstr(dummy.c_str(),",G")) {
708 size_t p = dummy.find(",G"); size_t q = dummy.find_first_of(',',p+1);
709 dum = dummy.substr(p,q-p);
710 if(dum.length()>2) sscanf(dum.c_str(),",G%f",&gain);
711 }
712 if(strstr(dummy.c_str(),",P")) {
713 size_t p = dummy.find(",P"); size_t q = dummy.find_first_of(',',p+1);
714 dum = dummy.substr(p,q-p);
715 if(dum.length()>2) sscanf(dum.c_str(),",P%f",&par6);
716 }
717 if(strstr(dummy.c_str(),",R")) {
718 size_t p = dummy.find(",R"); size_t q = dummy.find_first_of(',',p+1);
719 dum = dummy.substr(p,q-p);
720 if(dum.length()>2) sscanf(dum.c_str(),",R%f",&readout);
721 }
722 if(strstr(dummy.c_str(),",m")) {
723 size_t p = dummy.find(",m"); size_t q = dummy.find_first_of(',',p+1);
724 dum = dummy.substr(p,q-p);
725 if(dum.length()>2) sscanf(dum.c_str(),",m%f",&pixmin);
726 }
727 if(strstr(dummy.c_str(),",M")) {
728 size_t p = dummy.find(",M"); size_t q = dummy.find_first_of(',',p+1);
729 dum = dummy.substr(p,q-p);
730 if(dum.length()>2) sscanf(dum.c_str(),",M%f",&pixmax);
731 }
732 if(strstr(dummy.c_str(),",l")) { // niveau de print
733 size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
734 dum = dummy.substr(p,q-p);
735 float ab;
736 if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
737 if(ab<0) ab = 0.;
738 lp = (int) ab; lpg = (int) 10.*(ab-lp);
739 }
740 if(strstr(dummy.c_str(),",w")) { // niveau de print
741 size_t p = dummy.find(",w"); size_t q = dummy.find_first_of(',',p+1);
742 dum = dummy.substr(p,q-p);
743 if(dum.length()>2) sscanf(dum.c_str(),",w%d/%d",&zoom,&wmax);
744 }
745}
746
747dummy = "," + tokens[0]+ ",";
748if(strstr(dummy.c_str(),",G,")) typ_psf=PSF::psfGaussienne;
749else if(strstr(dummy.c_str(),",Gi,")) typ_psf=PSF::psfGaussInteg;
750else if(strstr(dummy.c_str(),",d,") ) typ_psf=PSF::psfDelGaussienne;
751else if(strstr(dummy.c_str(),",di,")) typ_psf=PSF::psfDelGaussinteg;
752else if(strstr(dummy.c_str(),",D,") ) {typ_psf=PSF::psfDel1Gaussienne; if(par6<=0.) par6=1.;}
753else if(strstr(dummy.c_str(),",Di,")) {typ_psf=PSF::psfDel1GaussInteg; if(par6<=0.) par6=1.;}
754else if(strstr(dummy.c_str(),",M,") ) {typ_psf=PSF::psfMoffat; if(par6<=1.) par6=2.5;}
755else if(strstr(dummy.c_str(),",Mi,")) {typ_psf=PSF::psfMoffInteg; if(par6<=1.) par6=2.5;}
756
757if(gain<=0.) gain=1.; if(readout<0.) readout=0.;
758if(zoom<1) zoom=1; if(zoom>10) zoom=10; // Il semble que le zoom soit limite a 10 !!!
759if(wmax>500) wmax=500;
760if(lp<0) lp=1; if(lpg<0) lpg=0;
761
762// Wigdet,piimage et image courants
763NamedObjMgr omg;
764PIStdImgApp* app = omg.GetImgApp();
765if(app == NULL)
766 {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
767PIBaseWdg* curw= app->CurrentBaseWdg() ;
768if(curw == NULL)
769 {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
770if(curw->kind() != PIImage::ClassId)
771 {cout<<"Current window not an Image"<<endl; return(-4);}
772PIImage* curpimg = (PIImage*)curw;
773if(curpimg == NULL)
774 {cout<<"Pas de PIImage*"<<endl; return(-5);}
775P2DArrayAdapter* img = curpimg->Image();
776if(img == NULL)
777 {cout<<"Pas de Image*"<<endl; return(-6);}
778string nom = "fitpsf_";
779
780// Definition du pave courant sur la piimage
781int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();
782int x0 = curpimg->XPave()-dx/2, y0 = curpimg->YPave()-dy/2;
783if(x0<0) x0 = 0; if(y0<0) y0 = 0;
784int x1 = x0+dx, y1 = y0+dy;
785if(x1>img->XSize()) x1=img->XSize(); if(y1>img->YSize()) y1=img->YSize();
786dx = x1-x0; dy = y1-y0;
787int dxy = (dx>dy)?dx:dy;
788ImageR4* pav = new ImageR4(dx, dy);
789{for(int j=0;j<dy;j++) for(int i=0;i<dx;i++) (*pav)(i,j)=(*img)(i+x0,j+y0);}
790
791// Bruit sur le pave
792DynCCD mdynccd(kPhotonNoise,pixmin,pixmax,gain,readout);
793ImageR4 *nois = NoiseImage(pav, &mdynccd);
794
795// Print
796int w=2*zoom*dxy+25; if(w>wmax){w=wmax;zoom=(wmax-25)/(2*dxy);if(zoom<1)zoom=1;}
797if(lp>1)
798 {cout<<"FitPSFImg: fun="<<tokens[0]<<","<<typ_psf
799 <<" Write res="<<wres<<" psf="<<wpsf<<" pav="<<wpav<<" nois="<<wnois
800 <<" Par6="<<par6<<" lp="<<lp<<","<<lpg<<endl;
801 mdynccd.Print();
802 cout<<"... dx,y="<<dxy<<" zoom="<<zoom<<" w="<<w<<endl;}
803
804// PSF et Fit
805PSF mypsf(typ_psf,100.0);
806int npar = mypsf.NParm();
807PSFStar mypsfstar(&mypsf);
808if(npar>7) mypsfstar.Parm(6)=(double)par6; // Cas moffat et Dl gauss avec coeff var.
809int rc = mypsfstar.IniFit(*pav,*nois,lpg);
810if(rc != 0)
811 {cout<<"PIDemoApp::FitPSFImg(), Echec IniFit (rc="<<rc<<")"<<endl;
812 delete nois; delete pav; return(-10); }
813rc = mypsfstar.Fit(*pav,*nois,lpg);
814if(rc != 0)
815 {cout<<"PIDemoApp::FitPSFImg(), Echec Fit (rc="<<rc<<")"<<endl;
816 delete nois; delete pav; return(-11); }
817if(lp) mypsfstar.Print(2);
818float smax,axisrat,tiltdeg;
819rc = paramga((float) mypsfstar.SigX(),(float) mypsfstar.SigY()
820 ,(float) mypsfstar.Rho(),&smax,&axisrat,&tiltdeg);
821if(rc==0 && lp)
822 printf("... sa=%g sb=%g (%g) teta=%g\n",smax,smax*axisrat,axisrat,tiltdeg);
823else printf("Erreur paramga \n");
824
825// Creation des pave resultats
826ImageR4* res = new ImageR4(dx, dy);
827ImageR4* psf = new ImageR4(dx, dy);
828CopieImage(*res,*pav); *res -= mypsfstar;
829*psf = 0.; *psf += mypsfstar;
830
831// Si on demande d'ecrire au moins un des resultats (pav,nois,res,psf) => menage!
832bool delobj = (wpav||wnois||wres||wpsf)?true:false;
833
834// Display des resultats, fenetre 2x2 {{pav,nois},{res,psf}}
835app->CreateGraphWin(2,2,w,w); app->SetZoomAtt(zoom);
836dummy = nom+"pav"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
837if(wpav) {
838 // Si on veut garder l'objet
839 omg.AddObj(pav,dummy);
840 omg.DisplayObj(omg.LastObjName());
841} else {
842 // "true" pour tuer (delete) de l'objet ImageR4 lors de la destruction de PIImage
843 app->DispImage(new ImageAdapter<r_4>(pav,true),dummy,Disp_Next);
844}
845dummy = nom+"nois"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
846if(wnois) {
847 omg.AddObj(nois,dummy);
848 omg.DisplayObj(omg.LastObjName());
849} else {
850 app->DispImage(new ImageAdapter<r_4>(nois,true),dummy,Disp_Next);
851}
852dummy = nom+"res"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
853if(wres) {
854 omg.AddObj(res,dummy);
855 omg.DisplayObj(omg.LastObjName());
856} else {
857 app->DispImage(new ImageAdapter<r_4>(res,true),dummy,Disp_Next);
858}
859dummy = nom+"psf"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
860if(wpsf) {
861 omg.AddObj(psf,dummy);
862 omg.DisplayObj(omg.LastObjName());
863} else {
864 app->DispImage(new ImageAdapter<r_4>(psf,true),dummy,Disp_Next);
865}
866
867return(0);
868}
869
870/* Nouvelle-Fonction cmv 16/10/98 */
871int Photom_Ouv(cmdtok& tokens)
872// Photometrie d'ouverture sur une Image.
873//| tokens[0] = Type de calcul du fond et du flux.
874//| : b = fond calcule sur les pixels du bord du pave.
875//| : a = fond sur les pixels autres que ceux servant au calcul du flux.
876//| : faa.bbb = fond fixe a aa.bbb
877//| defaut : "b"
878//| tokens[1] = options pour le fit ex: l1.0
879//| : ra = "a" rayon minimal de scan (def=0)
880//| : ga = "a" taille de la fenetre glissante pour les extrapolations (def=4)
881//| : la = "a" niveau de print
882//| tokens[2] = options graphiques
883//| Return: 0 si Ok, <0 sinon.
884{
885// Decodage des arguments
886if(tokens.size()<1) return(-1);
887int typ_fond=1, rmin=0, rmax, npt_gliss=4, lp=0;
888float fondimpos=0., fondbord=0., fondautre;
889string opt="";
890string dummy; string dum;
891
892if(tokens.size()>=3) opt = tokens[2];
893if(tokens.size()>=2) {
894 dummy = "," + tokens[1]+ ",";
895 if(strstr(dummy.c_str(),",r")) {
896 size_t p = dummy.find(",r"); size_t q = dummy.find_first_of(',',p+1);
897 dum = dummy.substr(p,q-p);
898 if(dum.length()>2) sscanf(dum.c_str(),",r%d",&lp);
899 }
900 if(strstr(dummy.c_str(),",g")) {
901 size_t p = dummy.find(",g"); size_t q = dummy.find_first_of(',',p+1);
902 dum = dummy.substr(p,q-p);
903 if(dum.length()>2) sscanf(dum.c_str(),",g%d",&npt_gliss);
904 }
905 if(strstr(dummy.c_str(),",l")) {
906 size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
907 dum = dummy.substr(p,q-p);
908 if(dum.length()>2) sscanf(dum.c_str(),",l%d",&lp);
909 }
910}
911
912dummy = "," + tokens[0]+ ",";
913if(strstr(dummy.c_str(),",b,")) typ_fond=1;
914else if(strstr(dummy.c_str(),",a,")) typ_fond=2;
915else if(strstr(dummy.c_str(),",f")) {
916 typ_fond=3;
917 size_t p = dummy.find(",f"); size_t q = dummy.find_first_of(',',p+1);
918 dum = dummy.substr(p,q-p);
919 if(dum.length()>2) sscanf(dum.c_str(),",f%f",&fondimpos);
920}
921
922if(rmin<0) rmin=0; if(npt_gliss<=0) npt_gliss=4; if(lp<0) lp=1;
923
924// Wigdet,piimage et image courants
925NamedObjMgr omg;
926PIStdImgApp* app = omg.GetImgApp();
927if(app == NULL)
928 {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
929PIBaseWdg* curw= app->CurrentBaseWdg() ;
930if(curw == NULL)
931 {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
932if(curw->kind() != PIImage::ClassId)
933 {cout<<"Current window not an Image"<<endl; return(-4);}
934PIImage* curpimg = (PIImage*)curw;
935if(curpimg == NULL)
936 {cout<<"Pas de PIImage*"<<endl; return(-5);}
937P2DArrayAdapter* img = curpimg->Image();
938if(img == NULL)
939 {cout<<"Pas de Image*"<<endl; return(-6);}
940string nom = "ouver_";
941
942// Definition du pave courant sur la piimage
943int xc = curpimg->XPave(), yc = curpimg->YPave();
944int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();
945int dxym = (dx<dy)?dx:dy; if(dxym%2==0) dxym++; rmax = dxym/2;
946if(dxym<3||dxym>img->XSize()||dxym>img->YSize())
947 {cout<<"Pave trop petit/grand: "<<dxym<<endl; return(-10);}
948if(rmin>=rmax)
949 {cout<<"rayon mini="<<rmin<<">= rmax="<<rmax<<endl; return(-11);}
950int x0=xc-rmax, y0=yc-rmax, x1=xc+rmax, y1=yc+rmax;
951if(x0<0||y0<0||x1>=img->XSize()||y1>=img->YSize())
952 {cout<<"Pave hors image"<<endl; return(-12);}
953int ic = rmax, jc = rmax;
954ImageR4* pav = new ImageR4(dxym, dxym);
955{for(int j=-rmax;j<=rmax;j++) for(int i=-rmax;i<=rmax;i++)
956 (*pav)(ic+i,jc+j)=(*img)(xc+i,yc+j);}
957
958// Print
959if(lp>1)
960 cout<<"photom_ouv[typ="<<typ_fond<<"]: pav["<<dx<<","<<dy<<" ("<<dxym
961 <<"] xc,yc="<<xc<<","<<yc<<"] ic,jc="<<ic<<","<<jc<<endl
962 <<"] (fond="<<fondimpos<<") pxmil="<<(*pav)(ic,ic)
963 <<" -> Rph de "<<rmin<<" a "<<rmax<<endl;
964
965// Ouverture des buffers et objets necessaires
966float *p = new float[dxym*dxym];
967r_4 xnt[4];
968NTuple *ntfl=NULL, *ntex=NULL;
969{char* name[4]={"npx","flux","fond","mag"}; ntfl = new NTuple(4,name);}
970{char* name[3]={"k","a","f0"}; ntex = new NTuple(3,name);}
971
972// Calcul du fond sur les pixels du bord du pave
973{
974int nf = 0;
975for(int i=0;i<dxym;i++) {p[nf]=(*pav)(i,0); nf++; p[nf]=(*pav)(i,dxym-1); nf++;}
976for(int j=1;j<dxym-1;j++) {p[nf]=(*pav)(0,j); nf++; p[nf]=(*pav)(dxym-1,j); nf++;}
977qsort(p,(size_t) nf,sizeof(float),qSort_Float);
978fondbord = (p[nf/2]+p[(nf-1)/2])/2.;
979if(lp>0) cout<<"Fond_Bord: "<<fondbord<<" (avec 1 px au bord="<<nf<<")"<<endl;
980}
981
982// Photometrie d'ouverture variable.
983int rfin = -1;
984for(int r=rmin;r<=rmax;r++) {
985 double sum=0.; int ns=0, nf=0;
986 for(int i=0;i<dxym;i++) for(int j=0;j<dxym;j++) {
987 if(ic-r<=i && i<=ic+r && jc-r<=j && j<=jc+r ) {sum += (*pav)(i,j); ns++;}
988 else {p[nf] = (*pav)(i,j); nf++;}
989 }
990 if(typ_fond==2) {
991 if(nf<2) continue;
992 qsort(p,(size_t) nf,sizeof(float),qSort_Float);
993 fondautre = (p[nf/2]+p[(nf-1)/2])/2.;
994 }
995 double fond=-1.e19;
996 if(typ_fond==3) fond=fondimpos; else if(typ_fond==2) fond=fondautre; else fond=fondbord;
997 rfin = r;
998 xnt[0] = ns;
999 xnt[1] = sum-fond*ns;
1000 xnt[2] = fond;
1001 xnt[3] = (xnt[1]>0.)? -2.5*log10((double)xnt[1]): -99.;
1002 if(lp>0)
1003 cout<<"Ouv[r="<<r<<",ns="<<xnt[0]<<",nf="<<nf<<"]: fd="<<xnt[2]
1004 <<" fl="<<xnt[1]<<" mag="<<xnt[3]<<endl;
1005 ntfl->Fill(xnt);
1006}
1007delete [] p;
1008int nr = ntfl->NEntry();
1009if(lp>1) cout<<"Nombre de calculs de flux: "<<nr<<endl;
1010if(nr<=0)
1011 {cout<<"Pas de flux calcules"<<endl; delete pav; delete ntfl; delete ntex; return(-13);}
1012
1013cout<<"Ouv[rouv="<<rfin<<",npix="<<xnt[0]<<"]: fond="<<xnt[2]
1014 <<" flux="<<xnt[1]<<" mag="<<xnt[3]<<endl;
1015
1016// Extrapolation: F_mes = F_vrai + npxflux*(Fd_vrai-Fd_mes)
1017// On fit F_mes=f(npxflux) sur une fenetre glissante de npt_gliss pts.
1018bool ex = false;
1019if(nr>=npt_gliss) {
1020 ex = true;
1021 double *npix = new double[nr], *flux = new double[nr], *err = new double[nr];
1022 for(int i=0;i<nr;i++) { ntfl->GetVec(i,xnt); npix[i]=xnt[0]; flux[i]=xnt[1]; err[i]=1.;}
1023 double a,b; int n;
1024 for(int k=0;k<=nr-npt_gliss;k++) {
1025 n = npt_gliss; FitLin(&npix[k],&flux[k],&err[k],&n,&b,&a);
1026 xnt[0]=k; xnt[1]=a; xnt[2]=b; ntex->Fill(xnt);
1027 if(lp>1)
1028 cout<<"Fit["<<k<<","<<k+npt_gliss-1<<"] f0="<<xnt[2]<<" pente="<<xnt[1]<<endl;
1029 }
1030 delete [] npix; delete [] flux; delete [] err;
1031}
1032
1033// Display des resultats, fenetre 2x2 {{flux%ns,fond%ns},{f0,pente}}
1034if(ex) app->CreateGraphWin(2,2); else app->CreateGraphWin(1,2);
1035dummy = nom+"*"; omg.DelObjects(dummy);
1036
1037dummy=nom+"pav"; omg.AddObj(pav,dummy);
1038
1039string xd=""; string yd=""; string zd,exd,eyd,ezd;
1040dummy=nom+"flux"; omg.AddObj(ntfl,dummy);
1041xd="npx"; yd="flux";
1042omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1043if(typ_fond==2)
1044 {xd="npx"; yd="fond"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
1045if(ex) {
1046 dummy=nom+"fit"; omg.AddObj(ntex,dummy);
1047 if(typ_fond!=2)
1048 {xd="a"; yd="f0"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
1049 xd = "k"; yd = "f0";
1050 omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1051 xd = "k"; yd = "a";
1052 omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
1053}
1054
1055return(0);
1056}
1057#endif
Note: See TracBrowser for help on using the repository browser.