source: Sophya/trunk/SophyaPI/ProgPI/sopiamodule.cc@ 3430

#include "sopnamsp.h"
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <ctype.h>
#include <math.h>

#include <typeinfo>

#include <vector>
#include <string>

#include "piacmd.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"
#include "servnobjm.h"
#include "tvector.h"
#include "pitvmaad.h"
#include "fftpserver.h"
#include "bruit.h"
#include "piscdrawwdg.h"
#include "ctimer.h"

#include "pidrawer.h"
#include "nomgadapter.h"
#include "nomskymapadapter.h"

extern "C" {
void sopiamodule_init();
void sopiamodule_end();
}

// --- fonctions pour faire FFT
void SophyaFFT_forw(string& nom, string& nomout, string dopt);
void SophyaFFT_back(string& nom, string& nomout, string dopt, bool forcez8=false);
// void SophyaFFT_crefilter(string& nomvec, string& nomfilter, string func, string dopt);
void SophyaFFT_filter(string& nom, string& filt, string& outvec, string dopt, bool fgvec=false);

//--- fonctions pour trace d'ellipse d'erreur
int  diag_ellipse(double s1, double s2, double c, 
                  double &a, double &b, double& tet, double& phi, int prt=0);

void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy, 
                 double scale=1., int nbpoints=101);

//--- Declaration du CmdExecutor du module
class sopiamoduleExecutor : public CmdExecutor {
public:
                sopiamoduleExecutor();
  virtual       ~sopiamoduleExecutor();
  virtual int   Execute(string& keyw, vector<string>& args, string& toks);
};

//---  Classe pour trace de grille en projection spherique Mollweide 
class MollweideGridDrawer : public PIDrawer {
public:
                    MollweideGridDrawer(int np, int nm, int nx, int ny);
                    MollweideGridDrawer(int np, int nm, double x_larg=0., double x_0=0., 
                                        double y_larg=0., double y_0=0.);
  virtual           ~MollweideGridDrawer();
  virtual void      Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax);
  virtual void      UpdateLimits();

protected:
  double xlarg, x0;
  double ylarg, y0;
  int nbparall, nbmerid;
  int nbpt;
  bool fgrevtheta;
};

/* --Methode-- */
MollweideGridDrawer::MollweideGridDrawer(int np, int nm, int nx, int ny)
{
  x0 = (double)nx/2.; 
  y0 = (double)ny/2.;
  xlarg = nx;  
  ylarg = ny;
  nbparall = np;
  nbmerid = nm;
  nbpt = ny/5;
  if ((nbpt % 20) == 0) nbpt++;  
  fgrevtheta = true;
}

/* --Methode-- */
MollweideGridDrawer::MollweideGridDrawer(int np, int nm, double x_larg, double x_0, 
                                         double y_larg, double y_0)
{
  if (x_larg > 0.) {
    xlarg = x_larg; x0 = x_0;
  }
  else {
    xlarg = M_PI*2.;  x0 = M_PI;
  }
  if (y_larg > 0.) {
    ylarg = y_larg; y0 = y_0;
  }
  else {
    ylarg = M_PI;  y0 = 0.;
  }
  nbparall = np;
  nbmerid = nm;
  nbpt = 101;
  fgrevtheta = false;
}

/* --Methode-- */
MollweideGridDrawer::~MollweideGridDrawer()
{
}

/* --Methode-- */
void MollweideGridDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
{
  PIGrCoord * xxl = new PIGrCoord[nbpt];
  PIGrCoord * xxr = new PIGrCoord[nbpt];
  PIGrCoord * yy = new PIGrCoord[nbpt];

  char buff[64];

  // Trace des paralleles
  g->DrawLine(xmin, y0, xmax, y0);
  g->DrawString(x0+xlarg/100, y0+ylarg/100, "0");

  //  PIGrCoord asc, PIGrCoord desc;
  int i,k;

  int nkp = (nbparall-1)/2;
  if (nkp > 0) {
    double dy = 0.5*ylarg/(double)(nkp+1);
    double dtd = 90./(double)(nkp+1);
    for(k=1; k<=nkp; k++) {
      double fx = sin(acos(2.*(double)k*dy/ylarg))*0.5;
      double yc = k*dy+y0;
      g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc); 
      double ctd = (fgrevtheta) ? -dtd*k : dtd*k;
      sprintf(buff, "%5.1f", ctd);
      g->DrawString(x0, yc, buff);
      g->DrawString(x0-0.3*xlarg, yc, buff);
      g->DrawString(x0+0.25*xlarg, yc, buff);
      yc = -k*dy+y0;
      g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc); 
      sprintf(buff, "%5.1f", -ctd);
      g->DrawString(x0, yc, buff);
      g->DrawString(x0-0.3*xlarg, yc, buff);
      g->DrawString(x0+0.25*xlarg, yc, buff);
    }
  }

  // Trace des meridiennes
  g->DrawLine(x0, ymin, x0, ymax);

  int nkm = (nbmerid-1)/2;
  if (nkm < 1)  nkm = 1;
  double dphi = M_PI/nkm;
  double dpd = 180./(double)(nkm);
  double dy = ylarg/(nbpt-1);
  for(k=0; k<nkm; k++) {
    double phi0 = k*dphi;
    double dphimil = fabs(phi0-M_PI);
    if (dphimil < 1.e-6)  continue;
    for(i=0; i<nbpt; i++) {
      double yc = (double)i*dy+y0-0.5*ylarg;
      double fx = (2.*fabs(yc-y0)/ylarg <= 1.) ? 
        sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5 : 0.;
      yy[i] = yc;
      xxl[i] = x0-fx*xlarg;
      xxr[i] = x0+fx*xlarg;
    }

    g->DrawPolygon(xxl, yy, nbpt, false);
    g->DrawPolygon(xxr, yy, nbpt, false);

    if (k == 0)  continue;
    for(i=-1; i<=1; i++) {
      double yc = 0.20*i*ylarg+y0;
      double fx = sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5;
      double xc = x0-fx*xlarg;
      sprintf(buff, "%4.1f", dpd*k);
      g->DrawString(xc, yc, buff);
      xc = x0+fx*xlarg;
      sprintf(buff, "%4.1f", 360.-dpd*k);
      g->DrawString(xc, yc, buff);
    }
  }

  delete [] xxl;
  delete [] xxr;
  delete [] yy;
}

/* --Methode-- */
void MollweideGridDrawer::UpdateLimits()
{
  SetLimits(x0-xlarg/2., x0+xlarg/2., y0-ylarg/2., y0+ylarg/2.);
}
// Cette fonction calcule les parametres de l'ellipse des erreurs
// a,b 1/2 grand, petit axe ; tet : Angle du grand axe / Ox 
// A partir de s1 = Sigma_1^2 , s2 = Sigma_2^2 , c = Covar_1_2 = Sigma_1 x Sigma_2 x Rho 

//---------------------------------------------------------------------------
//--- Fonctions de trace d'ellipses d'erreur 
//---------------------------------------------------------------------------

/* Nouvelle-Fonction */
int  diag_ellipse(double s1, double s2, double c, 
                  double &a, double &b, double& tet, double& phi, int prt)
{
  int rc = 0;
  // On calcule les valeurs propres de la matrice 
  // | s1  c  |
  // | c   s2 |
  a = b = 0.;  // Demi axes de l'ellipse
  tet = 0.;     // Cos theta de l'angle de l'axe de l'ellipse
  phi = 0.;
  double det = (s1-s2)*(s1-s2) + 4*c*c;
  if (det < 0.) {
    if (prt > 0) 
      cout << "ERREUR : diag_ellipse, ERREUR det = " << det << " <0. s1=" 
           << s1 << " s2=" << s2 << " c=" << c << endl;
    return 99;
  }
  // les deux valeurs propres l1, l2
  double l1 = 0.5*((s1+s2)+sqrt(det));
  double l2 = 0.5*((s1+s2)-sqrt(det));
  if (l1 < 0.) {
    if (prt > 1) 
      cout << "diag_ellipse/Warning l1 < 0. -> a = -sqrt(l1)" << endl; 
    a = -sqrt(-l1);
    rc += 2;
  }
  else   a = sqrt(l1);
  if (l2 < 0.) {
    if (prt > 1) 
      cout << "diag_ellipse/Warning l2 < 0. -> b = -sqrt(l2)" << endl; 
    b = -sqrt(-l2);
    rc += 4;
  }
  else b = sqrt(l2);
  
  if (fabs(c) > 1.e-9)  tet = atan2((l1-s1), c);
  else {
    if (prt > 2) 
      cout << "diag_ellipse/Warning c (Correlation) = 0. -> Tet=0 " << endl; 
  }
  // Calcul de tan(2 phi)
  double t2phi = 0.;
  if (fabs(s1-s2) > 1.e-9) {
    t2phi = 2.*c/(s1-s2);
    phi = atan(t2phi)/2.;
  }
  else {
    if (prt > 3) 
      cout << "diag_ellipse/Warning fabs(s1-s2)=" << fabs(s1-s2)
           << " -> Phi=Pi/4 ... " << endl; 
    phi = 0.25*M_PI;
  }
  if (prt > 0) 
    cout << "diag_ellipse/Info: s1=" << s1 << " s2=" << s2 
         << " c=" << c << " ---> l1=" << l1 << " l2=" << l2 
         << "\n a= " << a << " b=" << b << " Theta= " << tet 
         << "( " << tet*2/M_PI << " xPi/2)  Phi= " << phi 
         << "( " << phi*2/M_PI << " x Pi/2) " << endl;
  return rc;
}

/* Nouvelle-Fonction */
void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy, 
                 double scale, int nbpoints)
{
  a *= scale;
  b *= scale;
  
  int npt = nbpoints;
  double pas = M_PI*2./(npt-1);
  evx.SetSize(npt+1);
  evy.SetSize(npt+1);
  for(int i=0; i<=npt; i++) {
    double t = i*pas;
    double xr = a*cos(t);
    double yr = b*sin(t);
    double ar = atan2(yr, xr);
    double r = sqrt(xr*xr+yr*yr);
    evx(i) = r*cos(ar+tet);
    evy(i) = r*sin(ar+tet);
  }
  return;
}

//---------------------------------------------------------------------------
//--- Classe sopiamoduleExecutor
//---------------------------------------------------------------------------
/* --Methode-- */
sopiamoduleExecutor::sopiamoduleExecutor()
{

PIACmd * mpiac;
NamedObjMgr omg;
mpiac = omg.GetImgApp()->CmdInterpreter();

string kw,usage;
// ------------- Help group FFT 
string hgrp = "FFT";

kw = "fftforw";
usage = "FFT on a vector -> computes forward Fourier transform";
usage += "\n Usage: fftforw vecSig vecSpec [graphic_att] ";
usage += "\n vecSpec = FFTForward(vecSig) ";
usage += "\n vecSig: Input data vector of type r_8 or complex<r_8> ";
usage += "\n vecSpec: Output data vector of type complex<r_8> ";
usage += "\n See also : fftback fftfilter fftfuncfilter ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "fftback";
usage = "FFT on a vector -> computes backward Fourier transform ";
usage += "\n Usage: fftback vecSpec vecSig [graphic_att] [C/Z] ";
usage += "\n vecSig = FFTBackward(vecSpec) ";
usage += "\n vecSpec: Input data vector of type complex<r_8> ";
usage += "\n vecSig:  Output Data vector of type r_8 (default) ";
usage += "\n    or complex<r_8> if C/Z specified or ";
usage += "\n    vecSpec computed by fftforw on a complex vector";
usage += "\n See also : fftforw fftfilter fftfuncfilter ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "fftfilter";
usage = "Filter (multiply) vecSpec (vector complex<r_8>) by Filter (vector<r_8>)";
usage += "\n Usage: fftfilter vecSpec FilterVec vecFiltSpec [graphic_att] ";
usage += "\n   vecFiltSpec(i) = vecSpec(i) * complex(FilterVec(i),0) ";
usage += "\n See also : fftforw fftbackw fftfuncfilter ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "fftfuncfilter";
usage = "Filter (multiply) vecSpec (vector complex<r_8>) by FilterFunc(i) ";
usage += "\n Usage: fftfilter vecSpec FilterFunc vecFiltSpec [graphic_att] ";
usage += "\n   vecFiltSpec(i) = vecSpec(i) * complex(FilterFunc(i),0) ";
usage += "\n See also : fftforw fftbackw fftfilter ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

//--------------------------
hgrp = "SophyaCmd";

kw = "mollgridsph";
usage = "Creates a spherical coordinate grid in Molleweide projection ";
usage += "\n Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att] ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mollgrid";
usage = "Creates a spherical coordinate grid in Molleweide projection ";
usage += "\n Usage: mollgrid [Nb_Parallel Nb_Meridien graphic_att] ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "setprjmoldefval";
usage = "Set default value for Molleweide projection of spherical maps (outside maps)";
usage += "\n Usage: setprjmoldefval OutOfMapValue ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

//--- Trace d'ellipse d'erreur
hgrp = "SophyaCmd";
kw = "errorellipse";
usage = " Error ellipse plotter \n";
usage += "Usage: errellipse [-fill] xc yc sx2 sy2 cov [scale=1] [gr_opt] [npt=128] \n"; 
usage = "   Adds an error ellipse, specified by its center (xc,yc) \n";
usage += "  and error matrix (covariance) parameters  SigX^2,SigY^2,Covar \n" ;
usage += "  A global scaling parameters scale , graphic attributes and \n";
usage += "   nomber of points can optionnaly be specified\n";
usage += "  if the -fill flag specified, plots a filled ellipse \n" ;
mpiac->RegisterCommand(kw, usage, this, hgrp);

}

/* --Methode-- */
sopiamoduleExecutor::~sopiamoduleExecutor()
{
}

/* --Methode-- */
int sopiamoduleExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
{

if (kw == "powerspec" || kw == "fftforw") {
  if (tokens.size() < 2) {
    cout << "Usage: fftforw vecSig vecSpec [graphic_att] " << endl;
    return(0);
    }
  if (tokens.size() < 3)  tokens.push_back((string)"");
  SophyaFFT_forw(tokens[0], tokens[1], tokens[2]);
}
else if (kw == "fftback") {
  if (tokens.size() < 2) {
    cout << "Usage: fftback vecSpec vecSig [graphic_att] [C/Z]" << endl;
    return(0);
    }
  bool forcez8 = false;
  if (tokens.size() > 3) {
    char fcz8 = tolower(tokens[3][0]);
    if ( (fcz8 == 'c') || (fcz8 == 'z') ) forcez8 = true;
  }
  if (tokens.size() < 3)  tokens.push_back((string)"");
  SophyaFFT_back(tokens[0], tokens[1], tokens[2], forcez8);
}
else if ( (kw == "fftfilter") || (kw == "fftfuncfilter") ) {
  if (tokens.size() < 3) {
    cout << "Usage: fftfilter/fftfuncfilter vecSpec Filter vecFiltSpec [graphic_att]" << endl;
    return(0);
    }
  if (tokens.size() < 4)  tokens.push_back((string)"");
  bool fgvec = false;
  if (kw == "fftfilter") fgvec = true; 
  SophyaFFT_filter(tokens[0], tokens[1], tokens[2], tokens[3], fgvec);
}

else if (kw == "prjmoldefval") {
  if (tokens.size() < 1) {
    cout << "Usage: prjmoldefval OutOfMapValue" << endl;
    return(0);
    }
  Set_Project_Mol_DefVal(atof(tokens[0].c_str()) );
}

else if ( (kw == "mollgridsph") || (kw == "mollgrid") ) {
   NamedObjMgr omg;
   MollweideGridDrawer * mollgrid = NULL;
   string dopt="";
   if (kw == "mollgridsph") {
     if (tokens.size() < 1) {
       cout << "Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att]"
            << endl;
       return(0);
     }
     int np = 5;
     int nm = 5;
     dopt = "same,thinline";
     if (tokens.size() > 1)  np = atoi(tokens[1].c_str());
     if (tokens.size() > 2)  nm = atoi(tokens[2].c_str());
     if (tokens.size() > 3)  dopt = tokens[3];
     NObjMgrAdapter* obja = omg.GetObjAdapter(tokens[0]);
     if (obja == NULL) {
       cout << "mollgridsph Error , No object with name " << tokens[0] << endl; 
       return(0);
     }
     P2DArrayAdapter* arr = obja->Get2DArray(dopt);
     if (arr == NULL) {
       cout << "mollgridsph Error , object has non 2DArrayAdapter - Not a spherical map "
            << endl;
       return(0);       
     }
     int nx = arr->XSize();
     int ny = arr->YSize();
     delete arr;
     mollgrid = new MollweideGridDrawer(np, nm, nx, ny);
     cout << " mollgridsph: Creating MollweideGridDrawer() for object" << tokens[0] << endl;
   }
   else if (kw == "mollgrid") {
     int np = 5;
     int nm = 5;
     if (tokens.size() > 0)  np = atoi(tokens[0].c_str());
     if (tokens.size() > 1)  nm = atoi(tokens[1].c_str());
     if (tokens.size() > 2)  dopt = tokens[2];
     mollgrid = new MollweideGridDrawer(np, nm);
     cout << " mollgrid: Creating MollweideGridDrawer(" << np << "," << nm << ")" << endl;
   }

   if (mollgrid == NULL) {
     cout << "mollgridsph Error/BUG !!! mollgrid = NULL "  << endl; 
     return(0);
   }

   int wrsid = 0;
   //bool fgsr = true;
   
   string name = "mollgrid";
   wrsid = omg.GetImgApp()->DispScDrawer(mollgrid, name, dopt); 
 }
else if (kw == "errorellipse") {
  if ((tokens.size()<5) || ((tokens.size()>0)&&(tokens[0]=="-fill")&&(tokens.size()<6))) {
    cout << "Usage: errorellipse [-fill] xc yc sx2 sy2 cov [scale=1 gropt npt=128]" << endl;
    return(1);
  }
  bool fgfill = false;
  unsigned int offtok = 0;
  if (tokens[0] == "-fill") {
    offtok = 1;  fgfill = true;
  }
  double xc = atof(tokens[0+offtok].c_str());
  double yc = atof(tokens[1+offtok].c_str());
  double sx = atof(tokens[2+offtok].c_str());
  double sy = atof(tokens[3+offtok].c_str());
  double cov = atof(tokens[4+offtok].c_str());
  double scale = 1.;
  if (tokens.size() > (5+offtok)) scale = atof(tokens[5+offtok].c_str());
  string gropt;
  if (tokens.size() > (6+offtok)) gropt = tokens[6+offtok];
  int npt = 128;
  if (tokens.size() > (7+offtok)) npt = atoi(tokens[7+offtok].c_str());

  double a,b,tet,phi;
  diag_ellipse(sx, sy, cov, a, b, tet, phi, 1);
  Vector evx, evy;
  ellipse2vec(a,b,tet,evx, evy, scale, npt+1);
  vector<double> vx, vy;
  for(int k=0; k<evx.Size(); k++) {
    vx.push_back(evx(k)+xc);
    vy.push_back(evy(k)+yc);
  }
  NamedObjMgr omg;
  omg.GetImgApp()->AddPoly(vx, vy, gropt, fgfill, false);
}

return(0);

}

static sopiamoduleExecutor * piaerex = NULL;
/* Nouvelle-Fonction */
void sopiamodule_init()
{
// Fonction d'initialisation du module
// Appele par le gestionnaire de modules de piapp (PIACmd::LoadModule())
if (piaerex) delete piaerex;
piaerex = new sopiamoduleExecutor;
}

/* Nouvelle-Fonction */
void sopiamodule_end()
{
// Desactivation du module
if (piaerex) delete piaerex;
piaerex = NULL;
}


/*  Nouvelle-Fonction */
void SophyaFFT_forw(string& nom, string& nomout, string dopt)
{
//Timer tm("powerspec");
NamedObjMgr omg;
AnyDataObj* obj=omg.GetObj(nom);
if (obj == NULL) {
  cout<<"SophyaFFT_forw() Error , Pas d'objet de nom "<<nom<<endl;
  return;
}

TVector<r_8>* vin = dynamic_cast<TVector<r_8> *>(obj);
TVector< complex<r_8> >* vinc = dynamic_cast<TVector< complex<r_8> >*>(obj);
if (vin == NULL && vinc == NULL) {
  cout<<"SophyaFFT_forw() Error , Objet n'est pas un vecteur r_8 ou complex<r_8>"<<endl;
  return;
}
TVector< complex<r_8> > * vout = new TVector< complex<r_8> >  ;
FFTPackServer ffts;
if(vin!=NULL) {
  cout<<"SophyaFFT_forw() - Computing FFT of TVector<r_8> of size "<<vin->NElts()<< endl;
  ffts.FFTForward(*vin, *vout);
  cout<<"...Output TVector<complex<r_8>> of size: "<<vout->NElts()<<endl;
  vout->Info()["FFTTYPE"]="F:R8->Z8";
} else {
  cout<<"SophyaFFT_forw() - Computing FFT of vector TVector<complex<r_8>> of size "<<vinc->NElts()<< endl;
  ffts.FFTForward(*vinc, *vout);
  cout<<"...Output TVector<complex<r_8>> of size: "<<vout->NElts()<<endl;
  vout->Info()["FFTTYPE"] = "F:Z8->Z8";
}
// cout << " SophyaFFT - FFT done " << endl;
 // tm.Split(" FFT done ");

omg.AddObj(vout, nomout);
omg.DisplayObj(nomout, dopt);
return;

}


/*  Nouvelle-Fonction */
void SophyaFFT_back(string& nom, string& nomout, string dopt, bool forcez8)
{
  //Timer tm("powerspec");
NamedObjMgr omg;
AnyDataObj* obj=omg.GetObj(nom);
if (obj == NULL) {
  cout << "SophyaFFT_back() Error , Pas d'objet de nom " << nom << endl;
  return;
}

TVector< complex<r_8> >* vinc = dynamic_cast<TVector< complex<r_8> >*>(obj);
if (vinc == NULL) {
  cout << "SophyaFFT_back() Error , Objet n'est pas un vecteur complex<r_8>" <<  endl;
  return;
  }

TVector< r_8 > * vout = NULL;
TVector< complex<r_8> > * voutc = NULL;

if (forcez8 || (vinc->Info().GetS("FFTTYPE") == "F:Z8->Z8") ) {
  voutc = new TVector< complex<r_8> >;
} 
else {
  vout = new TVector< r_8 >;
}

FFTPackServer ffts;
cout << "SophyaFFT_back() - Computing Backward FFT of vector size " << vinc->NElts() << endl;
if(vout!=NULL) {
  ffts.FFTBackward(*vinc, *vout);
  cout<<"...Output TVector<r_8> of size: "<<vout->NElts()<<endl;
  vout->Info()["FFTTYPE"]="B:Z8->R8";
  omg.AddObj(vout, nomout);
} else {
  ffts.FFTBackward(*vinc, *voutc);
  cout<<"...Output TVector<complex<r_8>> of size: "<<voutc->NElts()<<endl;
  voutc->Info()["FFTTYPE"]="B:Z8->Z8";
  omg.AddObj(voutc, nomout);
}

omg.DisplayObj(nomout, dopt);
return;
}

/*  Nouvelle-Fonction */
void SophyaFFT_crefilter(string& nomvec, string& nomfilter, string func, string dopt)
{
NamedObjMgr omg;
AnyDataObj* obj=omg.GetObj(nomvec);
if(obj == NULL) {
  cout<<"SophyaFFT_crefilter() Error , Pas d'objet de nom "<<nomvec<<endl;
  return;
}

TVector< complex<r_8> >* vfft = dynamic_cast<TVector< complex<r_8> >*>(obj);
if(vfft == NULL) {
  cout<<"SophyaFFT_crefilter() Error , Objet "<<nomvec<<" n'est pas un Tvector<complex<r_8>>"<< endl;
  return;
}
sa_size_t n = vfft->Size();
if(n<=0) return;

cout<<"Creating filter "<<nomfilter<<"("<<n<<") from vec "<<nomvec<<" with fun "<<func<<endl;
omg.GetServiceObj()->PlotFunc(func,nomfilter,0.,(double)n,n,dopt); 
}

/*  Nouvelle-Fonction */
void SophyaFFT_filter(string& nom, string& filt, string& outvec, string dopt, bool fgvec)
{
NamedObjMgr omg;
AnyDataObj* obj=omg.GetObj(nom);
if(obj == NULL) {
  cout<<"SophyaFFT_filter() Error , Pas d'objet de nom "<<nom<<endl;
  return;
}
TVector< complex<r_8> >* vfft = dynamic_cast<TVector< complex<r_8> >*>(obj);
if(vfft == NULL) {
  cout<<"SophyaFFT_filter() Error , Objet "<<nom<<" n'est pas un Tvector<complex<r_8>>"<< endl;
  return;
}
int n = vfft->Size();
if(n<=0) {
  cout<<"SophyaFFT_filter() Error , vector size = 0 " << endl;
  return;
}

string nomfilter = "/autoc/fft_filter";
if (fgvec) {
  nomfilter = filt; 
  cout<<"SophyaFFT_filter(): using vector " << filt << " for filtering spectrum ..." << endl;   
}
else {
  cout<<"SophyaFFT_filter(): Creating filter vector "<<nomfilter
      <<" from function " << filt << endl; 
  string fcrdopt = "nodisp";
  omg.GetServiceObj()->PlotFunc(filt, nomfilter, 0., (double)n, n, fcrdopt); 
}
AnyDataObj* objfilter=omg.GetObj(nomfilter);
if(objfilter == NULL) {
  cout<<"SophyaFFT_filter() Error , Pas d'objet de nom "<<nomfilter<<endl;
  return;
}

TVector<r_8>* vfilter = dynamic_cast<TVector<r_8>*>(objfilter);
if(vfft == NULL) {
  cout<<"SophyaFFT_filter() Error , Objet "<<nomfilter<<" n'est pas un Tvector<r_8>"<<endl;
  return;
}

n = (vfilter->Size()<vfft->Size()) ? vfilter->Size(): vfft->Size();

cout<<" SophyaFFT_filter(): filtering from 0 to "<<n<<endl;
TVector< complex<r_8> > * filtfft = new TVector< complex<r_8> >(*vfft, false);
for(int i=0;i<n;i++) (*filtfft)(i) *= (*vfilter)(i);

cout<<"SophyaFFT_filter(): Adding obj: " << outvec << endl;   
omg.AddObj(filtfft, outvec);
omg.DisplayObj(outvec, dopt);
return;
}