#include "fourier.hpp"
#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 Δτ = 1e-3;
  Real τₘₐₓ = 1e3;
  Real τ₀ = 0;
  Real β = 0.5;
  unsigned iterations = 10;

  int opt;

  while ((opt = getopt(argc, argv, "d:T:t:b:I:")) != -1) {
    switch (opt) {
    case 'd':
      Δτ = atof(optarg);
      break;
    case 'T':
      τₘₐₓ = atof(optarg);
      break;
    case 't':
      τ₀ = atof(optarg);
      break;
    case 'b':
      β = atof(optarg);
      break;
    case 'I':
      iterations = (unsigned)atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  Real Γ₀ = 1;
  Real μ = 1;
  if (τ₀ > 0) {
    μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀);
  }

  Real τ = 0;
  unsigned N = τₘₐₓ / Δτ + 1;
  std::vector<Real> C(N);
  std::vector<Real> R(N);
  std::vector<Real> Γ(N);

  for (unsigned i = 0; i < N; i++) {
    Real τ = i * Δτ;
    if (τ₀ > 0) {
      C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2));
      Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀);
    } else {
      C[i] = (Γ₀ / μ) * exp(-μ * τ);
    }
    R[i] = exp(-μ * τ);
  }

  for (unsigned it = 0; it < iterations; it++) {
    /* 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] + pow(β, 2) * I;
      R₊[i] = R₊[i - 1] + dR * Δτ;
    }

    /* Second step: integrate C from R */
  }

  τ = 0;
  for (Real Ci : C) {
    std::cout << τ << " " << Ci << std::endl;
    τ += Δτ;
  }

  std::cerr << - 2 * β / Γ₀ * energy(C, R, λ, p, s, Δτ) << std::endl;

  return 0;
}