diff options
| -rw-r--r-- | get_energy.cpp | 48 |
1 files changed, 47 insertions, 1 deletions
diff --git a/get_energy.cpp b/get_energy.cpp index 503823d..bb01095 100644 --- a/get_energy.cpp +++ b/get_energy.cpp @@ -34,6 +34,47 @@ inline Real ddf(Real λ, unsigned p, unsigned s, Real q) { return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q); } +class FourierTransform { +private: + std::vector<Real> a; + std::vector<Complex> â; + fftw_plan plan_r2c; + fftw_plan plan_c2r; + Real Δω; + Real Δτ; +public: + FourierTransform(unsigned n, Real Δω, Real Δτ) : a(2 * n), â(n + 1), Δω(Δω), Δτ(Δτ) { + plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a.data(), reinterpret_cast<fftw_complex*>(â.data()), 0); + plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â.data()), a.data(), 0); + } + + ~FourierTransform() { + fftw_destroy_plan(plan_r2c); + fftw_destroy_plan(plan_c2r); + fftw_cleanup(); + } + + std::vector<Complex> fourier(const std::vector<Real>& c) { + a = c; + fftw_execute(plan_r2c); + std::vector<Complex> ĉ(â.size()); + for (unsigned i = 0; i < â.size(); i++) { + ĉ[i] = â[i] * (Δτ * M_PI); + } + return ĉ; + } + + std::vector<Real> inverse(const std::vector<Complex>& ĉ) { + â = ĉ; + fftw_execute(plan_c2r); + std::vector<Real> c(a.size()); + for (unsigned i = 0; i < a.size(); i++) { + c[i] = a[i] * (Δω / (2 * M_PI)); + } + return c; + } +}; + int main(int argc, char* argv[]) { unsigned p = 2; unsigned s = 2; @@ -82,9 +123,12 @@ int main(int argc, char* argv[]) { unsigned n = pow(2, log2n); Real Δτ = τₘₐₓ / M_PI / n; + Real Δω = M_PI / τₘₐₓ; Real y = y₀; + FourierTransform fft(n, Δω, Δτ); + while (y += Δy, y <= yₘₐₓ) { std::vector<Real> C(2 * n); std::vector<Real> R(2 * n); @@ -103,7 +147,9 @@ int main(int argc, char* argv[]) { energy += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ; } - std::cerr << y << " " << energy << std::endl; + std::vector<Complex> Ct = fft.fourier(C); + + std::cerr << y << " " << energy << " " << Ct[0].real() << std::endl; } return 0; |