Removed some more duplicate code.

4 files changed, 10 insertions, 7 deletions

diff --git a/fourier.cpp b/fourier.cpp
index 7808989..8943521 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -71,7 +71,7 @@ Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p,
 
 std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) {
   std::vector<Real> RddfC(C.size());
-  for (unsigned i = 0; i < C.size(); i++) {
+  for (unsigned i = 0; i < C.size() / 2; i++) {
     RddfC[i] = R[i] * ddf(λ, p, s, C[i]);
   }
   std::vector<Complex> RddfCt = fft.fourier(RddfC);
@@ -84,3 +84,9 @@ std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransfo
 
   return {RddfCt, dfCt};
 }
+
+Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real y) {
+  auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
+
+  return ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+}
diff --git a/fourier.hpp b/fourier.hpp
index 86335f7..3925058 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -30,3 +30,4 @@ public:
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
 Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ);
 std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);
+Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real y);
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 86c3b24..2fe7e49 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -103,9 +103,7 @@ int main(int argc, char* argv[]) {
     Ct = fft.fourier(C);
     Rt = fft.fourier(R);
 
-    auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
-    z = ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+    z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, y);
   }
 
   while (y += Δy, y <= yₘₐₓ) {
diff --git a/get_energy.cpp b/get_energy.cpp
index a388e0c..b1527dc 100644
--- a/get_energy.cpp
+++ b/get_energy.cpp
@@ -75,9 +75,7 @@ int main(int argc, char* argv[]) {
       std::vector<Complex> Ct = fft.fourier(C);
       std::vector<Complex> Rt = fft.fourier(R);
 
-      auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
-      Real z = ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+      Real z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, y);
 
       std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl;
     }