1 #ifndef CNTR_MATSUBARA_IMPL_H 2 #define CNTR_MATSUBARA_IMPL_H 14 template <
typename F,
typename I> F get_tau(
const I tau_idx,
const F beta,
const I ntau) {
15 return tau_idx * (beta / ntau);
19 template <
typename F,
typename I> F get_omega(
const I m,
const F beta,
const I sig) {
23 const I matsub_one = ((1 - sig) >> 1);
27 return PI / beta * (2 * m + matsub_one);
70 template <
typename T,
int SIZE1>
71 void set_first_order_tail(std::complex<T> *xmat, std::complex<T> *coeff, T beta,
int sg,
72 int ntau,
int sig,
int size1) {
78 std::complex<T> *z1 =
new std::complex<T>[sg];
80 for (
int r = 0; r <= ntau; r++) {
82 element_set<T, SIZE1>(size1, xmat + r * sg, coeff);
83 element_smul<T, SIZE1>(size1, xmat + r * sg, -0.5);
87 double tau = r * (beta / ntau);
88 element_set<T, SIZE1>(size1, z1, coeff);
89 element_smul<T, SIZE1>(size1, z1, tau / beta);
90 element_incr<T, SIZE1>(size1, xmat + r * sg, z1);
133 template <
typename T,
class GG,
int SIZE1>
135 typedef std::complex<T>
cplx;
136 int ntau, r, m, sg, size1 = G.size1();
138 cplx *z, *z1, expfac;
139 sg = G.element_size();
143 one = (sig == -1 ? 1.0 : 0.0);
144 for (m = 0; m <= ntau; m++) {
145 element_set_zero<T, SIZE1>(size1, z);
146 for (r = 0; r <= ntau; r++) {
147 arg = ((2 * m + one) * r * PI) / ntau;
148 expfac =
cplx(cos(arg), sin(arg));
149 element_set<T, SIZE1>(size1, z1, G.matptr(r));
150 element_smul<T, SIZE1>(size1, z1, expfac);
151 element_incr<T, SIZE1>(size1, z, z1);
153 element_set<T, SIZE1>(size1, mdft + m * sg, z);
190 template <
typename T,
class GG,
int SIZE1>
191 void matsubara_ft(std::complex<T> *result,
int m, GG &G, std::complex<T> *mdft,
int sig,
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];
198 one = (sig == -1 ? 1.0 : 0.0);
200 sg = G.element_size();
206 arg = (2 * m + one) * PI / ntau;
210 }
else if (order == 1) {
212 }
else if (order == 3) {
215 std::cerr <<
"matsubara_ft: unknown order" << std::endl;
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));
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);
236 for (l = 0; l < sg; l++)
237 result[l] = dtau * (z1[l] + dft[l]);
248 #endif // CNTR_MATSUBARA_IMPL_H
std::complex< double > cplx
Integrator< T > & I(int k)
void get_dftcorr_linear(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction term needed for linearly corrected DFT...
void matsubara_dft(std::complex< T > *mdft, GG &G, int sig)
Computes the Fourier series coefficients of a Matsubara function by plain DFT.
void 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...
void get_dftcorr_cubic(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction terms needed for cubically corrected DFT...