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().

                                     {
    typedef std::complex<T> cplx;
    double corfac, arg, dtau, one;
    int ntau, m1, r, l, sg, size1 = G.size1();
    cplx *z1, *z2, *dft, bcorr[4];

    one = (sig == -1 ? 1.0 : 0.0);
    ntau = G.ntau();
    sg = G.element_size();
    dtau = beta / ntau;
    z1 = new cplx[sg];
    z2 = new cplx[sg];
    dft = new cplx[sg];

    arg = (2 * m + one) * PI / ntau; /*arg=omn*dtau, omn=(2m+1)Pi/beta*/
    if (order == 0) {
        corfac = 1.0;
        bcorr[0] = 0.0;
    } else if (order == 1) {
        fourier::get_dftcorr_linear(arg, &corfac, bcorr);
    } else if (order == 3) {
        fourier::get_dftcorr_cubic(arg, &corfac, bcorr);
    } else {
        std::cerr << "matsubara_ft: unknown order" << std::endl;
        abort();
    }
    m1 = m; /*shift m1 to interval [0,ntau-1]*/
    while (m1 < 0)
        m1 += ntau;
    while (m1 >= ntau)
        m1 -= ntau;
    for (l = 0; l < sg; l++)
        dft[l] = mdft[m1 * sg + l] * corfac;
    element_set_zero<T, SIZE1>(size1, z1);
    for (r = 0; r <= order; r++) {
        element_set<T, SIZE1>(size1, z2, G.matptr(r));
        element_smul<T, SIZE1>(size1, z2, bcorr[r]);
        element_incr<T, SIZE1>(size1, z1, z2);
        element_set<T, SIZE1>(size1, z2, G.matptr(ntau - r));
        // element_smul<T,SIZE1>(size1,z2,-cplx(bcorr[r].real(),-bcorr[r].imag()));
        element_smul<T, SIZE1>(size1, z2,
                               (1. * sig) * cplx(bcorr[r].real(), -bcorr[r].imag()));
        element_incr<T, SIZE1>(size1, z1, z2);
    }
    for (l = 0; l < sg; l++)
        result[l] = dtau * (z1[l] + dft[l]);

    delete[] z1;
    delete[] z2;
    delete[] dft;
}

Here is the call graph for this function: