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;