NESSi  v1.0.2 The NonEquilibrium Systems Simulation Library

◆ matsubara_ft()

template<typename T , class GG , int SIZE1>
 void cntr::matsubara_ft ( std::complex< T > * result, int m, GG & G, std::complex< T > * mdft, int sig, T beta, int order )

Computes the Fourier series coefficients of a Matsubara function by cubically corrected DFT.

Purpose

Evaluates the Fourier integral $$G^\mathrm{M}(i \omega_m) = \int^\beta_0 d\tau\,e^{i \omega_m \tau} G^\mathrm{M}(\tau)$$ by linearly or cubically corrected DFT.

Parameters
 result [complex] on return, pointer storing the Fourier coeffient $$G^\mathrm{M}(i \omega_m)$$ m [int] index of Matsubara frequency $$\omega_m$$ G [GG] Matsubara function to be Fourier transformed mdft [complex] pointer storing the Fourier coeffients computed by matsubara_dft in element representation sig [int] sig=-1 for fermions, sig=+1 for bosons beta [double] inverse temperature order [int] order=1 / order=3 for linearly/cubically corrected DFT

Definition at line 191 of file cntr_matsubara_impl.hpp.

References fourier::get_dftcorr_cubic(), and fourier::get_dftcorr_linear().

192  {
193  typedef std::complex<T> cplx;
194  double corfac, arg, dtau, one;
195  int ntau, m1, r, l, sg, size1 = G.size1();
196  cplx *z1, *z2, *dft, bcorr[4];
197
198  one = (sig == -1 ? 1.0 : 0.0);
199  ntau = G.ntau();
200  sg = G.element_size();
201  dtau = beta / ntau;
202  z1 = new cplx[sg];
203  z2 = new cplx[sg];
204  dft = new cplx[sg];
205
206  arg = (2 * m + one) * PI / ntau; /*arg=omn*dtau, omn=(2m+1)Pi/beta*/
207  if (order == 0) {
208  corfac = 1.0;
209  bcorr[0] = 0.0;
210  } else if (order == 1) {
211  fourier::get_dftcorr_linear(arg, &corfac, bcorr);
212  } else if (order == 3) {
213  fourier::get_dftcorr_cubic(arg, &corfac, bcorr);
214  } else {
215  std::cerr << "matsubara_ft: unknown order" << std::endl;
216  abort();
217  }
218  m1 = m; /*shift m1 to interval [0,ntau-1]*/
219  while (m1 < 0)
220  m1 += ntau;
221  while (m1 >= ntau)
222  m1 -= ntau;
223  for (l = 0; l < sg; l++)
224  dft[l] = mdft[m1 * sg + l] * corfac;
225  element_set_zero<T, SIZE1>(size1, z1);
226  for (r = 0; r <= order; r++) {
227  element_set<T, SIZE1>(size1, z2, G.matptr(r));
228  element_smul<T, SIZE1>(size1, z2, bcorr[r]);
229  element_incr<T, SIZE1>(size1, z1, z2);
230  element_set<T, SIZE1>(size1, z2, G.matptr(ntau - r));
231  // element_smul<T,SIZE1>(size1,z2,-cplx(bcorr[r].real(),-bcorr[r].imag()));
232  element_smul<T, SIZE1>(size1, z2,
233  (1. * sig) * cplx(bcorr[r].real(), -bcorr[r].imag()));
234  element_incr<T, SIZE1>(size1, z1, z2);
235  }
236  for (l = 0; l < sg; l++)
237  result[l] = dtau * (z1[l] + dft[l]);
238
239  delete[] z1;
240  delete[] z2;
241  delete[] dft;
242 }
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 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
Here is the call graph for this function: