#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 yₘₐₓ = 0.5;
  Real Δy = 0.05;

  unsigned log2n = 8;
  Real τₘₐₓ = 20;

  unsigned maxIterations = 1000;
  Real ε = 1e-14;
  Real γ = 1;

  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':
      yₘₐₓ = atof(optarg);
      break;
    case 'd':
      Δy = atof(optarg);
      break;
    case 'I':
      maxIterations = (unsigned)atof(optarg);
      break;
    case 'g':
      γ = atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  unsigned N = pow(2, log2n);

  Real Δτ = 1e-2;
  Real Δω = 1e-2;
  Real Δs = 1e-2;
  Real k = 0.1;

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

  std::vector<Real> C(N);
  std::vector<Real> R(N);

  LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs);

  for (unsigned n = 0; n < N; n++) {
    if (τ₀ > 0) {
      C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
    } else {
      C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
    }
    R[n] = exp(-μ * fft.t(n));
  }

  std::vector<Complex> Ct = fft.fourier(C, true);
  std::vector<Complex> Rt = fft.fourier(R, false);

  /*
  for (unsigned n = 0; n < N; n++) {
    std::cout << fft.t(n) << " " << C[n] << std::endl;
  }
  */

  for (unsigned n = 0; n < N; n++) {
    std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl;
  }

  /*
  // start from the exact solution for τ₀ = 0
  for (unsigned i = 0; i < N + 1; i++) {
    Real ω = i * Δω;
    Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2));
    Rt[i] = 1.0 / (μ + 1i * ω);
    Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2));
  }
  C = fft.inverse(Ct);
  R = fft.inverse(Rt);
  */

  return 0;
}