NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library

◆ get_dftcorr_linear()

void fourier::get_dftcorr_linear ( double  th,
double *  corfac,
cplx endcor 
)

Returns correction factor \(W(\theta)\) and boundary correction term \(\alpha_0(\theta)\) needed for linearly 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 linear interpolation of \(f(t)\), the \(I(\omega)\) can be computed by by linearlt 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

get_dftcorr_linear computes the required correction factor \(W(\theta)\) and the boundary term \(\alpha_0(\theta)\), where \(\theta = \omega h\).

Parameters
th

[double] \(\theta = \omega h\) as explained above.

corfac

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

endcor

[complex] on return, boundary term \(\alpha_0(\theta)\)

Definition at line 117 of file fourier.cpp.

Referenced by cntr::matsubara_ft().

118 {
119  double ai,ar;
120  double th2,th4,th6,cth,sth;
121 
122  if (fabs(th) < 5.0e-2) {
123  th2=th*th;
124  th4=th2*th2;
125  th6=th4*th2;
126  *corfac=1.0-(1.0/12.0)*th2+(1.0/360.0)*th4-(1.0/20160.0)*th6;
127  ar=-0.5+th2/24.0-(1.0/720.0)*th4+(1.0/40320.0)*th6;
128  ai=th*(1.0/6.0-(1.0/120.0)*th2+(1.0/5040.0)*th4-(1.0/362880.0)*th6);
129  } else {
130  cth=cos(th);
131  sth=sin(th);
132  th2=th*th;
133  *corfac=2.0*(1.0-cth)/th2;
134  ar=(cth-1.0)/th2;
135  ai=(th-sth)/th2;
136  }
137  endcor[0]=cplx(ar,ai);
138  endcor[1]=0;
139 }
std::complex< double > cplx
Definition: fourier.cpp:11
+ Here is the caller graph for this function: