Force monotonicity of C and R

1 files changed, 16 insertions, 1 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 04dccfd..2d7d1f9 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -145,7 +145,22 @@ int main(int argc, char* argv[]) {
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-      Real C₀ = (Cₜ₊₁[0] + Cₜ₊₁[1]) / 2;
+      Real Cmin = Cₜ₊₁[N-1];
+      for (unsigned i = 0; i < N; i++) {
+        if (Cₜ₊₁[N-1-i] > Cmin) {
+          Cmin = Cₜ₊₁[N-1-i];
+        }
+        Cₜ₊₁[N-1-i] = Cmin;
+      }
+      Real Rmin = Rₜ₊₁[N-1];
+      for (unsigned i = 0; i < N; i++) {
+        if (Rₜ₊₁[N-1-i] > Rmin) {
+          Rmin = Rₜ₊₁[N-1-i];
+        }
+        Rₜ₊₁[N-1-i] = Rmin;
+      }
+
+      Real C₀ = Cₜ₊₁[0];
 
       μₜ *= pow(tanh(C₀-1)+1, x);