#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 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);

  while (y += Δy, 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;
}