NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
fourier.hpp
Go to the documentation of this file.
1 #ifndef _FOURIER_H
2 #define _FOURIER_H
3 
4 
5 #include <cmath>
6 #include <cassert>
7 #include <iostream>
8 #include <complex>
9 #include <vector>
10 #include <stdlib.h>
11 
12 
13 
14 namespace fourier {
15 
16 void get_dftcorr_linear(double th,double *corfac,std::complex<double> *endcor);
17 void get_dftcorr_cubic(double th,double *corfac,std::complex<double> *endcor);
18 
19 void complex_dftcor_cubic(double w, double delta,double a,double b,
20  std::complex<double> *endpts,std::complex<double> *endcor,double *corfac);
21 void dft_cplx(double w,int n,double a,double b,std::complex<double> *f,
22  std::complex<double> &res, std::complex<double> &err);
23 
24 
25 #define PI 3.14159265358979323846
26 #define ADFT_MINPTS 16
27 
42 class adft_func{
43  public:
44  typedef std::complex<double> cplx;
46  data_= 0;
47  size_=0;
48  tot_=0;
49  nmax_=0;
50  n_=0;
51  a_=0;
52  b_=0;
53  len_=0;
54  f_=0;
55  }
57 
74  adft_func(int size,int nmax){
75  assert(size>0 && nmax>0);
76  data_=new cplx [size];
77  size_=size;
78  tot_=0;
79  nmax_=nmax;
80  n_=0;
81  a_=new double [nmax];
82  b_=new double [nmax];
83  len_=new int [nmax];
84  f_=new cplx* [nmax];
85  }
87  ~adft_func(){
88  delete [] a_;
89  delete [] b_;
90  delete [] len_;
91  delete [] data_;
92  delete [] f_;
93  }
95  adft_func &operator=(const adft_func &a){
96  std::cerr << "= not implemented for adft" << std::endl;
97  abort();
98  }
100  adft_func(const adft_func &a){
101  std::cerr << " copy constructor not implemented for for adft" << std::endl;
102  abort();
103  }
105 
122  void resize(int size,int nmax){ // a resize always cancels all entries !!!
123  if(size<0 || nmax<0) {std::cerr << "adft: size<0 or nmax<0" << std::endl;abort();}
124  delete [] data_;
125  delete [] a_;
126  delete [] b_;
127  delete [] len_;
128  delete [] f_;
129  data_=new cplx [size];
130  size_=size;
131  tot_=0;
132  nmax_=nmax;
133  n_=0;
134  a_=new double [nmax];
135  b_=new double [nmax];
136  len_=new int [nmax];
137  f_=new cplx* [nmax];
138  }
140 
155  cplx *interval_data(int i){
156  if(i>=n_) {std::cerr << "interval i="<<i<<" >n="<< n_<<" in adft" << std::endl;abort();}
157  return f_[i];
158  }
160 
181  template < class function>
182  void init(double a,double b,int n,function &fz){
183  int i;
184  cplx *f;
185  double x,dx;
186 
187  if(a>=b){std::cerr << "adft_func: a>=b" << std::endl;abort();}
188  if(n<1){std::cerr << "adft_func: n<1" << std::endl;abort();}
189  if(nmax_<1 || size_<n) resize(n+1,1);
190  n_=1;
191  tot_=n+1;
192  len_[0]=n;
193  a_[0]=a;
194  b_[0]=b;
195  f=data_;
196  dx=(b-a)/((double) n);
197  for(i=0;i<=n;i++){
198  x=a+i*dx;
199  f[i]=fz(x);
200  }
201  f_[0]=data_;
202  }
204 
221  template < class function>
222  void split(int interv,function &fz){
223  // split interval interv by doubling points
224  int j,j0,j1,jm,len0,n0=n_,i1=interv+1,tot0;
225  double a0,b0,mid,x,dx;
226  cplx *f0,*f1;
227 
228 
229  if(interv>=n0 || interv<0) {std::cerr << "adft_cplx_split: interv out of range" << std::endl; abort();}
230  n0=n_;
231  tot0=tot_;
232  a0=a_[interv];
233  b0=b_[interv];
234  len0=len_[interv];
235 
236  n_++;
237  tot_ += len_[interv] + 1 ;
238  if(n_ > nmax_ || tot_ > size_ ) {std::cerr << "adft_cplx_split: something too small" << std::endl; abort();}
239 
240  /*shift all intervals after interv*/
241  for(j=n0-1;j>interv;j--){
242  a_[j+1]=a_[j];
243  b_[j+1]=b_[j];
244  len_[j+1]=len_[j];
245  f_[j+1]=f_[j];
246  }
247  /*new interval-edges*/
248  mid=(a0+b0)/2.0;
249  b_[interv]=mid;
250  b_[i1]=b0;
251  a_[i1]=mid;
252  a_[interv]=a0;
253  len_[i1]=len0;
254  /*reshuffle data and recalculate*/
255  f0=interval_data(interv);
256  f1=data_+tot0;
257  f_[i1]=f1;
258  dx=(b0-a0)/(2.0*len0);
259  if(len0%2==0){
260  f1[len0]=f0[len0];
261  jm=len0/2;
262  j1=len0;
263  for(j0=len0-1;j0>=jm;j0--){
264  j1--;
265  x=mid+dx*j1;
266  f1[j1]=fz(x);
267  j1--;
268  f1[j1]=f0[j0];
269  }
270  j1=len0;
271  for(j0=jm;j0>0;j0--){
272  f0[j1]=f0[j0];
273  j1--;
274  x=a0+dx*j1;
275  f0[j1]=fz(x);
276  j1--;
277  }
278  }else{
279  jm=(len0+1)/2;
280  j1=len0;
281  for(j0=len0;j0>=jm;j0--){
282  f1[j1]=f0[j0];
283  j1--;
284  x=mid+dx*j1;
285  f1[j1]=fz(x);
286  j1--;
287  }
288  j1=len0;
289  for(j0=jm-1;j0>=0;j0--){
290  x=a0+dx*j1;
291  f0[j1]=fz(x);
292  j1--;
293  f0[j1]=f0[j0];
294  j1--;
295  }
296  }
297  }
317  void dft(double w,cplx &result,cplx &err){
318  int i,n_i;
319  double a,b;
320  cplx res_i,err_i,*f_i;
321  if(n_<=0) {std::cerr << "adft_cplx: no interval in f" << std::endl; abort();}
322  result=0;
323  err=0;
324 
325  for(i=0;i<n_;i++){
326  f_i=interval_data(i);
327  a=a_[i];
328  b=b_[i];
329  n_i=len_[i];
330  fourier::dft_cplx(w,n_i,a,b,f_i,res_i,err_i);
331  result += res_i;
332  err += err_i;
333  }
334  }
335  // template <class T>
336  // void sample(double w,double a,double b,T (*fz)(double),int n,int limit){
337 
365  template <class function>
366  void sample(double w,double a,double b,function &fz,int n,int limit){
367  // n: number of points in each interval
368  // limit: max number of intervals
369  int i,n_cur,isplit;
370  double abserr,abserr_i;
371  cplx z,res_i,err_i,res_l,err_l,res_r,err_r;
372  std::vector<cplx> results(limit),errors(limit);
373 
374  if(n<ADFT_MINPTS || n%2!=0){ std::cerr << "adft_cplx_sample: n<16 || n%2!=0" << std::endl; abort();}
375  if(limit<1){std::cerr << "adft_cplx_sample: limit<1" << std::endl; abort();}
376  if(a>=b){ std::cerr << "adft_cplx_sample: a>=b" << std::endl;abort();}
377 
378  if(limit != nmax_ || (n+1)*limit != size_ ) resize((n+1)*limit,limit);
379  // Sample f in first interval
380  init(a,b,n,fz);
381  if(limit==1) return;
382  // split first interval
383  split(0,fz);
384  // do the first fourier transform
385  fourier::dft_cplx(w,n,a_[0],b_[0],interval_data(0),results[0],errors[0]);
386  fourier::dft_cplx(w,n,a_[1],b_[1],interval_data(1),results[1],errors[1]);
387  while(n_+1<limit){
388  // look for largest error in res-table
389  n_cur=n_;
390  abserr=-1.0;
391  isplit=-1;
392  for(i=0;i<n_cur;i++){
393  z=errors[i];
394  abserr_i=sqrt(z.real()*z.real()+z.imag()*z.imag());
395  if(abserr_i>abserr){
396  isplit=i;
397  abserr=abserr_i;
398  }
399  }
400  // split interval isplit and integrate partial intervals
401 
402  //std::cout << "isplit= "<< isplit << "n= "<< n_ << std::endl;
403 // std::cout << "n " << n_ << "i " << isplit << "[ " << a_[isplit] << "," << b_[isplit] << "]" << std::endl;
404  split(isplit,fz);
405 
406  fourier::dft_cplx(w,n,a_[isplit],b_[isplit],interval_data(isplit),res_l,err_l);
407  fourier::dft_cplx(w,n,a_[isplit+1],b_[isplit+1],interval_data(isplit+1),res_r,err_r);
408  // update res-table
409  for(i=n_cur-1;i>isplit;i--){
410  results[i+1]=results[i];
411  errors[i+1]=errors[i];
412  }
413  results[isplit+1]=res_r;
414  errors[isplit+1]=err_r;
415  results[isplit]=res_l;
416  errors[isplit]=err_l;
417  }
418  }
419  private:
421 
422  cplx *data_;
424 
425  int size_;
427 
428  int tot_;
430 
431  int nmax_;
433 
434  int n_;
436 
437  double *a_;
439 
440  double *b_;
442 
443  int *len_;
444  cplx **f_;
445 };
446 
447 
448 
449 } // namespace
450 
451 
452 #endif // fourier
453 
std::complex< double > cplx
Definition: fourier.hpp:44
Class adft_func contains functions for computing accurate Fourier transforms by adaptive splitting o...
Definition: fourier.hpp:42
std::complex< double > cplx
Definition: fourier.cpp:11
void get_dftcorr_linear(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction term needed for linearly corrected DFT...
Definition: fourier.cpp:117
void sample(double w, double a, double b, function &fz, int n, int limit)
Fourier samples the given function to determine the individual intervals.
Definition: fourier.hpp:366
void dft_cplx(double w, int n, double a, double b, cplx *f, cplx &res, cplx &err)
Computes the Fourier integral using cubically corrected DFT.
Definition: fourier.cpp:283
void complex_dftcor_cubic(double w, double delta, double a, double b, cplx *endpts, cplx *endcor, double *corfac)
Returns correction factor and evaluate boundary correction terms needed for cubically corrected DFT...
Definition: fourier.cpp:182
void dft(double w, cplx &result, cplx &err)
Computes the cubically corrected DFT.
Definition: fourier.hpp:317
void get_dftcorr_cubic(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction terms needed for cubically corrected DFT...
Definition: fourier.cpp:43