#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 Δτ = τₘₐₓ / M_PI / n;
  Real Δω = M_PI / τₘₐₓ;

  Real y = y₀;

  FourierTransform fft(n, Δω, Δτ, FFTW_ESTIMATE);

  while (y += Δy, y <= yₘₐₓ) {
    std::vector<Real> C(2 * n);
    std::vector<Real> R(2 * n);

    std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
    if (cfile.is_open()) {
      cfile.read((char*)(C.data()), C.size() * sizeof(Real));
      cfile.close();

      std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
      rfile.read((char*)(R.data()), R.size() * 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);

      auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);

      Real z = ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();

      std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl;
    }
  }

  return 0;
}