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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <iostream.h>
#include <string>
#include <vector>
#include <map>
#if defined(__KCC__)
using std::string ;
#include <vector.h>
#include <map.h>
#endif

#include "ctimer.h"

#include "nobjmgr.h"
#include "piacmd.h"
#include "pistdimgapp.h"

#include "histos.h"
#include "histos2.h"
#include "hisprof.h"
#include "ntuple.h"
#include "generaldata.h"


//  Pour le decoupage des commandes en lignes 
typedef vector<string> cmdtok;

//  Fonction pour faire le Fit de PSF sur image   
static int FitPSFImg(cmdtok& tokens);
static int Photom_Ouv(cmdtok& tokens);


/* --Methode-- */
PIACmd::PIACmd(NamedObjMgr* omg, PIStdImgApp* app)
{
if (omg)  { mOmg= false; mObjMgr = omg; }
else { mObjMgr = new NamedObjMgr;  mOmg = true; }
mImgApp = app;
system("cp history.pic hisold.pic");
hist.open("history.pic");
trace = false;   timing = false;
gltimer = NULL;
dynlink = NULL;
}

/* --Methode-- */
PIACmd::~PIACmd()
{
if (mOmg)  delete mObjMgr;
hist.close();
if (gltimer) { delete gltimer;  gltimer = NULL; }
}

/* --Methode-- */
int PIACmd::Exec(string& file)
{
char line_buff[512];
FILE *fip;

if ( (fip = fopen(file.c_str(),"r")) == NULL ) {
  cerr << "PIACmd::Exec() Error opening file " << file << endl;
  hist << "##! PIACmd::Exec() Error opening file " << file << endl;
  return(0);
  }
 
hist << "### Executing commands from " << file << endl;
if (trace) { 
  mImgApp->GetConsole()->AddStr("### Executing commands from ", PIVA_Magenta);
  mImgApp->GetConsole()->AddStr(file.c_str(), PIVA_Magenta);
  mImgApp->GetConsole()->AddStr("\n", PIVA_Magenta);
  }

while (fgets(line_buff,511,fip) != NULL)
  {
  if (trace) mImgApp->GetConsole()->AddStr(line_buff, PIVA_Magenta);
  line_buff[strlen(line_buff)-1] = '\0';   /*  LF/CR de la fin */
  string line(line_buff);
  Do(line);
  }
hist << "### End of Exec( " << file << " ) " << endl;
if (trace) { 
  mImgApp->GetConsole()->AddStr("### End of Exec( ", PIVA_Magenta);
  mImgApp->GetConsole()->AddStr(file.c_str(), PIVA_Magenta);
  mImgApp->GetConsole()->AddStr(" ) \n", PIVA_Magenta);
  }

return(0);
}

/* --Methode-- */
int PIACmd::Do(string& s)
{

cmdtok tokens;
if (s.length() < 1)  return(0);

hist << s << endl;   // On enregistre les commandes 

if (s[0] == '#') {
  cout << "PIACmd::Do() Comment-Line:" << s << endl;
  return(0); 
  } 
string toks,kw;
size_t p = s.find_first_not_of(" ");
s = s.substr(p);
p = 0;
size_t q = s.find_first_of(" ");
size_t l = s.length();

if (q < l)
  {  kw = s.substr(p,q-p);  toks = s.substr(q, l-q); }
else { kw = s.substr(p,l-p);  toks = ""; }

q = 0;
while (q < l)  {
  p = toks.find_first_not_of(" ",q+1); // au debut d'un token
  if (p>=l) break;
  q = toks.find_first_of(" ",p); // la fin du token;
  string token = toks.substr(p,q-p);
  tokens.push_back(token);
  }

for(int k=0; k<tokens.size(); k++) { // On remplace les $varname par la valeur de la variable
 if ((tokens[k])[0] != '$')  continue; 
 CmdVarList::iterator it = mVars.find(tokens[k].substr(1));
 if (it != mVars.end())  tokens[k] = (*it).second;
}

// cout << "PIACmd::Do() DBG  KeyW= " << kw << " NbArgs= " << tokens.size() << endl;
//  for(int ii=0; ii<tokens.size(); ii++)
//  cout << "arg[ " << ii << " ] : " << tokens[ii] << endl;

// >>>>>>>>>>> Commande d'interpreteur
if (kw == "help") {
  cout << "\n PIACmd::Do() ----- Help : Liste des commandes ------------ \n";
  cout << "help  shelp  exec  set  listvars  traceon  traceoff \n"; 
  cout << "timingon  timingoff  shell \n"; 
  cout << "link  call    zone  gratt  stacknext \n" ;
  cout << "openppf  saveall  listobjs openfits savefits \n"; 
  cout << "print rename  del delobjs \n";
  cout << "newh1d  newh2d  newprof newgfd func  func2d \n" ; 
  cout << "disp  surf  nt2d  nt3d  gfd2d  gfd3d \n" ;
  cout << "plot2d plot3d  projh1d  projh2d  projprof \n";  
  cout << "fillvec  fillnt fillgd1 fillgd2  fft fit fitpsf aper\n";
  cout << "*** Commande sans args -> liste d'arguments\n";
  cout << "*** dopt: DisplayOption = next/same/win/stack,color,... \n";
  cout << "*** $varname remplace par la valeur de la variable varname\n" << endl; 
  }

else if (kw == "shelp" )   {
  cout << "\n PIACmd::Do() ------------- SHelp ------------ " << endl;
  cout << "exec filename   --   set varname string  --  clr varname  --  listvars \n"; 
  cout << "traceon/off  --  timingon/off  --  shell cmdline\n";
  cout << "link fnameso f1 [f2 f3]  --   call userf [arg1 arg2 ...]\n";
  cout << "zone nx ny  --  gratt attributes_list  --  stacknext \n";
  cout << "openfits file -- openppf file -- saveall file -- listobjs --\n";
  cout << "print nomobj -- rename nomobj nomnew -- del nomobj -- delobjs pattern (*,?) \n";
  cout << "newh1d nomh xmin xmax nbin  --  newh2d nomh xmin xmax nbinx ymin ymax nbiny\n";
  cout << "newprof nomh xmin xmax nbin [ymin ymax] \n";
  cout << "newgfd nomgfd nvar nalloc [errx(0/1)]\n";
  cout << "func expfunc xmin xmax npt opt\n";
  cout << "func2d expfunc xmin xmax nptx ymin ymax npty opt\n";
  cout << "disp nomobj opt  --  surf nomobj opt\n";
  cout << "nt2d nomobj varx vary errx erry opt\n";
  cout << "nt3d nomobj varx vary varz errx erry errz dopt\n";
  cout << "gfd2d nomobj numvarx erreur=(x y xy) opt\n";
  cout << "gfd3d nomobj numvarx numvary erreur=(x y z xy xz yz xyz) opt\n";
  cout << "plot2d nomobj expx expy [experrx experry] expcut  dopt\n";
  cout << "plot3d nomobj expx expy expz expcut dopt\n";
  cout << "projh1d nomobj expx expwt expcut nomh1 dopt\n";
  cout << "projh2d nomobj expx expy expwt expcut nomh2 dopt\n";
  cout << "projprof nomobj expx expy expwt expcut nomprof dopt\n";
  cout << "fillnt nomobj expx expy expz expt  expcut nomnt dopt\n";
  cout << "fillvec nomobj expx expcut nomvec dopt\n";
  cout << "fillgd1 nomobj expx expy experry expcut nomgfd\n";
  cout << "fillgd2 nomobj expx expy expz experrz expcut nomgfd\n";
  cout << "fft nomobj dopt\n";
  cout << "fit nomobj func [p:p1,.,pn s:s1,.,sn m:m1,.,mn M:M1,.,Mn I:i1,i2 o:... o:...]\n";
  cout << "fitpsf psf_typ [options_fit]\n";
  cout << "aper fond_typ [options_calcul]\n";
}

else if (kw == "set") {
  if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: set varname string" << endl;  return(0); }
  mVars[tokens[0]] = tokens[1];
  }
else if (kw == "clr") {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: clr varname" << endl;  return(0); }
  CmdVarList::iterator it = mVars.find(tokens[0]);
  if (it != mVars.end())  mVars.erase(it);
  else cerr << "PIACmd::Do() Pas de variable de nom " << tokens[0] << endl;
  }
else if (kw == "listvars") {
  cout << "PIACmd::Do()  Variable List , VarName = Value \n";
  CmdVarList::iterator it;
  for(it = mVars.begin(); it != mVars.end(); it++)  
    cout << (*it).first << " = " <<  (*it).second << "\n";
  cout << endl;
  }
else if (kw == "traceon")  { cout << "PIACmd::Do()  -> Trace ON mode " << endl; trace = true; }
else if (kw == "traceoff") { cout << "PIACmd::Do()  -> Trace OFF mode " << endl; trace = false; }
else if (kw == "timingon") { 
  cout << "PIACmd::Do()  -> Timing ON mode " << endl; 
  if (gltimer)   delete gltimer;   gltimer = new Timer("PIA-CmdInterpreter ");   timing = true; }
else if (kw == "timingoff") { 
  cout << "PIACmd::Do()  -> Timing OFF mode " << endl; 
  if (gltimer)  delete gltimer;  gltimer = NULL;  timing = false; }
else if (kw == "exec") {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: exec filename" << endl;  return(0); }
  Exec(tokens[0]);
  }
else if (kw == "shell") {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: shell cmdline" << endl;  return(0); }
  system(toks.c_str());
  }

// >>>>>>>>>>> Fenetre graphique , changement d'attributs graphiques 
else if (kw == "zone") {
  if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: zone nx ny" << endl;  return(0); }
  int nx, ny;
  nx = ny = 1;
  nx = atoi(tokens[0].c_str());    ny = atoi(tokens[1].c_str());
  mObjMgr->SetGraphicWinZone(nx, ny);
  }  
else if (kw == "gratt") {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: gratt attributes_list (att=def->defaut)" << endl;  return(0); }
  mObjMgr->SetGraphicAttributes(tokens[0]); 
  }
else if (kw == "stacknext") mImgApp->StackWinNext(); 

// >>>>>>>>>>> Link dynamique de fonctions C++
else if (kw == "link" ) {
  if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: link fnameso f1 [f2 f3]" << endl;  return(0); }
  string sph = "";
  for(int gg=0; gg<5; gg++)   tokens.push_back(sph);
  int rc = LinkUserFuncs(tokens[0], tokens[1], tokens[2], tokens[3]);
  if (rc == 0)  cout << "PIAcmd::Do(): Link from " << tokens[0] << " OK " << endl;
}
else if (kw == "call" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: call userf [arg1 arg2 ...]" << endl;  return(0); }
  UsFmap::iterator it = usfmap.find(tokens[0]);
  if (it == usfmap.end()) {
    cerr << "PIAcmd::Do() No User Function " << tokens[0] << endl;
    return(0);
  }
  cout << "PIAcmd::Do() Call " << tokens[0] << "( ... )" << endl;
// on est oblige de faire un cast  etant donne qu'on 
// utilise donc des DlFunction  (Reza 20/08/98)  voir commentaire ds .h (pb g++)
  DlUserProcFunction fuf = (DlUserProcFunction)(*it).second;
  TRY {
    tokens.erase(tokens.begin());
    fuf(tokens);
  }  CATCH(merr) {
    fflush(stdout); 
    cout << endl; 
    cerr << endl;
    string es = PeidaExc(merr);
    cerr << "PIACmd::Do() Call UserFunc  Exception :" << merr << es;
    }
}

// >>>>>>>>>>> lecture/ecriture des objets, gestion des objets 
else if (kw == "openfits" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: openfits file " << endl;  return(0); }
  else mObjMgr->ReadFits(tokens[0]);
}
else if (kw == "savefits" ) {
  if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: savefits nomobj filename " << endl;  return(0); }
  else mObjMgr->SaveFits(tokens[0], tokens[1]);
}
else if (kw == "openppf" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: openppf file " << endl; return(0); }
  mObjMgr->ReadAll(tokens[0]);  
}
else if (kw == "saveall" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: saveall file " << endl; return(0); }
  mObjMgr->SaveAll(tokens[0]);  
}
else if (kw == "print" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: print nomobj " << endl; return(0); }
  mObjMgr->PrintObj(tokens[0]);  
}
else if (kw == "rename" ) {
  if (tokens.size() < 2) { cout << "PIACmd::Do() Usage: rename nomobj nomnew" << endl; return(0); }
  mObjMgr->RenameObj(tokens[0], tokens[1]);  
}
else if (kw == "del" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: del nomobj " << endl; return(0); }
  mObjMgr->DelObj(tokens[0]);  
}
else if (kw == "delobjs" ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: delobjs nomobjpattern (*,?) " << endl; return(0); }
  mObjMgr->DelObjects(tokens[0]);  
}
else if (kw == "listobjs")   mObjMgr->ListObjs();

// >>>>>>>>>>> Creation d'histos 1D-2D
else if (kw == "newh1d") {
  if (tokens.size() < 4) { cout << "PIACmd::Do() Usage: newh1d nom xmin xmax nbin" << endl; return(0); }
  int nbx;
  float xmin, xmax;
  nbx = 100;
  xmin = 0.;  xmax = 1.;
  nbx = atoi(tokens[3].c_str());
  xmin = atof(tokens[1].c_str());   xmax = atof(tokens[2].c_str());
  Histo* h = new Histo(xmin, xmax, nbx);
  mObjMgr->AddObj(h, tokens[0]);
  }
else if (kw == "newh2d") {
  if (tokens.size() < 7) { 
    cout << "PIACmd::Do() Usage:newh2d nomh xmin xmax nbinx ymin ymax nbiny" << endl; 
    return(0); 
    }
  int nbx, nby;
  float xmin, xmax;
  float ymin, ymax;
  nbx = nby = 50;
  xmin = 0.;  xmax = 1.;
  ymin = 0.;  ymax = 1.;
  nbx = atoi(tokens[3].c_str());
  nby = atoi(tokens[6].c_str());
  xmin = atof(tokens[1].c_str());   xmax = atof(tokens[2].c_str());
  ymin = atof(tokens[4].c_str());   ymax = atof(tokens[5].c_str());
  Histo2D* h = new Histo2D(xmin, xmax, nbx, ymin, ymax, nby);
  mObjMgr->AddObj(h, tokens[0]);
  }
else if (kw == "newprof") {
  if (tokens.size() < 4)
    { cout << "PIACmd::Do() Usage: newprof nom xmin xmax nbin [ymin ymax]" << endl; return(0); }
  int nbx;
  float xmin, xmax, ymin = 1., ymax = -1.;
  if(tokens.size() > 5)
    {ymin = atof(tokens[4].c_str());   ymax = atof(tokens[5].c_str());}
  nbx = 100;
  xmin = 0.;  xmax = 1.;
  nbx = atoi(tokens[3].c_str());
  xmin = atof(tokens[1].c_str());   xmax = atof(tokens[2].c_str());
  HProf* h = new HProf(xmin, xmax, nbx, ymin, ymax);
  mObjMgr->AddObj(h, tokens[0]);
  }
else if (kw == "newgfd") {
  if (tokens.size() < 3)
    { cout << "PIACmd::Do() Usage: newgfd nvar nalloc [errx(0/1)]" << endl; return(0); }
  int nvar, nalloc, errx=0;
  if (tokens.size() > 3)
    { errx = atoi(tokens[3].c_str()); if(errx>0) errx=1; else errx = 0;}
  nvar = atoi(tokens[1].c_str()); nalloc = atoi(tokens[2].c_str());
  if(nvar>0 && nalloc>0) {
    GeneralFitData* gfd = new GeneralFitData(nvar,nalloc,errx);
    mObjMgr->AddObj(gfd, tokens[0]);
    }
  }

// >>>>>>>>>>>  Affichage des objets 
else if ( (kw == "disp") || (kw == "surf") ) {
  if (tokens.size() < 1) { cout << "PIACmd::Do() Usage: disp nomobj opt" << endl; return(0); }
  string opt = "n";
  if (tokens.size() > 1)  opt = tokens[1];
  if (kw == "disp") mObjMgr->DisplayObj(tokens[0], opt);
  else if (kw == "surf")  mObjMgr->DisplaySurf3D(tokens[0], opt);
  }

else if (kw == "nt2d") {
  if (tokens.size() < 5) { cout << "PIACmd::Do() Usage: nt2d nomobj varx vary errx erry opt" << endl; return(0); }
  string opt = "n";
  if (tokens.size() > 5)  opt = tokens[5];
  string ph = "";
  mObjMgr->DisplayNT(tokens[0],tokens[1],tokens[2], ph, tokens[3], tokens[4], ph, opt);
  }
else if (kw == "nt3d") {
  if (tokens.size() < 7) { cout << "PIACmd::Do() Usage: nt3d nomobj varx vary varz errx erry errz opt" << endl; return(0); }
  string opt = "n";
  if (tokens.size() > 7)  opt = tokens[7];
  mObjMgr->DisplayNT(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], tokens[6], opt);
  }

else if (kw == "gfd2d") {
  if(tokens.size()<2)
    {cout<<"PIACmd::Do() Usage: gfd2d nomobj numvarx erreur=(x y xy) opt"<<endl;
     return(0);}
  string numvary = "";
  string err = "";
  string opt = "n";
  if(tokens.size()>2) err = tokens[2];
  if(tokens.size()>3) opt = tokens[3];
  mObjMgr->DisplayGFD(tokens[0],tokens[1],numvary,err,opt);
  }
else if (kw == "gfd3d") {
  if(tokens.size()<3)
    {cout<<"PIACmd::Do() Usage: gfd3d nomobj numvarx numvary erreur=(x y z xy xz yz xyz) opt"<<endl;
     return(0);}
  string err = "";
  string opt = "n";
  if(tokens.size()>3) err = tokens[3];
  if(tokens.size()>4) opt = tokens[4];
  mObjMgr->DisplayGFD(tokens[0],tokens[1],tokens[2],err,opt);
  }

// >>>>>>>>>>>  Trace de fonctions 
else if ( (kw == "func") ) {
  if (tokens.size() < 4) { cout << "PIACmd::Do() Usage: func expr xmin xmax npt opt" << endl; return(0); }
  string opt = "n";
  if (tokens.size() > 4)  opt = tokens[4];
  int np;
  float xmin, xmax;
  np = 100;
  xmin = 0.;  xmax = 1.;
  np = atoi(tokens[3].c_str());
  xmin = atof(tokens[1].c_str());   xmax = atof(tokens[2].c_str());
  mObjMgr->PlotFunc(tokens[0], xmin, xmax, np, opt); 
  }
else if ( (kw == "func2d") ) {
  if (tokens.size() < 7) {
    cout << "PIACmd::Do() Usage: func2d expr xmin xmax nptx ymin ymax npty opt" << endl;
    return(0);
    }
  int npx, npy;
  float xmin, xmax;
  float ymin, ymax;
  npx = npy = 50;
  xmin = 0.;  xmax = 1.;
  ymin = 0.;  ymax = 1.;
  npx = atoi(tokens[3].c_str());
  npy = atoi(tokens[6].c_str());
  xmin = atof(tokens[1].c_str());   xmax = atof(tokens[2].c_str());
  ymin = atof(tokens[4].c_str());   ymax = atof(tokens[5].c_str());
  string opt = "n";
  if (tokens.size() > 7)  opt = tokens[7];
  mObjMgr->PlotFunc2D(tokens[0], xmin, xmax, ymin, ymax, npx, npy, opt);
  }

// >>>>>>>>>>>  Trace d'expressions de N_Tuple, StarList, etc ...
else if (kw == "plot2d" ) {
  if (tokens.size() < 4) {
    cout << "PIACmd::Do() Usage: plot2d nomobj expx expy [experrx experry] expcut opt" << endl;
    return(0); 
    }
  string errx = ""; string erry = ""; string ecut = "1";
  string opt = "n";
  if (tokens.size() < 6) { // Plot sans les erreurs
    ecut = tokens[3];
    if (tokens.size() > 4)  opt = tokens[4];
    }
  else {                 // Plot avec les erreurs
    errx = tokens[3]; erry = tokens[4]; ecut = tokens[5];
    if (tokens.size() > 6)  opt = tokens[6];
    }
  mObjMgr->DisplayPoints2D(tokens[0],tokens[1],tokens[2],errx,erry,ecut,opt);
  }

else if (kw == "plot3d" ) {
  if (tokens.size() < 5) { 
    cout << "PIACmd::Do() Usage: plot3d nomobj expx expy expz expcut opt" << endl;
    return(0); 
    }
  string opt = "n";
  if (tokens.size() > 5)  opt = tokens[5];
  mObjMgr->DisplayPoints3D(tokens[0],tokens[1],tokens[2],tokens[3], tokens[4], opt);
  }

else if (kw == "projh1d" ) {
  if (tokens.size() < 5) { 
    cout << "PIACmd::Do() Usage: projh1d nomobj expx expwt expcut nomh1 opt" << endl;
    return(0);
    }
  string opt = "n";
  if (tokens.size() > 5)  opt = tokens[5];
  mObjMgr->ProjectH1(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], opt);
  }

else if (kw == "projh2d" ) {
  if (tokens.size() < 6) { 
    cout << "PIACmd::Do() Usage: projh2 nomobj expx expy expwt expcut nomh2 opt" << endl;
    return(0);
    }
  string opt = "n";
  if (tokens.size() > 6)  opt = tokens[6];

  mObjMgr->ProjectH2(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], opt);
  }

else if (kw == "projprof" ) {
  if (tokens.size() < 6) { 
    cout << "PIACmd::Do() Usage: projprof nomobj expx expy expwt expcut nomprof opt" << endl;
    return(0);
    }
  string opt = "n";
  if (tokens.size() > 6)  opt = tokens[6];

  mObjMgr->ProjectHProf(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], opt);
  }

else if (kw == "fillnt" ) {
  if (tokens.size() < 7) { 
    cout << "PIACmd::Do() Usage: fillnt nomobj expx expy expz expt expcut nomnt" << endl;
    return(0);
    }
  mObjMgr->FillNT(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5], tokens[6] );
  }

else if (kw == "fillvec" ) {
  if (tokens.size() < 4) { 
    cout << "PIACmd::Do() Usage: fillvec nomobj expx expcut nomvec opt" << endl;
    return(0);
    }
  string opt = "n";
  if (tokens.size() > 4)  opt = tokens[4];
  mObjMgr->FillVect(tokens[0],tokens[1],tokens[2], tokens[3], opt);
  }

else if (kw == "fillgd1" ) {
  if (tokens.size() < 5) {
    cout << "PIACmd::Do() Usage: fillgd1 nomobj expx expy experry expcut nomgfd" << endl;
    return(0);
    }
  string nomgfd = "";
  if (tokens.size() > 5) nomgfd = tokens[5];
  string expy = "";
  mObjMgr->FillGFD(tokens[0],tokens[1], expy, tokens[2], tokens[3], tokens[4],nomgfd);
  }

else if (kw == "fillgd2" ) {
  if (tokens.size() < 6) {
    cout << "PIACmd::Do() Usage: fillgd2 nomobj expx expy expz experrz expcut nomgfd" << endl;
    return(0);
    }
  string nomgfd = "";
  if (tokens.size() > 6) nomgfd = tokens[6];
  mObjMgr->FillGFD(tokens[0],tokens[1],tokens[2], tokens[3], tokens[4], tokens[5],nomgfd);
  }

// Calcul TF 
else if (kw == "fft" ) {
  if (tokens.size() < 1) { 
    cout << "PIACmd::Do() Usage: fft nomobj opt" << endl;
    return(0);
    }
  string opt = "n";
  if (tokens.size() > 1)  opt = tokens[1];
  mObjMgr->FFT(tokens[0], opt);
  }

// Fit 1D sur objets 1D. Egalement Fit 2D sur objets 2D.
else if (kw == "fit") {
  if (tokens.size() < 2) {
    cout <<"PIACmd::Do() Usage:fit nomobj func \n"
         <<" [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]\n";
    return(0);
  }
  string p=""; string s=""; string m=""; string M=""; string O="";
  if (tokens.size()>2)
    for(int ip=2;ip<tokens.size();ip++) {
      if(tokens[ip].length()<=2) continue;
      const char *c = tokens[ip].c_str();
      if(c[1]!=':') continue;
      if(c[0]=='p')      p=c+2;
      else if(c[0]=='s') s=c+2;
      else if(c[0]=='m') m=c+2;
      else if(c[0]=='M') M=c+2;
      else if(c[0]=='o') {O += ","; O += c+2;}
    }
  mObjMgr->Fit12D(tokens[0],tokens[1],p,s,m,M,O);
}

#ifndef SANS_StarReco 
// FITPSF sur image
else if (kw == "fitpsf" ) {
  if (tokens.size() < 1) { 
    cout << "PIACmd::Do() Usage: fitpsf psf_typ [options_fit]" << endl;
    return(0);
    }
  FitPSFImg(tokens);
  }

// Photometrie d'ouvertue sur image 
else if (kw == "aper" ) {
  if (tokens.size() < 1) { 
    cout << "PIACmd::Do() Usage: aper fond_typ [opt_calcul opt_disp]" << endl;
    return(0);
    }
  Photom_Ouv(tokens);
  }
#endif

else  { 
  cerr << "PIACmd::Do() Erreur - Commande " << kw << " inconuue ! " << endl;    
  return(-1);
  }
if (timing)  gltimer->Split();

return(0);
}



/* --Methode-- */
int PIACmd::LinkUserFuncs(string& fnameso, string& func1, string& func2, string& func3)
//                          string& func4, string& func5)
{
string cmd;
int rc;

if (dynlink) delete dynlink;    dynlink = NULL;
usfmap.clear();

dynlink = new PDynLinkMgr(fnameso, true);
if (dynlink == NULL) { 
  string sn = fnameso; 
  cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur ouverture SO " << sn << endl; 
  return(2); 
  }

int nok=0;
// on utilise donc des DlFunction  (Reza 20/08/98)  voir commentaire ds .h (pb g++)
// DlUserProcFunction f = NULL;
DlFunction f = NULL;
if ((func1.length() < 1) || (func1 == "-") || (func1 == ".") )  goto fin;
// f = (DlUserProcFunction)  dlsym(dlhandle, func1.c_str());
f = dynlink->GetFunction(func1);
if (f) { nok++;  usfmap[func1] = f; }
else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func1 << endl;

if ((func2.length() < 1) || (func2 == "-") || (func2 == ".") )  goto fin;
// f = (DlUserProcFunction)  dlsym(dlhandle, func2.c_str());
f = dynlink->GetFunction(func2);
if (f) { nok++;  usfmap[func2] = f; }
else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func2 << endl;

if ((func3.length() < 1) || (func3 == "-") || (func3 == ".") )  goto fin;
// f = (DlUserProcFunction)  dlsym(dlhandle, func3.c_str());
f = dynlink->GetFunction(func3);
if (f) { nok++;  usfmap[func3] = f; }
else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func3 << endl;

/*  Pb compile g++ 2.7.2 
if ((func4.length() < 1) || (func4 == "-") || (func4 == ".") )  goto fin;
// f = (DlUserProcFunction)  dlsym(dlhandle, func4.c_str());
f = dynlink->GetFunction(func4);
if (f) { nok++;  usfmap[func4] = f; }
else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func4 << endl;

if ((func5.length() < 1) || (func5 == "-") || (func5 == ".") )  goto fin;
// f = (DlUserProcFunction)  dlsym(dlhandle, func5.c_str());
f = dynlink->GetFunction(func5);
if (f) { nok++;  usfmap[func5] = f; }
else cerr << "PIACmd/LinkUserFuncs_Erreur: Erreur linking " << func5 << endl;
*/
fin:
if (nok < 1) { if (dynlink) delete dynlink;    dynlink = NULL;  return(3); }
else return(0);
}

/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////// Fonctions ///////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
#ifndef SANS_StarReco 
#include "psf.h"
#include "dynccd.h"
#include "pimgadapter.h"
#include "nbmath.h"
#include "nbtri.h"

/* Nouvelle-Fonction    rz + cmv 15/10/98 */
int FitPSFImg(cmdtok& tokens)
//      Fit de PSF sur une Image.
//| tokens[0] = Type de PSF ()
//|      : G   = fit gaussienne+fond (volume)
//|      : Gi  = fit gaussienne+fond integree (volume)
//|      : d   = fit DL de gaussienne+fond (volume)
//|      : di  = fit DL de gaussienne+fond integree (volume)
//|      : D   = fit DL de gaussienne+fond avec coeff variable p6 (volume)
//|      : Di  = fit DL de gaussienne+fond integree avec coeff variable p6 (volume)
//|      : M   = fit Moffat+fond (expos=p6) (volume)
//|      : Mi  = fit Moffat+fond integree (expos=p6) (volume)
//| tokens[1] = options pour le fit ex: l1.0,G2.,R7.,m10.,M65000.,P3.,f,n,r,i,w10/450
//|      : r = ecriture de l'image des residus (def=non)
//|      : f = ecriture de l'image de de la PSF fitee (def=non)
//|      : i = ecriture du pave initial selectionne pour le fit (def=non)
//|      : n = ecriture de l'image du bruit (def=non)
//|      : Gaa.bbb = gain en e-/ADU (def=1)
//|      : Raa.bbb = readout noise en e- (def=0)
//|      : maa.bbb = limite basse dynamique en ADU (def=-9.e19)
//|      : Maa.bbb = limite haute dynamique en ADU (def=9.e19)
//|      : Paa.bbb = valeur du parametre supplementaire pour 
//|      : la.a = print a=FitPSFImg, b=PSF.Fit PSF.IniFit
//|      : wa/b = "a" est le facteur de zoom de la fenetre de representation
//|               des resultats (def=6) et b sa taille maximale (def=500).
//| Return: 0 si Ok, <0 sinon.
{
// Decodage des arguments
if(tokens.size()<1) return(-1);
int typ_psf=PSF::psfGaussienne;
bool wres=false, wpsf=false, wpav=false, wnois=false;
float gain=1., readout=0., pixmin=-9.e19, pixmax=9.e19;
int zoom=6, wmax=500;
float par6=-1.;
int lp=1, lpg=0;
string dummy,dum;

if(tokens.size()>=2) {
  dummy = "," + tokens[1]+ ",";
  if(strstr(dummy.c_str(),",r,")) wres=true;
  if(strstr(dummy.c_str(),",f,")) wpsf=true;
  if(strstr(dummy.c_str(),",i"))  wpav=true;
  if(strstr(dummy.c_str(),",n"))  wnois=true;
  if(strstr(dummy.c_str(),",G")) {
    size_t p = dummy.find(",G"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",G%f",&gain);
  }
  if(strstr(dummy.c_str(),",P")) {
    size_t p = dummy.find(",P"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",P%f",&par6);
  }
  if(strstr(dummy.c_str(),",R")) {
    size_t p = dummy.find(",R"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",R%f",&readout);
  }
  if(strstr(dummy.c_str(),",m")) {
    size_t p = dummy.find(",m"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",m%f",&pixmin);
  }
  if(strstr(dummy.c_str(),",M")) {
    size_t p = dummy.find(",M"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",M%f",&pixmax);
  }
  if(strstr(dummy.c_str(),",l")) { // niveau de print
    size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    float ab;
    if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
    if(ab<0) ab = 0.;
    lp = (int) ab; lpg = (int) 10.*(ab-lp);
  }
  if(strstr(dummy.c_str(),",w")) { // niveau de print
    size_t p = dummy.find(",w"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",w%d/%d",&zoom,&wmax);
  }
}

dummy = "," + tokens[0]+ ",";
if(strstr(dummy.c_str(),",G,"))        typ_psf=PSF::psfGaussienne;
else if(strstr(dummy.c_str(),",Gi,"))  typ_psf=PSF::psfGaussInteg;
else if(strstr(dummy.c_str(),",d,") )  typ_psf=PSF::psfDelGaussienne;
else if(strstr(dummy.c_str(),",di,"))  typ_psf=PSF::psfDelGaussinteg;
else if(strstr(dummy.c_str(),",D,") ) {typ_psf=PSF::psfDel1Gaussienne; if(par6<=0.) par6=1.;}
else if(strstr(dummy.c_str(),",Di,")) {typ_psf=PSF::psfDel1GaussInteg; if(par6<=0.) par6=1.;}
else if(strstr(dummy.c_str(),",M,") ) {typ_psf=PSF::psfMoffat;         if(par6<=1.) par6=2.5;}
else if(strstr(dummy.c_str(),",Mi,")) {typ_psf=PSF::psfMoffInteg;      if(par6<=1.) par6=2.5;}

if(gain<=0.) gain=1.; if(readout<0.) readout=0.;
if(zoom<1) zoom=1; if(zoom>10) zoom=10; // Il semble que le zoom soit limite a 10 !!!
if(wmax>500) wmax=500;
if(lp<0) lp=1; if(lpg<0) lpg=0;

// Wigdet,piimage et image courants
NamedObjMgr omg;
PIStdImgApp* app = omg.GetImgApp();
if(app == NULL)
  {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
PIBaseWdg* curw= app->CurrentBaseWdg() ;
if(curw == NULL)
  {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
if(curw->kind() != PIImage::ClassId)
  {cout<<"Current window not an Image"<<endl; return(-4);}
PIImage* curpimg = (PIImage*)curw;
if(curpimg == NULL)
  {cout<<"Pas de PIImage*"<<endl; return(-5);}
P2DArrayAdapter* img = curpimg->Image();
if(img == NULL)
  {cout<<"Pas de Image*"<<endl; return(-6);}
string nom = "fitpsf_";

// Definition du pave courant sur la piimage
int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();  
int x0 = curpimg->XPave()-dx/2, y0 = curpimg->YPave()-dy/2;
if(x0<0) x0 = 0; if(y0<0) y0 = 0;
int x1 = x0+dx, y1 = y0+dy;
if(x1>img->XSize()) x1=img->XSize(); if(y1>img->YSize()) y1=img->YSize();
dx = x1-x0;  dy = y1-y0;
int dxy = (dx>dy)?dx:dy;
ImageR4* pav = new ImageR4(dx, dy); 
{for(int j=0;j<dy;j++) for(int i=0;i<dx;i++) (*pav)(i,j)=(*img)(i+x0,j+y0);}

// Bruit sur le pave
DynCCD mdynccd(kPhotonNoise,pixmin,pixmax,gain,readout);
ImageR4 *nois = NoiseImage(pav, &mdynccd);

// Print
int w=2*zoom*dxy+25; if(w>wmax){w=wmax;zoom=(wmax-25)/(2*dxy);if(zoom<1)zoom=1;}
if(lp>1)
  {cout<<"FitPSFImg: fun="<<tokens[0]<<","<<typ_psf
       <<"  Write res="<<wres<<" psf="<<wpsf<<" pav="<<wpav<<" nois="<<wnois
       <<"  Par6="<<par6<<"  lp="<<lp<<","<<lpg<<endl;
   mdynccd.Print();
   cout<<"... dx,y="<<dxy<<" zoom="<<zoom<<" w="<<w<<endl;}

// PSF et Fit
PSF mypsf(typ_psf,100.0);
int npar = mypsf.NParm();
PSFStar mypsfstar(&mypsf);
if(npar>7) mypsfstar.Parm(6)=(double)par6; // Cas moffat et Dl gauss avec coeff var.
int rc = mypsfstar.IniFit(*pav,*nois,lpg);
if(rc != 0)
  {cout<<"PIDemoApp::FitPSFImg(), Echec IniFit (rc="<<rc<<")"<<endl;
   delete nois;  delete pav; return(-10); }
rc = mypsfstar.Fit(*pav,*nois,lpg);
if(rc != 0)
  {cout<<"PIDemoApp::FitPSFImg(), Echec Fit (rc="<<rc<<")"<<endl;
   delete nois;  delete pav; return(-11); }
if(lp) mypsfstar.Print(2);
float smax,axisrat,tiltdeg;
rc = paramga((float) mypsfstar.SigX(),(float) mypsfstar.SigY()
            ,(float) mypsfstar.Rho(),&smax,&axisrat,&tiltdeg);
if(rc==0 && lp)
  printf("...   sa=%g sb=%g (%g) teta=%g\n",smax,smax*axisrat,axisrat,tiltdeg);
else printf("Erreur paramga \n");

// Creation des pave resultats
ImageR4* res = new ImageR4(dx, dy);
ImageR4* psf = new ImageR4(dx, dy);
CopieImage(*res,*pav); *res -= mypsfstar;
*psf = 0.; *psf += mypsfstar;

// Si on demande d'ecrire au moins un des resultats (pav,nois,res,psf) => menage!
bool delobj = (wpav||wnois||wres||wpsf)?true:false;

// Display des resultats, fenetre 2x2 {{pav,nois},{res,psf}}
app->CreateGraphWin(2,2,w,w); app->SetZoomAtt(zoom);
dummy = nom+"pav"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
if(wpav) {
  // Si on veut garder l'objet
  omg.AddObj(pav,dummy);
  omg.DisplayObj(omg.LastObjName());
} else {
  // "true" pour tuer (delete) de l'objet ImageR4 lors de la destruction de PIImage
  app->DispImage(new ImageAdapter<r_4>(pav,true),dummy,Disp_Next);
}
dummy = nom+"nois"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
if(wnois) {
  omg.AddObj(nois,dummy);
  omg.DisplayObj(omg.LastObjName());
} else {
  app->DispImage(new ImageAdapter<r_4>(nois,true),dummy,Disp_Next);
}
dummy = nom+"res"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
if(wres) {
  omg.AddObj(res,dummy);
  omg.DisplayObj(omg.LastObjName());
} else {
  app->DispImage(new ImageAdapter<r_4>(res,true),dummy,Disp_Next);
}
dummy = nom+"psf"; dum = dummy+"*"; if(delobj) omg.DelObjects(dum);
if(wpsf) {
  omg.AddObj(psf,dummy);
  omg.DisplayObj(omg.LastObjName());
} else {
  app->DispImage(new ImageAdapter<r_4>(psf,true),dummy,Disp_Next);
}

return(0);
}

/* Nouvelle-Fonction    cmv 16/10/98 */
int Photom_Ouv(cmdtok& tokens)
//      Photometrie d'ouverture sur une Image.
//| tokens[0] = Type de calcul du fond et du flux.
//|   : b = fond calcule sur les pixels du bord du pave.
//|   : a = fond sur les pixels autres que ceux servant au calcul du flux.
//|   : faa.bbb = fond fixe a aa.bbb
//|   defaut : "b"
//| tokens[1] = options pour le fit ex: l1
//|      : ra = "a" rayon minimal de scan (def=0)
//|      : ga = "a" taille de la fenetre glissante pour les extrapolations (def=4)
//|      : la = "a" niveau de print
//|      : c  = la surface pour la photometrie est un carre (defaut=cercle)
//|      : C  = recentrage du pave sur le pixel le + brillant  (defaut=non)
//| tokens[2] = options graphiques
//| Return: 0 si Ok, <0 sinon.
{
// Decodage des arguments
if(tokens.size()<1) return(-1);
int typ_fond=1, rmin=0, rmax, npt_gliss=4, lp=0;
float fondimpos=0., fondbord=0., fondautre;
bool cercle=true, centre=false;
string opt="";
string dummy; string dum;

if(tokens.size()>=3) opt = tokens[2];
if(tokens.size()>=2) {
  dummy = "," + tokens[1]+ ",";
  if(strstr(dummy.c_str(),",r")) {
    size_t p = dummy.find(",r"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",r%d",&rmin);
  }
  if(strstr(dummy.c_str(),",g")) {
    size_t p = dummy.find(",g"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",g%d",&npt_gliss);
  }
  if(strstr(dummy.c_str(),",l")) {
    size_t p = dummy.find(",l"); size_t q = dummy.find_first_of(',',p+1);
    dum = dummy.substr(p,q-p);
    if(dum.length()>2) sscanf(dum.c_str(),",l%d",&lp);
  }
  if(strstr(dummy.c_str(),",c")) cercle = false;
  if(strstr(dummy.c_str(),",C")) centre = true;
}

dummy = "," + tokens[0]+ ",";
if(strstr(dummy.c_str(),",b,"))       typ_fond=1;
else if(strstr(dummy.c_str(),",a,"))  typ_fond=2;
else if(strstr(dummy.c_str(),",f"))  {
  typ_fond=3;
  size_t p = dummy.find(",f"); size_t q = dummy.find_first_of(',',p+1);
  dum = dummy.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",f%f",&fondimpos);
}

if(rmin<0) rmin=0; if(npt_gliss<=0) npt_gliss=4; if(lp<0) lp=1;

// Wigdet,piimage et image courants
NamedObjMgr omg;
PIStdImgApp* app = omg.GetImgApp();
if(app == NULL)
  {cout<<"Pas de PIStdImgApp*"<<endl; return(-2);}
PIBaseWdg* curw= app->CurrentBaseWdg() ;
if(curw == NULL)
  {cout<<"Pas de PIBaseWdg*"<<endl; return(-3);}
if(curw->kind() != PIImage::ClassId)
  {cout<<"Current window not an Image"<<endl; return(-4);}
PIImage* curpimg = (PIImage*)curw;
if(curpimg == NULL)
  {cout<<"Pas de PIImage*"<<endl; return(-5);}
P2DArrayAdapter* img = curpimg->Image();
if(img == NULL)
  {cout<<"Pas de Image*"<<endl; return(-6);}
string nom = "ouver_";

// Definition du pave courant sur la piimage
int xc = curpimg->XPave(), yc = curpimg->YPave();
int dx = curpimg->XSzPave(), dy = curpimg->YSzPave();
int dxym = (dx<dy)?dx:dy; if(dxym%2==0) dxym++; rmax = dxym/2;
if(dxym<3||dxym>img->XSize()||dxym>img->YSize())
  {cout<<"Pave trop petit/grand: "<<dxym<<endl; return(-10);}
if(rmin>=rmax)
  {cout<<"rayon mini="<<rmin<<">= rmax="<<rmax<<endl; return(-11);}
if(centre) {  // Recentrage autour du pixel le + brillant?
  int xcnew=-1, ycnew=-1; float pxnew=-1.e20,pxold=-1.e20,px;
  for(int j=yc-rmax;j<=yc+rmax;j++) {
    if(j<0 || j>=img->YSize()) continue;
    for(int i=xc-rmax;i<=xc+rmax;i++) {
      if(i<0 || i>=img->XSize()) continue;
      px = (*img)(i,j);
      if(i==xc && j==yc) pxold = px;
      if((xcnew==-1 && ycnew==-1) || px>=pxnew) {xcnew=i; ycnew=j; pxnew=px;}
      // cout<<"xcnew,ycnew,pxnew="<<xcnew<<" "<<ycnew<<" "<<pxnew<<endl;
  } }
  if(lp>1) cout<<"Re-centrage(x,y,px): old=("<<xc<<","<<yc<<","<<pxold<<")"
               <<" -> new=("<<xcnew<<","<<ycnew<<","<<pxnew<<")"<<endl;
  if(xcnew!=-1 && ycnew!=-1 && pxnew>pxold) {xc=xcnew; yc=ycnew;}
}
int x0=xc-rmax, y0=yc-rmax, x1=xc+rmax, y1=yc+rmax;
if(x0<0||y0<0||x1>=img->XSize()||y1>=img->YSize())
  {cout<<"Pave hors image"<<endl; return(-12);}
int ic = rmax, jc = rmax;
ImageR4* pav = new ImageR4(dxym, dxym);
{for(int j=-rmax;j<=rmax;j++) for(int i=-rmax;i<=rmax;i++)
    (*pav)(ic+i,jc+j)=(*img)(xc+i,yc+j);}

// Print
if(lp>1)
  cout<<"photom_ouv[cercle="<<cercle<<",typ="<<typ_fond<<"]: pav["
      <<dx<<","<<dy<<","<<dxym
      <<"] xc,yc=["<<xc<<","<<yc<<"] ic,jc="<<ic<<","<<jc<<endl
      <<"          (fond="<<fondimpos<<") pxmil="<<(*pav)(ic,ic)
      <<" -> Rph de "<<rmin<<" a "<<rmax<<endl;

// Ouverture des buffers et objets necessaires
float *p = new float[dxym*dxym];
r_4 xnt[4];
NTuple *ntfl=NULL, *ntex=NULL;
{char* name[4]={"npx","flux","fond","mag"}; ntfl = new NTuple(4,name);}
{char* name[3]={"k","a","f0"};              ntex = new NTuple(3,name);}

// Calcul du fond sur les pixels du bord du pave
{
int nf = 0;
for(int i=0;i<dxym;i++)   {p[nf]=(*pav)(i,0); nf++; p[nf]=(*pav)(i,dxym-1); nf++;}
for(int j=1;j<dxym-1;j++) {p[nf]=(*pav)(0,j); nf++; p[nf]=(*pav)(dxym-1,j); nf++;}
qsort(p,(size_t) nf,sizeof(float),qSort_Float);
fondbord = (p[nf/2]+p[(nf-1)/2])/2.;
if(lp>0) cout<<"Fond_Bord: "<<fondbord<<" (avec 1 px au bord="<<nf<<")"<<endl;
}

// Photometrie d'ouverture variable.
int rfin = -1;
for(int r=rmin;r<=rmax;r++) {
  double sum=0.; int ns=0, nf=0;
  for(int i=0;i<dxym;i++) for(int j=0;j<dxym;j++) {
    double ray2 = (i-ic)*(i-ic)+(j-jc)*(j-jc)-0.1;
    if( (cercle && ray2<=r*r) ||
        (!cercle && ic-r<=i && i<=ic+r && jc-r<=j && j<=jc+r) )
          {sum += (*pav)(i,j); ns++;}
    else {p[nf] = (*pav)(i,j); nf++;}
  }
  if(typ_fond==2) {
    if(nf<2) continue;
    qsort(p,(size_t) nf,sizeof(float),qSort_Float);
    fondautre = (p[nf/2]+p[(nf-1)/2])/2.;
  }
  double fond=-1.e19;
  if(typ_fond==3) fond=fondimpos; else if(typ_fond==2) fond=fondautre; else fond=fondbord;
  rfin = r;
  xnt[0] = ns;
  xnt[1] = sum-fond*ns;
  xnt[2] = fond;
  xnt[3] = (xnt[1]>0.)? -2.5*log10((double)xnt[1]): -99.;
  if(lp>0)
    cout<<"Ouv[r="<<r<<",ns="<<xnt[0]<<",nf="<<nf<<"]: fd="<<xnt[2]
        <<" fl="<<xnt[1]<<" mag="<<xnt[3]<<endl;
  ntfl->Fill(xnt);
}
delete [] p;
int nr = ntfl->NEntry();
if(lp>1) cout<<"Nombre de calculs de flux: "<<nr<<endl;
if(nr<=0)
  {cout<<"Pas de flux calcules"<<endl; delete pav; delete ntfl; delete ntex; return(-13);}

cout<<"Ouv[rouv="<<rfin<<",npix="<<xnt[0]<<"]: fond="<<xnt[2]
    <<" flux="<<xnt[1]<<" mag="<<xnt[3]<<endl;

// Extrapolation: F_mes = F_vrai + npxflux*(Fd_vrai-Fd_mes)
// On fit F_mes=f(npxflux) sur une fenetre glissante de npt_gliss pts.
bool ex = false;
if(nr>=npt_gliss) {
  ex = true;
  double *npix = new double[nr], *flux = new double[nr], *err  = new double[nr];
  for(int i=0;i<nr;i++) { ntfl->GetVec(i,xnt); npix[i]=xnt[0]; flux[i]=xnt[1]; err[i]=1.;}
  double a,b; int n;
  for(int k=0;k<=nr-npt_gliss;k++) {
    n = npt_gliss; FitLin(&npix[k],&flux[k],&err[k],&n,&b,&a);
    xnt[0]=k; xnt[1]=a; xnt[2]=b; ntex->Fill(xnt);
    if(lp>1)
      cout<<"Fit["<<k<<","<<k+npt_gliss-1<<"] f0="<<xnt[2]<<" pente="<<xnt[1]<<endl;
  }
  delete [] npix; delete [] flux; delete [] err;
}

// Display des resultats, fenetre 2x2 {{flux%ns,fond%ns},{f0,pente}}
if(ex) app->CreateGraphWin(2,2); else app->CreateGraphWin(1,2);
dummy = nom+"*"; omg.DelObjects(dummy);

dummy=nom+"pav"; omg.AddObj(pav,dummy);

string xd=""; string yd=""; string zd,exd,eyd,ezd;
dummy=nom+"flux"; omg.AddObj(ntfl,dummy);
xd="npx"; yd="flux";
omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
if(typ_fond==2)
  {xd="npx"; yd="fond"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
if(ex) {
  dummy=nom+"fit"; omg.AddObj(ntex,dummy);
  if(typ_fond!=2)
    {xd="a"; yd="f0"; omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);}
  xd = "k"; yd = "f0";
  omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
  xd = "k"; yd = "a";
  omg.DisplayNT(omg.LastObjName(),xd,yd,zd,exd,eyd,ezd,opt);
}

return(0);
}
#endif