NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library

◆ complex_dftcor_cubic()

void fourier::complex_dftcor_cubic ( double  w,
double  delta,
double  a,
double  b,
cplx endpts,
cplx endcor,
double *  corfac 
)

Returns correction factor \(W(\theta)\) and evaluate boundary correction terms needed for cubically corrected DFT.

Purpose

For a function \(f(t)\), represented on an equidistant grid \(t_n= n h\), \(n=0,\dots,N\), we would like to compute the Fourier integral \(I(\omega) = \int^b_a dt\, e^{i \omega t}f(t)\). Using piecewie cubic interpolation of \(f(t)\), the \(I(\omega)\) can be computed by by cubically corrected discrete Fourier transformation (DFT). The algorithm is explained in

W. H. Press, S. A. Teukolosky, W. T. Vetterling, B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing, Cambridge University Press, 2007, chapter 13.9

complex_dftcor_cubic computes the required correction factor \(W(\theta)\) and the boundary terms \(\alpha_j(\theta)f_j\), \(j=0,\dots,3\), and \( e^{i\omega(b-a)}\alpha^*_j(\theta)f_{N-j} \) where \(\theta = \omega h\).

Parameters
w

[double] frequency \(\omega\)

delta

[double] grid spacing \(h\)

a

[double] starting point of interval

b

[double] end point of interval

endpts

[complex] 'endpts[j]', 'j=0,1,2,3', corresponds to \(f_j\), while 'endpts[7-j]', 'j=0,1,2,3', corresponds to \(f_{N-j}\)

endcor

[complex] on return, boundary correction terms, stored in the order as 'endpts'

corfac

[double] on return, correction factor \(W(\theta)\)

Definition at line 182 of file fourier.cpp.

Referenced by dft_cplx().

183 {
184  double ai[4],ar[4],arg,t;
185  double t2,t4,t6;
186  double cth,ctth,spth2,sth,sth4i,stth,th,th2,th4,tmth2,tth4i;
187  cplx expfac,cr,cl,cor,al;
188  int j;
189 
190  th=w*delta;
191  if (fabs(th) < 5.0e-2) {
192  t=th;
193  t2=t*t;
194  t4=t2*t2;
195  t6=t4*t2;
196  *corfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6;
197  ar[0]=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6;
198  ar[1]=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6;
199  ar[2]=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0;
200  ar[3]=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0;
201  ai[0]=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6);
202  ai[1]=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6);
203  ai[2]=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6);
204  ai[3]=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6);
205  } else {
206  cth=cos(th);
207  sth=sin(th);
208  ctth=cth*cth-sth*sth;
209  stth=2.0e0*sth*cth;
210  th2=th*th;
211  th4=th2*th2;
212  tmth2=3.0e0-th2;
213  spth2=6.0e0+th2;
214  sth4i=1.0/(6.0e0*th4);
215  tth4i=2.0e0*sth4i;
216  *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth);
217  ar[0]=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth));
218  ai[0]=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth);
219  ar[1]=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth);
220  ai[1]=sth4i*(30.0e0*th-5.0e0*spth2*sth);
221  ar[2]=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth);
222  ai[2]=tth4i*(-12.0e0*th+2.0e0*spth2*sth);
223  ar[3]=sth4i*(2.0e0*tmth2-spth2*cth);
224  ai[3]=sth4i*(6.0e0*th-spth2*sth);
225  }
226  cl=0;
227  for(j=0;j<=3;j++){
228  al=cplx(ar[j],ai[j]);
229  cor=endpts[j]*al;
230  cl+=cor;
231  }
232  cr=0;
233  for(j=0;j<=3;j++){
234  al=cplx(ar[j],-ai[j]);
235  cor=endpts[7-j]*al;
236  cr += cor;
237  }
238  arg=w*(b-a);
239  expfac=cplx(cos(arg),sin(arg));
240  cr *= expfac;
241  *endcor = cl+cr;
242 }
std::complex< double > cplx
Definition: fourier.cpp:11
+ Here is the caller graph for this function: