Rewrote the original integrator to use a new approach.HEADmaster

1 files changed, 19 insertions, 0 deletions

diff --git a/fourier.cpp b/fourier.cpp
index 07f8fdc..bc4b633 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -64,6 +64,25 @@ std::vector<Complex> FourierTransform::fourier() {
   return ĉ;
 }
 
+std::vector<Real> FourierTransform::convolve(const std::vector<Real>& Γh, const std::vector<Real>& R) {
+  a[0] = 0;
+  for (unsigned i = 1; i < n; i++) {
+    a[i] = R[i];
+    a[2 * n - i] = -R[i];
+  }
+  fftw_execute(plan_r2c);
+  for (unsigned i = 1; i < n + 1; i++) {
+    â[i] *= Γh[i] * (Δτ * M_PI);
+  }
+  fftw_execute(plan_c2r);
+  std::vector<Real> dC(n);
+  for (unsigned i = 0; i < n; i++) {
+    dC[i] = a[i] * (Δω / (2 * M_PI));
+  }
+
+  return dC;
+}
+
 std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) {
   for (unsigned i = 0; i < n + 1; i++) {
     â[i] = ĉ[i];