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

Last change on this file since 3430 was 3430, checked in by ansari, 18 years ago

Ajout commande errorellipse (fait pour these de S. Bargot) a sopiamodule.cc - Reza 11/12/2007
File size: 20.6 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <stdio.h>
4#include <stdlib.h>
5#include <iostream>
6#include <ctype.h>
7#include <math.h>
8
9#include <typeinfo>
10
11#include <vector>
12#include <string>
13
14#include "piacmd.h"
15#include "nobjmgr.h"
16#include "pistdimgapp.h"
17#include "servnobjm.h"
18#include "tvector.h"
19#include "pitvmaad.h"
20#include "fftpserver.h"
21#include "bruit.h"
22#include "piscdrawwdg.h"
23#include "ctimer.h"
24
25#include "pidrawer.h"
26#include "nomgadapter.h"
27#include "nomskymapadapter.h"
28
29extern "C" {
30void sopiamodule_init();
31void sopiamodule_end();
32}
33
34// --- fonctions pour faire FFT
35void SophyaFFT_forw(string& nom, string& nomout, string dopt);
36void SophyaFFT_back(string& nom, string& nomout, string dopt, bool forcez8=false);
37// void SophyaFFT_crefilter(string& nomvec, string& nomfilter, string func, string dopt);
38void SophyaFFT_filter(string& nom, string& filt, string& outvec, string dopt, bool fgvec=false);
39
40//--- fonctions pour trace d'ellipse d'erreur
41int diag_ellipse(double s1, double s2, double c,
42 double &a, double &b, double& tet, double& phi, int prt=0);
43
44void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy,
45 double scale=1., int nbpoints=101);
46
47//--- Declaration du CmdExecutor du module
48class sopiamoduleExecutor : public CmdExecutor {
49public:
50 sopiamoduleExecutor();
51 virtual ~sopiamoduleExecutor();
52 virtual int Execute(string& keyw, vector<string>& args, string& toks);
53};
54
55//--- Classe pour trace de grille en projection spherique Mollweide
56class MollweideGridDrawer : public PIDrawer {
57public:
58 MollweideGridDrawer(int np, int nm, int nx, int ny);
59 MollweideGridDrawer(int np, int nm, double x_larg=0., double x_0=0.,
60 double y_larg=0., double y_0=0.);
61 virtual ~MollweideGridDrawer();
62 virtual void Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax);
63 virtual void UpdateLimits();
64
65protected:
66 double xlarg, x0;
67 double ylarg, y0;
68 int nbparall, nbmerid;
69 int nbpt;
70 bool fgrevtheta;
71};
72
73/* --Methode-- */
74MollweideGridDrawer::MollweideGridDrawer(int np, int nm, int nx, int ny)
75{
76 x0 = (double)nx/2.;
77 y0 = (double)ny/2.;
78 xlarg = nx;
79 ylarg = ny;
80 nbparall = np;
81 nbmerid = nm;
82 nbpt = ny/5;
83 if ((nbpt % 20) == 0) nbpt++;
84 fgrevtheta = true;
85}
86
87/* --Methode-- */
88MollweideGridDrawer::MollweideGridDrawer(int np, int nm, double x_larg, double x_0,
89 double y_larg, double y_0)
90{
91 if (x_larg > 0.) {
92 xlarg = x_larg; x0 = x_0;
93 }
94 else {
95 xlarg = M_PI*2.; x0 = M_PI;
96 }
97 if (y_larg > 0.) {
98 ylarg = y_larg; y0 = y_0;
99 }
100 else {
101 ylarg = M_PI; y0 = 0.;
102 }
103 nbparall = np;
104 nbmerid = nm;
105 nbpt = 101;
106 fgrevtheta = false;
107}
108
109/* --Methode-- */
110MollweideGridDrawer::~MollweideGridDrawer()
111{
112}
113
114/* --Methode-- */
115void MollweideGridDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
116{
117 PIGrCoord * xxl = new PIGrCoord[nbpt];
118 PIGrCoord * xxr = new PIGrCoord[nbpt];
119 PIGrCoord * yy = new PIGrCoord[nbpt];
120
121 char buff[64];
122
123 // Trace des paralleles
124 g->DrawLine(xmin, y0, xmax, y0);
125 g->DrawString(x0+xlarg/100, y0+ylarg/100, "0");
126
127 // PIGrCoord asc, PIGrCoord desc;
128 int i,k;
129
130 int nkp = (nbparall-1)/2;
131 if (nkp > 0) {
132 double dy = 0.5*ylarg/(double)(nkp+1);
133 double dtd = 90./(double)(nkp+1);
134 for(k=1; k<=nkp; k++) {
135 double fx = sin(acos(2.*(double)k*dy/ylarg))*0.5;
136 double yc = k*dy+y0;
137 g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc);
138 double ctd = (fgrevtheta) ? -dtd*k : dtd*k;
139 sprintf(buff, "%5.1f", ctd);
140 g->DrawString(x0, yc, buff);
141 g->DrawString(x0-0.3*xlarg, yc, buff);
142 g->DrawString(x0+0.25*xlarg, yc, buff);
143 yc = -k*dy+y0;
144 g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc);
145 sprintf(buff, "%5.1f", -ctd);
146 g->DrawString(x0, yc, buff);
147 g->DrawString(x0-0.3*xlarg, yc, buff);
148 g->DrawString(x0+0.25*xlarg, yc, buff);
149 }
150 }
151
152 // Trace des meridiennes
153 g->DrawLine(x0, ymin, x0, ymax);
154
155 int nkm = (nbmerid-1)/2;
156 if (nkm < 1) nkm = 1;
157 double dphi = M_PI/nkm;
158 double dpd = 180./(double)(nkm);
159 double dy = ylarg/(nbpt-1);
160 for(k=0; k<nkm; k++) {
161 double phi0 = k*dphi;
162 double dphimil = fabs(phi0-M_PI);
163 if (dphimil < 1.e-6) continue;
164 for(i=0; i<nbpt; i++) {
165 double yc = (double)i*dy+y0-0.5*ylarg;
166 double fx = (2.*fabs(yc-y0)/ylarg <= 1.) ?
167 sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5 : 0.;
168 yy[i] = yc;
169 xxl[i] = x0-fx*xlarg;
170 xxr[i] = x0+fx*xlarg;
171 }
172
173 g->DrawPolygon(xxl, yy, nbpt, false);
174 g->DrawPolygon(xxr, yy, nbpt, false);
175
176 if (k == 0) continue;
177 for(i=-1; i<=1; i++) {
178 double yc = 0.20*i*ylarg+y0;
179 double fx = sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5;
180 double xc = x0-fx*xlarg;
181 sprintf(buff, "%4.1f", dpd*k);
182 g->DrawString(xc, yc, buff);
183 xc = x0+fx*xlarg;
184 sprintf(buff, "%4.1f", 360.-dpd*k);
185 g->DrawString(xc, yc, buff);
186 }
187 }
188
189 delete [] xxl;
190 delete [] xxr;
191 delete [] yy;
192}
193
194/* --Methode-- */
195void MollweideGridDrawer::UpdateLimits()
196{
197 SetLimits(x0-xlarg/2., x0+xlarg/2., y0-ylarg/2., y0+ylarg/2.);
198}
199// Cette fonction calcule les parametres de l'ellipse des erreurs
200// a,b 1/2 grand, petit axe ; tet : Angle du grand axe / Ox
201// A partir de s1 = Sigma_1^2 , s2 = Sigma_2^2 , c = Covar_1_2 = Sigma_1 x Sigma_2 x Rho
202
203//---------------------------------------------------------------------------
204//--- Fonctions de trace d'ellipses d'erreur
205//---------------------------------------------------------------------------
206
207/* Nouvelle-Fonction */
208int diag_ellipse(double s1, double s2, double c,
209 double &a, double &b, double& tet, double& phi, int prt)
210{
211 int rc = 0;
212 // On calcule les valeurs propres de la matrice
213 // | s1 c |
214 // | c s2 |
215 a = b = 0.; // Demi axes de l'ellipse
216 tet = 0.; // Cos theta de l'angle de l'axe de l'ellipse
217 phi = 0.;
218 double det = (s1-s2)*(s1-s2) + 4*c*c;
219 if (det < 0.) {
220 if (prt > 0)
221 cout << "ERREUR : diag_ellipse, ERREUR det = " << det << " <0. s1="
222 << s1 << " s2=" << s2 << " c=" << c << endl;
223 return 99;
224 }
225 // les deux valeurs propres l1, l2
226 double l1 = 0.5*((s1+s2)+sqrt(det));
227 double l2 = 0.5*((s1+s2)-sqrt(det));
228 if (l1 < 0.) {
229 if (prt > 1)
230 cout << "diag_ellipse/Warning l1 < 0. -> a = -sqrt(l1)" << endl;
231 a = -sqrt(-l1);
232 rc += 2;
233 }
234 else a = sqrt(l1);
235 if (l2 < 0.) {
236 if (prt > 1)
237 cout << "diag_ellipse/Warning l2 < 0. -> b = -sqrt(l2)" << endl;
238 b = -sqrt(-l2);
239 rc += 4;
240 }
241 else b = sqrt(l2);
242
243 if (fabs(c) > 1.e-9) tet = atan2((l1-s1), c);
244 else {
245 if (prt > 2)
246 cout << "diag_ellipse/Warning c (Correlation) = 0. -> Tet=0 " << endl;
247 }
248 // Calcul de tan(2 phi)
249 double t2phi = 0.;
250 if (fabs(s1-s2) > 1.e-9) {
251 t2phi = 2.*c/(s1-s2);
252 phi = atan(t2phi)/2.;
253 }
254 else {
255 if (prt > 3)
256 cout << "diag_ellipse/Warning fabs(s1-s2)=" << fabs(s1-s2)
257 << " -> Phi=Pi/4 ... " << endl;
258 phi = 0.25*M_PI;
259 }
260 if (prt > 0)
261 cout << "diag_ellipse/Info: s1=" << s1 << " s2=" << s2
262 << " c=" << c << " ---> l1=" << l1 << " l2=" << l2
263 << "\n a= " << a << " b=" << b << " Theta= " << tet
264 << "( " << tet*2/M_PI << " xPi/2) Phi= " << phi
265 << "( " << phi*2/M_PI << " x Pi/2) " << endl;
266 return rc;
267}
268
269/* Nouvelle-Fonction */
270void ellipse2vec(double a, double b, double tet, Vector& evx, Vector& evy,
271 double scale, int nbpoints)
272{
273 a *= scale;
274 b *= scale;
275
276 int npt = nbpoints;
277 double pas = M_PI*2./(npt-1);
278 evx.SetSize(npt+1);
279 evy.SetSize(npt+1);
280 for(int i=0; i<=npt; i++) {
281 double t = i*pas;
282 double xr = a*cos(t);
283 double yr = b*sin(t);
284 double ar = atan2(yr, xr);
285 double r = sqrt(xr*xr+yr*yr);
286 evx(i) = r*cos(ar+tet);
287 evy(i) = r*sin(ar+tet);
288 }
289 return;
290}
291
292//---------------------------------------------------------------------------
293//--- Classe sopiamoduleExecutor
294//---------------------------------------------------------------------------
295/* --Methode-- */
296sopiamoduleExecutor::sopiamoduleExecutor()
297{
298
299PIACmd * mpiac;
300NamedObjMgr omg;
301mpiac = omg.GetImgApp()->CmdInterpreter();
302
303string kw,usage;
304// ------------- Help group FFT
305string hgrp = "FFT";
306
307kw = "fftforw";
308usage = "FFT on a vector -> computes forward Fourier transform";
309usage += "\n Usage: fftforw vecSig vecSpec [graphic_att] ";
310usage += "\n vecSpec = FFTForward(vecSig) ";
311usage += "\n vecSig: Input data vector of type r_8 or complex<r_8> ";
312usage += "\n vecSpec: Output data vector of type complex<r_8> ";
313usage += "\n See also : fftback fftfilter fftfuncfilter ";
314mpiac->RegisterCommand(kw, usage, this, hgrp);
315
316kw = "fftback";
317usage = "FFT on a vector -> computes backward Fourier transform ";
318usage += "\n Usage: fftback vecSpec vecSig [graphic_att] [C/Z] ";
319usage += "\n vecSig = FFTBackward(vecSpec) ";
320usage += "\n vecSpec: Input data vector of type complex<r_8> ";
321usage += "\n vecSig: Output Data vector of type r_8 (default) ";
322usage += "\n or complex<r_8> if C/Z specified or ";
323usage += "\n vecSpec computed by fftforw on a complex vector";
324usage += "\n See also : fftforw fftfilter fftfuncfilter ";
325mpiac->RegisterCommand(kw, usage, this, hgrp);
326
327kw = "fftfilter";
328usage = "Filter (multiply) vecSpec (vector complex<r_8>) by Filter (vector<r_8>)";
329usage += "\n Usage: fftfilter vecSpec FilterVec vecFiltSpec [graphic_att] ";
330usage += "\n vecFiltSpec(i) = vecSpec(i) * complex(FilterVec(i),0) ";
331usage += "\n See also : fftforw fftbackw fftfuncfilter ";
332mpiac->RegisterCommand(kw, usage, this, hgrp);
333
334kw = "fftfuncfilter";
335usage = "Filter (multiply) vecSpec (vector complex<r_8>) by FilterFunc(i) ";
336usage += "\n Usage: fftfilter vecSpec FilterFunc vecFiltSpec [graphic_att] ";
337usage += "\n vecFiltSpec(i) = vecSpec(i) * complex(FilterFunc(i),0) ";
338usage += "\n See also : fftforw fftbackw fftfilter ";
339mpiac->RegisterCommand(kw, usage, this, hgrp);
340
341//--------------------------
342hgrp = "SophyaCmd";
343
344kw = "mollgridsph";
345usage = "Creates a spherical coordinate grid in Molleweide projection ";
346usage += "\n Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att] ";
347mpiac->RegisterCommand(kw, usage, this, hgrp);
348
349kw = "mollgrid";
350usage = "Creates a spherical coordinate grid in Molleweide projection ";
351usage += "\n Usage: mollgrid [Nb_Parallel Nb_Meridien graphic_att] ";
352mpiac->RegisterCommand(kw, usage, this, hgrp);
353
354kw = "setprjmoldefval";
355usage = "Set default value for Molleweide projection of spherical maps (outside maps)";
356usage += "\n Usage: setprjmoldefval OutOfMapValue ";
357mpiac->RegisterCommand(kw, usage, this, hgrp);
358
359//--- Trace d'ellipse d'erreur
360hgrp = "SophyaCmd";
361kw = "errorellipse";
362usage = " Error ellipse plotter \n";
363usage += "Usage: errellipse [-fill] xc yc sx2 sy2 cov [scale=1] [gr_opt] [npt=128] \n";
364usage = " Adds an error ellipse, specified by its center (xc,yc) \n";
365usage += " and error matrix (covariance) parameters SigX^2,SigY^2,Covar \n" ;
366usage += " A global scaling parameters scale , graphic attributes and \n";
367usage += " nomber of points can optionnaly be specified\n";
368usage += " if the -fill flag specified, plots a filled ellipse \n" ;
369mpiac->RegisterCommand(kw, usage, this, hgrp);
370
371}
372
373/* --Methode-- */
374sopiamoduleExecutor::~sopiamoduleExecutor()
375{
376}
377
378/* --Methode-- */
379int sopiamoduleExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
380{
381
382if (kw == "powerspec" || kw == "fftforw") {
383 if (tokens.size() < 2) {
384 cout << "Usage: fftforw vecSig vecSpec [graphic_att] " << endl;
385 return(0);
386 }
387 if (tokens.size() < 3) tokens.push_back((string)"");
388 SophyaFFT_forw(tokens[0], tokens[1], tokens[2]);
389}
390else if (kw == "fftback") {
391 if (tokens.size() < 2) {
392 cout << "Usage: fftback vecSpec vecSig [graphic_att] [C/Z]" << endl;
393 return(0);
394 }
395 bool forcez8 = false;
396 if (tokens.size() > 3) {
397 char fcz8 = tolower(tokens[3][0]);
398 if ( (fcz8 == 'c') || (fcz8 == 'z') ) forcez8 = true;
399 }
400 if (tokens.size() < 3) tokens.push_back((string)"");
401 SophyaFFT_back(tokens[0], tokens[1], tokens[2], forcez8);
402}
403else if ( (kw == "fftfilter") || (kw == "fftfuncfilter") ) {
404 if (tokens.size() < 3) {
405 cout << "Usage: fftfilter/fftfuncfilter vecSpec Filter vecFiltSpec [graphic_att]" << endl;
406 return(0);
407 }
408 if (tokens.size() < 4) tokens.push_back((string)"");
409 bool fgvec = false;
410 if (kw == "fftfilter") fgvec = true;
411 SophyaFFT_filter(tokens[0], tokens[1], tokens[2], tokens[3], fgvec);
412}
413
414else if (kw == "prjmoldefval") {
415 if (tokens.size() < 1) {
416 cout << "Usage: prjmoldefval OutOfMapValue" << endl;
417 return(0);
418 }
419 Set_Project_Mol_DefVal(atof(tokens[0].c_str()) );
420}
421
422else if ( (kw == "mollgridsph") || (kw == "mollgrid") ) {
423 NamedObjMgr omg;
424 MollweideGridDrawer * mollgrid = NULL;
425 string dopt="";
426 if (kw == "mollgridsph") {
427 if (tokens.size() < 1) {
428 cout << "Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att]"
429 << endl;
430 return(0);
431 }
432 int np = 5;
433 int nm = 5;
434 dopt = "same,thinline";
435 if (tokens.size() > 1) np = atoi(tokens[1].c_str());
436 if (tokens.size() > 2) nm = atoi(tokens[2].c_str());
437 if (tokens.size() > 3) dopt = tokens[3];
438 NObjMgrAdapter* obja = omg.GetObjAdapter(tokens[0]);
439 if (obja == NULL) {
440 cout << "mollgridsph Error , No object with name " << tokens[0] << endl;
441 return(0);
442 }
443 P2DArrayAdapter* arr = obja->Get2DArray(dopt);
444 if (arr == NULL) {
445 cout << "mollgridsph Error , object has non 2DArrayAdapter - Not a spherical map "
446 << endl;
447 return(0);
448 }
449 int nx = arr->XSize();
450 int ny = arr->YSize();
451 delete arr;
452 mollgrid = new MollweideGridDrawer(np, nm, nx, ny);
453 cout << " mollgridsph: Creating MollweideGridDrawer() for object" << tokens[0] << endl;
454 }
455 else if (kw == "mollgrid") {
456 int np = 5;
457 int nm = 5;
458 if (tokens.size() > 0) np = atoi(tokens[0].c_str());
459 if (tokens.size() > 1) nm = atoi(tokens[1].c_str());
460 if (tokens.size() > 2) dopt = tokens[2];
461 mollgrid = new MollweideGridDrawer(np, nm);
462 cout << " mollgrid: Creating MollweideGridDrawer(" << np << "," << nm << ")" << endl;
463 }
464
465 if (mollgrid == NULL) {
466 cout << "mollgridsph Error/BUG !!! mollgrid = NULL " << endl;
467 return(0);
468 }
469
470 int wrsid = 0;
471 //bool fgsr = true;
472
473 string name = "mollgrid";
474 wrsid = omg.GetImgApp()->DispScDrawer(mollgrid, name, dopt);
475 }
476else if (kw == "errorellipse") {
477 if ((tokens.size()<5) || ((tokens.size()>0)&&(tokens[0]=="-fill")&&(tokens.size()<6))) {
478 cout << "Usage: errorellipse [-fill] xc yc sx2 sy2 cov [scale=1 gropt npt=128]" << endl;
479 return(1);
480 }
481 bool fgfill = false;
482 unsigned int offtok = 0;
483 if (tokens[0] == "-fill") {
484 offtok = 1; fgfill = true;
485 }
486 double xc = atof(tokens[0+offtok].c_str());
487 double yc = atof(tokens[1+offtok].c_str());
488 double sx = atof(tokens[2+offtok].c_str());
489 double sy = atof(tokens[3+offtok].c_str());
490 double cov = atof(tokens[4+offtok].c_str());
491 double scale = 1.;
492 if (tokens.size() > (5+offtok)) scale = atof(tokens[5+offtok].c_str());
493 string gropt;
494 if (tokens.size() > (6+offtok)) gropt = tokens[6+offtok];
495 int npt = 128;
496 if (tokens.size() > (7+offtok)) npt = atoi(tokens[7+offtok].c_str());
497
498 double a,b,tet,phi;
499 diag_ellipse(sx, sy, cov, a, b, tet, phi, 1);
500 Vector evx, evy;
501 ellipse2vec(a,b,tet,evx, evy, scale, npt+1);
502 vector<double> vx, vy;
503 for(int k=0; k<evx.Size(); k++) {
504 vx.push_back(evx(k)+xc);
505 vy.push_back(evy(k)+yc);
506 }
507 NamedObjMgr omg;
508 omg.GetImgApp()->AddPoly(vx, vy, gropt, fgfill, false);
509}
510
511return(0);
512
513}
514
515static sopiamoduleExecutor * piaerex = NULL;
516/* Nouvelle-Fonction */
517void sopiamodule_init()
518{
519// Fonction d'initialisation du module
520// Appele par le gestionnaire de modules de piapp (PIACmd::LoadModule())
521if (piaerex) delete piaerex;
522piaerex = new sopiamoduleExecutor;
523}
524
525/* Nouvelle-Fonction */
526void sopiamodule_end()
527{
528// Desactivation du module
529if (piaerex) delete piaerex;
530piaerex = NULL;
531}
532
533
534/* Nouvelle-Fonction */
535void SophyaFFT_forw(string& nom, string& nomout, string dopt)
536{
537//Timer tm("powerspec");
538NamedObjMgr omg;
539AnyDataObj* obj=omg.GetObj(nom);
540if (obj == NULL) {
541 cout<<"SophyaFFT_forw() Error , Pas d'objet de nom "<<nom<<endl;
542 return;
543}
544
545TVector<r_8>* vin = dynamic_cast<TVector<r_8> *>(obj);
546TVector< complex<r_8> >* vinc = dynamic_cast<TVector< complex<r_8> >*>(obj);
547if (vin == NULL && vinc == NULL) {
548 cout<<"SophyaFFT_forw() Error , Objet n'est pas un vecteur r_8 ou complex<r_8>"<<endl;
549 return;
550}
551TVector< complex<r_8> > * vout = new TVector< complex<r_8> > ;
552FFTPackServer ffts;
553if(vin!=NULL) {
554 cout<<"SophyaFFT_forw() - Computing FFT of TVector<r_8> of size "<<vin->NElts()<< endl;
555 ffts.FFTForward(*vin, *vout);
556 cout<<"...Output TVector<complex<r_8>> of size: "<<vout->NElts()<<endl;
557 vout->Info()["FFTTYPE"]="F:R8->Z8";
558} else {
559 cout<<"SophyaFFT_forw() - Computing FFT of vector TVector<complex<r_8>> of size "<<vinc->NElts()<< endl;
560 ffts.FFTForward(*vinc, *vout);
561 cout<<"...Output TVector<complex<r_8>> of size: "<<vout->NElts()<<endl;
562 vout->Info()["FFTTYPE"] = "F:Z8->Z8";
563}
564// cout << " SophyaFFT - FFT done " << endl;
565 // tm.Split(" FFT done ");
566
567omg.AddObj(vout, nomout);
568omg.DisplayObj(nomout, dopt);
569return;
570
571}
572
573
574/* Nouvelle-Fonction */
575void SophyaFFT_back(string& nom, string& nomout, string dopt, bool forcez8)
576{
577 //Timer tm("powerspec");
578NamedObjMgr omg;
579AnyDataObj* obj=omg.GetObj(nom);
580if (obj == NULL) {
581 cout << "SophyaFFT_back() Error , Pas d'objet de nom " << nom << endl;
582 return;
583}
584
585TVector< complex<r_8> >* vinc = dynamic_cast<TVector< complex<r_8> >*>(obj);
586if (vinc == NULL) {
587 cout << "SophyaFFT_back() Error , Objet n'est pas un vecteur complex<r_8>" << endl;
588 return;
589 }
590
591TVector< r_8 > * vout = NULL;
592TVector< complex<r_8> > * voutc = NULL;
593
594if (forcez8 || (vinc->Info().GetS("FFTTYPE") == "F:Z8->Z8") ) {
595 voutc = new TVector< complex<r_8> >;
596}
597else {
598 vout = new TVector< r_8 >;
599}
600
601FFTPackServer ffts;
602cout << "SophyaFFT_back() - Computing Backward FFT of vector size " << vinc->NElts() << endl;
603if(vout!=NULL) {
604 ffts.FFTBackward(*vinc, *vout);
605 cout<<"...Output TVector<r_8> of size: "<<vout->NElts()<<endl;
606 vout->Info()["FFTTYPE"]="B:Z8->R8";
607 omg.AddObj(vout, nomout);
608} else {
609 ffts.FFTBackward(*vinc, *voutc);
610 cout<<"...Output TVector<complex<r_8>> of size: "<<voutc->NElts()<<endl;
611 voutc->Info()["FFTTYPE"]="B:Z8->Z8";
612 omg.AddObj(voutc, nomout);
613}
614
615omg.DisplayObj(nomout, dopt);
616return;
617}
618
619/* Nouvelle-Fonction */
620void SophyaFFT_crefilter(string& nomvec, string& nomfilter, string func, string dopt)
621{
622NamedObjMgr omg;
623AnyDataObj* obj=omg.GetObj(nomvec);
624if(obj == NULL) {
625 cout<<"SophyaFFT_crefilter() Error , Pas d'objet de nom "<<nomvec<<endl;
626 return;
627}
628
629TVector< complex<r_8> >* vfft = dynamic_cast<TVector< complex<r_8> >*>(obj);
630if(vfft == NULL) {
631 cout<<"SophyaFFT_crefilter() Error , Objet "<<nomvec<<" n'est pas un Tvector<complex<r_8>>"<< endl;
632 return;
633}
634sa_size_t n = vfft->Size();
635if(n<=0) return;
636
637cout<<"Creating filter "<<nomfilter<<"("<<n<<") from vec "<<nomvec<<" with fun "<<func<<endl;
638omg.GetServiceObj()->PlotFunc(func,nomfilter,0.,(double)n,n,dopt);
639}
640
641/* Nouvelle-Fonction */
642void SophyaFFT_filter(string& nom, string& filt, string& outvec, string dopt, bool fgvec)
643{
644NamedObjMgr omg;
645AnyDataObj* obj=omg.GetObj(nom);
646if(obj == NULL) {
647 cout<<"SophyaFFT_filter() Error , Pas d'objet de nom "<<nom<<endl;
648 return;
649}
650TVector< complex<r_8> >* vfft = dynamic_cast<TVector< complex<r_8> >*>(obj);
651if(vfft == NULL) {
652 cout<<"SophyaFFT_filter() Error , Objet "<<nom<<" n'est pas un Tvector<complex<r_8>>"<< endl;
653 return;
654}
655int n = vfft->Size();
656if(n<=0) {
657 cout<<"SophyaFFT_filter() Error , vector size = 0 " << endl;
658 return;
659}
660
661string nomfilter = "/autoc/fft_filter";
662if (fgvec) {
663 nomfilter = filt;
664 cout<<"SophyaFFT_filter(): using vector " << filt << " for filtering spectrum ..." << endl;
665}
666else {
667 cout<<"SophyaFFT_filter(): Creating filter vector "<<nomfilter
668 <<" from function " << filt << endl;
669 string fcrdopt = "nodisp";
670 omg.GetServiceObj()->PlotFunc(filt, nomfilter, 0., (double)n, n, fcrdopt);
671}
672AnyDataObj* objfilter=omg.GetObj(nomfilter);
673if(objfilter == NULL) {
674 cout<<"SophyaFFT_filter() Error , Pas d'objet de nom "<<nomfilter<<endl;
675 return;
676}
677
678TVector<r_8>* vfilter = dynamic_cast<TVector<r_8>*>(objfilter);
679if(vfft == NULL) {
680 cout<<"SophyaFFT_filter() Error , Objet "<<nomfilter<<" n'est pas un Tvector<r_8>"<<endl;
681 return;
682}
683
684n = (vfilter->Size()<vfft->Size()) ? vfilter->Size(): vfft->Size();
685
686cout<<" SophyaFFT_filter(): filtering from 0 to "<<n<<endl;
687TVector< complex<r_8> > * filtfft = new TVector< complex<r_8> >(*vfft, false);
688for(int i=0;i<n;i++) (*filtfft)(i) *= (*vfilter)(i);
689
690cout<<"SophyaFFT_filter(): Adding obj: " << outvec << endl;
691omg.AddObj(filtfft, outvec);
692omg.DisplayObj(outvec, dopt);
693return;
694}
Note: See TracBrowser for help on using the repository browser.