Removed unused integrator code, cleaned up Makefile and .gitignore

7 files changed, 7 insertions, 626 deletions

diff --git a/.gitignore b/.gitignore
index 377d382..56f8120 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,13 +1,10 @@
 *.dat
 *.o
+vgcore.*
+fftw.wisdom
+walk
 correlations
-fourier_integrator
 log-fourier_integrator
 log-fourier_integrator_long
-vgcore.*
-walk
-get_energy
-integrator
-fftw.wisdom
 log_get_energy
 log_get_energy_long
diff --git a/Makefile b/Makefile
index adb9d3d..91081f5 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy log_get_energy_long
+all: walk correlations log-fourier_integrator log-fourier_integrator_long log_get_energy log_get_energy_long
 
 CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
 
@@ -8,36 +8,24 @@ walk: walk.cpp
 correlations: correlations.cpp
 	$(CC) correlations.cpp -o correlations
 
-integrator: integrator.cpp fourier.o p-spin.o fourier.hpp
-	$(CC) integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o integrator
-
-fourier.o: fourier.cpp fourier.hpp types.hpp p-spin.o
-	$(CC) fourier.cpp -c -o fourier.o
-
 p-spin.o: p-spin.cpp p-spin.hpp types.hpp
 	$(CC) p-spin.cpp -c -o p-spin.o
 
-log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp
-	$(CC) log-fourier.cpp -c -o log-fourier.o
-
 p-spin_long.o: p-spin.cpp p-spin.hpp types.hpp
 	$(CC) p-spin.cpp -D QUAD_PRECISION -c -o p-spin_long.o
 
+log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp
+	$(CC) log-fourier.cpp -c -o log-fourier.o
+
 log-fourier_long.o: log-fourier.cpp log-fourier.hpp types.hpp
 	$(CC) log-fourier.cpp -D QUAD_PRECISION -c -o log-fourier_long.o
 
-fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o
-	$(CC) fourier_integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o fourier_integrator
-
 log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o
 	$(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator
 
 log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
 	$(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long
 
-get_energy: get_energy.cpp fourier.hpp fourier.o p-spin.o
-	$(CC) get_energy.cpp fourier.o p-spin.o -lfftw3 -lgsl -o get_energy
-
 log_get_energy: log_get_energy.cpp log-fourier.hpp log-fourier.o p-spin.o
 	$(CC) log_get_energy.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log_get_energy
 
diff --git a/fourier.cpp b/fourier.cpp
deleted file mode 100644
index 3821623..0000000
--- a/fourier.cpp
+++ /dev/null
@@ -1,125 +0,0 @@
-#include "fourier.hpp"
-#include "p-spin.hpp"
-#include <fftw3.h>
-
-FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : n(n), Δω(Δω), Δτ(Δτ) {
-  a = fftw_alloc_real(2 * n);
-  â = reinterpret_cast<Complex*>(fftw_alloc_complex(n + 1));
-//  fftw_init_threads();
-//  fftw_plan_with_nthreads(FFTW_THREADS);
-  fftw_import_wisdom_from_filename("fftw.wisdom");
-  plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a, reinterpret_cast<fftw_complex*>(â), flags);
-  plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags);
-  fftw_export_wisdom_to_filename("fftw.wisdom");
-}
-
-FourierTransform::~FourierTransform() {
-  fftw_destroy_plan(plan_r2c);
-  fftw_destroy_plan(plan_c2r);
-  fftw_free(a);
-  fftw_free(â);
-  fftw_cleanup();
-}
-
-std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) {
-  for (unsigned i = 0; i < 2 * n; i++) {
-    a[i] = c[i];
-  }
-  fftw_execute(plan_r2c);
-  std::vector<Complex> ĉ(n + 1);
-  for (unsigned i = 0; i < n + 1; i++) {
-    ĉ[i] = â[i] * (Δτ * M_PI);
-  }
-  return ĉ;
-}
-
-std::vector<Complex> FourierTransform::fourier() {
-  fftw_execute(plan_r2c);
-  std::vector<Complex> ĉ(n+1);
-  for (unsigned i = 0; i < n+1; i++) {
-    ĉ[i] = â[i] * (Δτ * M_PI);
-  }
-  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];
-  }
-  fftw_execute(plan_c2r);
-  std::vector<Real> c(2*n);
-  for (unsigned i = 0; i < 2*n; i++) {
-    c[i] = a[i] * (Δω / (2 * M_PI));
-  }
-  return c;
-}
-
-void FourierTransform::writeToA(unsigned i, Real ai) {
-  a[i] = ai;
-}
-
-std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) {
-  return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
-}
-
-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 n = 0; n < C.size() / 4 -1; n++) {
-    Real h₂ₙ   = M_PI * Δτ;
-    Real h₂ₙ₊₁ = M_PI * Δτ;
-    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 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 λ) {
-  for (unsigned i = 0; i < C.size() / 2; i++) {
-    fft.writeToA(i,  R[i] * ddf(λ, p, s, C[i]));
-  }
-  for (unsigned i = C.size() / 2; i < C.size(); i++) {
-    fft.writeToA(i, 0);
-  }
-  std::vector<Complex> RddfCt = fft.fourier();
-
-  for (unsigned i = 0; i < C.size(); i++) {
-    fft.writeToA(i, df(λ, p, s, C[i]));
-  }
-  std::vector<Complex> dfCt = fft.fourier();
-
-  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 τ₀, Real y) {
-  auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-  Real Γ₀ = 1 + τ₀ / 2;
-
-  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
deleted file mode 100644
index 6be0931..0000000
--- a/fourier.hpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#pragma once
-#include "types.hpp"
-
-#include <cmath>
-#include <fftw3.h>
-#include <vector>
-#include <tuple>
-
-#ifndef FFTW_THREADS
-#define FFTW_THREADS 1
-#endif
-
-class FourierTransform {
-private:
-  Real* a;
-  Complex* â;
-  fftw_plan plan_r2c;
-  fftw_plan plan_c2r;
-  unsigned n;
-  Real Δω;
-  Real Δτ;
-public:
-  FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0);
-  ~FourierTransform();
-  std::vector<Complex> fourier(const std::vector<Real>& c);
-  std::vector<Complex> fourier();
-  std::vector<Real> inverse(const std::vector<Complex>& ĉ);
-  void writeToA(unsigned i, Real ai);
-  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
-};
-
-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 τ₀, Real y);
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
deleted file mode 100644
index 96fbffd..0000000
--- a/fourier_integrator.cpp
+++ /dev/null
@@ -1,213 +0,0 @@
-#include "fourier.hpp"
-#include <getopt.h>
-#include <iostream>
-#include <fstream>
-
-int main(int argc, char* argv[]) {
-  unsigned p = 2;
-  unsigned s = 2;
-  Real λ = 0.5;
-  Real τ₀ = 0;
-  Real β₀ = 0;
-  Real βₘₐₓ = 0.5;
-  Real Δβ = 0.05;
-
-  unsigned log2n = 8;
-  Real τₘₐₓ = 20;
-
-  unsigned maxIterations = 1000;
-  Real ε = 1e-14;
-  Real γ = 1;
-
-  bool loadData = false;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
-    switch (opt) {
-    case 'p':
-      p = atoi(optarg);
-      break;
-    case 's':
-      s = atoi(optarg);
-      break;
-    case '2':
-      log2n = atoi(optarg);
-      break;
-    case 'T':
-      τₘₐₓ = atof(optarg);
-      break;
-    case 't':
-      τ₀ = atof(optarg);
-      break;
-    case '0':
-      β₀ = atof(optarg);
-      break;
-    case 'y':
-      βₘₐₓ = atof(optarg);
-      break;
-    case 'd':
-      Δβ = atof(optarg);
-      break;
-    case 'I':
-      maxIterations = (unsigned)atof(optarg);
-      break;
-    case 'g':
-      γ = atof(optarg);
-      break;
-    case 'l':
-      loadData = true;
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  unsigned n = pow(2, log2n);
-
-  Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
-  Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
-
-  Real Γ₀ = 1.0;
-  Real μ = Γ₀;
-  if (τ₀ > 0) {
-    μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
-  }
-
-  std::vector<Real> C(2 * n);
-  std::vector<Real> R(2 * n);
-
-  FourierTransform fft(n, Δω, Δτ);
-  std::vector<Complex> Ct;
-  std::vector<Complex> Rt;
-
-  Real β = β₀;
-
-  if (!loadData) {
-    // start from the exact solution for τ₀ = 0
-    for (unsigned i = 0; i < n; i++) {
-      Real τ = i * Δτ * M_PI;
-      if (τ₀ > 0) {
-        C[i] = Γ₀ * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
-      } else {
-        C[i] = Γ₀ * exp(-μ * τ) / μ;
-      }
-      if (i > 0) {
-        C[2 * n - i] = C[i];
-      }
-      R[i] = exp(-μ * τ);
-    }
-    Ct = fft.fourier(C);
-    Rt = fft.fourier(R);
-  } else {
-    std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
-    cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
-    cfile.close();
-    for (unsigned i = 1; i < n; i++) {
-      C[2 * n - i] = C[i];
-    }
-    std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
-    rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
-    rfile.close();
-
-    Ct = fft.fourier(C);
-    Rt = fft.fourier(R);
-
-    μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);
-  }
-
-  std::vector<Real> Cb = C;
-  std::vector<Real> Rb = R;
-  std::vector<Complex> Ctb = Ct;
-  std::vector<Complex> Rtb = Rt;
-  Real μb = μ;
-
-  Real fac = 1;
-  while (β <= βₘₐₓ) {
-    Real ΔC = 1;
-    Real ΔCprev = 1000;
-    unsigned it = 0;
-    while (sqrt(2 * ΔC / C.size()) > ε) {
-      it++;
-      auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
-      for (unsigned i = 0; i < Rt.size(); i++) {
-        Real ω = i * Δω;
-        Rt[i] = (1.0 + pow(β, 2) * RddfCt[i] * Rt[i]) / (μ + 1i * ω);
-      }
-
-      for (unsigned i = 0; i < Ct.size(); i++) {
-        Real ω = i * Δω;
-        Ct[i] = (2 * Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(β, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (μ + 1i * ω);
-      }
-
-      std::vector<Real> Cnew = fft.inverse(Ct);
-      std::vector<Real> Rnew = fft.inverse(Rt);
-      for (unsigned i = n; i < 2 * n; i++) {
-        Rnew[i] = 0;
-      }
-
-      μ *= pow(tanh(C[0]-1)+1, 0.05);
-
-      ΔC = 0;
-      for (unsigned i = 0; i < Cnew.size() / 2; i++) {
-        ΔC += pow(Cnew[i] - C[i], 2);
-        ΔC += pow(Rnew[i] - R[i], 2);
-      }
-
-      for (unsigned i = 0; i < Cnew.size(); i++) {
-        C[i] += γ * (Cnew[i] - C[i]);
-      }
-
-      for (unsigned i = 0; i < Rnew.size() / 2; i++) {
-        R[i] += γ * (Rnew[i] - R[i]);
-      }
-
-      if (it % maxIterations == 0) {
-        if (ΔC < ΔCprev) {
-          γ = std::min(1.1 * γ, 1.0);
-        } else {
-          γ /= 2;
-        }
-
-        ΔCprev = ΔC;
-      }
-
-      std::cerr << "\x1b[2K" << "\r";
-      std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];
-    }
-
-    if (std::isnan(C[0])) {
-      Δβ /= 2;
-      β -= Δβ;
-      C = Cb;
-      R = Rb;
-      Ct = Ctb;
-      Rt = Rtb;
-      μ = μb;
-    } else {
-
-    std::cerr << "\x1b[2K" << "\r";
-
-    Real e = energy(C, R, p, s, λ, β, Δτ);
-
-    std::cerr << "y " << β << " " << e << " " << μ << std::endl;
-
-    std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real));
-    outfile.close();
-
-    std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
-    outfileR.close();
-    β += Δβ;
-    Cb = C;
-    Rb = R;
-    Rtb = Rt;
-    Ctb = Ct;
-    μb = μ;
-    }
-  }
-
-  return 0;
-}
diff --git a/get_energy.cpp b/get_energy.cpp
deleted file mode 100644
index deee26b..0000000
--- a/get_energy.cpp
+++ /dev/null
@@ -1,91 +0,0 @@
-#include "fourier.hpp"
-#include <getopt.h>
-#include <iostream>
-#include <fstream>
-#include <iomanip>
-
-int main(int argc, char* argv[]) {
-  unsigned p = 2;
-  unsigned s = 2;
-  Real λ = 0.5;
-  Real τ₀ = 0;
-  Real y₀ = 0;
-  Real yₘₐₓ = 0.5;
-  Real Δy = 0.05;
-
-  unsigned log2n = 8;
-  Real τₘₐₓ = 20;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -1) {
-    switch (opt) {
-    case 'p':
-      p = atoi(optarg);
-      break;
-    case 's':
-      s = atoi(optarg);
-      break;
-    case '2':
-      log2n = atoi(optarg);
-      break;
-    case 'T':
-      τₘₐₓ = atof(optarg);
-      break;
-    case 't':
-      τ₀ = atof(optarg);
-      break;
-    case '0':
-      y₀ = atof(optarg);
-      break;
-    case 'y':
-      yₘₐₓ = atof(optarg);
-      break;
-    case 'd':
-      Δy = atof(optarg);
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  unsigned n = pow(2, log2n);
-
-  Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
-  Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
-
-  Real y = y₀;
-
-  FourierTransform fft(n, Δω, Δτ);
-
-  std::vector<Real> C(2 * n);
-  std::vector<Real> R(2 * n);
-
-  std::cout << std::setprecision(15);
-
-  while (y = std::round(1e6 * (y + Δy)) / 1e6, y <= yₘₐₓ) {
-    std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
-    if (cfile.is_open()) {
-      cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
-      cfile.close();
-      for (unsigned i = 1; i < n; i++) {
-        C[2 * n - i] = C[i];
-      }
-
-      std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
-      rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
-      rfile.close();
-
-      Real e = energy(C, R, p, s, λ, y, Δτ);
-
-      std::vector<Complex> Ct = fft.fourier(C);
-      std::vector<Complex> Rt = fft.fourier(R);
-
-      Real z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
-
-      std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl;
-    }
-  }
-
-  return 0;
-}
diff --git a/integrator.cpp b/integrator.cpp
deleted file mode 100644
index 4f8246d..0000000
--- a/integrator.cpp
+++ /dev/null
@@ -1,140 +0,0 @@
-#include "fourier.hpp"
-#include "p-spin.hpp"
-#include <fstream>
-#include <fftw3.h>
-#include <getopt.h>
-#include <iostream>
-
-Real energy(const std::vector<Real>& C, std::vector<Real>& R, Real λ, unsigned p, unsigned s, Real Δτ) {
-  Real I = 0;
-  for (unsigned σ = 0; σ < C.size(); σ++) {
-    I += Δτ * df(λ, p, s, C[σ]) * R[σ];
-  }
-  return I;
-}
-
-int main(int argc, char* argv[]) {
-  unsigned p = 3;
-  unsigned s = 4;
-  Real λ = 0.5;
-  Real τₘₐₓ = 1e3;
-  Real τ₀ = 0;
-  Real β₀ = 0;
-  Real βₘₐₓ = 1;
-  Real Δβ = 1e-2;
-  unsigned iterations = 10;
-  unsigned log2n = 8;
-  Real ε = 1e-14;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "T:2:t:0:b:d:I:")) != -1) {
-    switch (opt) {
-    case 'T':
-      τₘₐₓ = atof(optarg);
-      break;
-    case '2':
-      log2n = atof(optarg);
-      break;
-    case 't':
-      τ₀ = atof(optarg);
-      break;
-    case '0':
-      β₀ = atof(optarg);
-      break;
-    case 'b':
-      βₘₐₓ = atof(optarg);
-      break;
-    case 'd':
-      Δβ = atof(optarg);
-      break;
-    case 'I':
-      iterations = (unsigned)atof(optarg);
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  unsigned N = pow(2, log2n);
-
-  Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / N;
-  Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
-
-  FourierTransform fft(N, Δω, Δτ, FFTW_ESTIMATE);
-
-  Real Γ₀ = 1;
-  Real μ = 1;
-  if (τ₀ > 0) {
-    μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀);
-  }
-
-  Real τ = 0;
-  std::vector<Real> C(N);
-  std::vector<Real> R(N);
-  std::vector<Real> Γ(N);
-  std::vector<Real> Γh(N+1);
-
-  Γh[0] = Γ₀;
-
-  for (unsigned i = 0; i < N; i++) {
-    Real τ = i * Δτ;
-    Real ω = (i + 1) * Δω * M_PI;
-    if (τ₀ > 0) {
-      C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2));
-      Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀);
-    } else {
-      C[i] = (Γ₀ / μ) * exp(-μ * τ);
-    }
-    Γh[i+1] = Γ₀ / (1 + pow(ω * τ₀, 2));
-    R[i] = exp(-μ * τ);
-  }
-
-  for (Real β = β₀; β < βₘₐₓ; β += Δβ) {
-    Real Rerr = 100;
-    while (sqrt(Rerr / N) > ε) {
-      /* First step: integrate R from C */
-      std::vector<Real> R₊(N);
-      R₊[0] = 1;
-      for (unsigned i = 1; i < N; i++) {
-        Real I = 0;
-        for (unsigned j = 0; j <= i; j++) {
-          I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ;
-        }
-        Real dR = -μ * R₊[i - 1] + pow(β, 2) * I;
-        R₊[i] = R₊[i - 1] + dR * Δτ;
-      }
-
-      Rerr = 0;
-      for (unsigned i = 0; i < N; i++) {
-        Rerr += pow(R[i] - R₊[i], 2);
-      }
-
-      R = R₊;
-
-      /* Second step: integrate C from R */
-      std::vector<Real> dC = fft.convolve(Γh, R);
-      Real Cₜ₊₁ = 0;
-      for (unsigned i = 0; i < N; i++) {
-        Real Cₜ = Cₜ₊₁ + dC[N - i - 1] * Δτ;
-        C[N - i - 1] = Cₜ;
-        Cₜ₊₁ = Cₜ;
-      }
-
-      /* Third step: adjust μ */
-      μ *= C[0];
-
-      std::cerr << β << " " << sqrt(Rerr / N) << std::endl;
-    }
-
-    std::ofstream outfile(fourierFile("Ci", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfile.write((const char*)(C.data()), N * sizeof(Real));
-    outfile.close();
-
-    std::ofstream outfileR(fourierFile("Ri", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
-    outfileR.write((const char*)(R.data()), N * sizeof(Real));
-    outfileR.close();
-  }
-
-  return 0;
-}