NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library

◆ get_dftcorr_cubic()

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

Returns correction factor \(W(\theta)\) and boundary correction terms \(\alpha_j(\theta)\) 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

get_dftcorr_cubic computes the required correction factor \(W(\theta)\) and the boundary terms \(\alpha_i(\theta)\), \(i=0,\dots,3\), 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_j(\theta)\), stored in 'endcor[j]'

Definition at line 43 of file fourier.cpp.

Referenced by cntr::matsubara_ft().

44 {
45  double ai[4],ar[4],t;
46  double t2,t4,t6;
47  double cth,ctth,spth2,sth,sth4i,stth,th2,th4,tmth2,tth4i;
48  int j;
49 
50  if (fabs(th) < 5.0e-2) {
51  t=th;
52  t2=t*t;
53  t4=t2*t2;
54  t6=t4*t2;
55  *corfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6;
56  ar[0]=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6;
57  ar[1]=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6;
58  ar[2]=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0;
59  ar[3]=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0;
60  ai[0]=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6);
61  ai[1]=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6);
62  ai[2]=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6);
63  ai[3]=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6);
64  } else {
65  cth=cos(th);
66  sth=sin(th);
67  ctth=cth*cth-sth*sth;
68  stth=2.0e0*sth*cth;
69  th2=th*th;
70  th4=th2*th2;
71  tmth2=3.0e0-th2;
72  spth2=6.0e0+th2;
73  sth4i=1.0/(6.0e0*th4);
74  tth4i=2.0e0*sth4i;
75  *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth);
76  ar[0]=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth));
77  ai[0]=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth);
78  ar[1]=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth);
79  ai[1]=sth4i*(30.0e0*th-5.0e0*spth2*sth);
80  ar[2]=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth);
81  ai[2]=tth4i*(-12.0e0*th+2.0e0*spth2*sth);
82  ar[3]=sth4i*(2.0e0*tmth2-spth2*cth);
83  ai[3]=sth4i*(6.0e0*th-spth2*sth);
84  }
85  for(j=0;j<=3;j++) endcor[j]=cplx(ar[j],ai[j]);
86 }
std::complex< double > cplx
Definition: fourier.cpp:11
+ Here is the caller graph for this function: