From 0af6d40ac9de7214458222e22c5566f95aa98ce4 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 19 Apr 2025 20:57:33 -0300
Subject: Use Simpson integrator for the regular fourier case too
---
fourier.cpp | 15 ++++++++++++---
1 file changed, 12 insertions(+), 3 deletions(-)
diff --git a/fourier.cpp b/fourier.cpp
index 95d0025..d637985 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -84,11 +84,20 @@ std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Rea
Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ) {
Real e = 0;
- for (unsigned i = 0; i < C.size() / 2; i++) {
- e += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
+ for (unsigned n = 0; n < C.size() / 4 -1; n++) {
+ Real h₂ₙ = Δτ;
+ Real h₂ₙ₊₁ = Δτ;
+ Real f₂ₙ = R[2*n] * df(λ, p, s, C[2*n]);
+ Real f₂ₙ₊₁ = R[2*n+1] * df(λ, p, s, C[2*n+1]);
+ Real f₂ₙ₊₂ = R[2*n+2] * df(λ, p, s, C[2*n+2]);
+ e += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
}
- return e;
+ return y * e;
}
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 λ) {
--
cgit v1.2.3-70-g09d2