#include <chrono>
#include <fstream>
#include <iostream>

#include "animation.hpp"
#include "space_wolff.hpp"
#include "spheres.hpp"
#include "torus_symmetries.hpp"

const unsigned D = 2;
typedef Model<double, D, TorusGroup<double, D>, Dimer<double, 2>> model;

int main(int argc, char* argv[]) {
  const unsigned D = 2;

  double L = 32;
  unsigned N = 1000;
  double T = 2.0 / log(1.0 + sqrt(2.0));
  double H = 1.0;
  unsigned n = 25;
  unsigned wait = 1000;

  double k = 1e2;
  double a = 0.1;

  int opt;

  while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:w:")) != -1) {
    switch (opt) {
    case 'n':
      n = (unsigned)atof(optarg);
      break;
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 'L':
      L = atof(optarg);
      break;
    case 'T':
      T = atof(optarg);
      break;
    case 'H':
      H = atof(optarg);
      break;
    case 'a':
      a = atof(optarg);
      break;
    case 'k':
      k = atof(optarg);
      break;
    case 'w':
      wait = atoi(optarg);
      break;
    default:
      exit(1);
    }
  }

  auto zSingle = zSpheresTorus<D>(L, a, k);

  std::function<double(const Spin<double, D, Dimer<double, D>>&,
                       const Spin<double, D, Dimer<double, D>>&)>
      Z = [L, zSingle, a, k](const Spin<double, D, Dimer<double, D>>& s1,
                          const Spin<double, D, Dimer<double, D>>& s2) -> double {
    Spin<double, D, Radius> s11 = {.x = s1.x + s1.s.relativePosition, .s = s1.s.radius};
    Spin<double, D, Radius> s12 = {.x = s1.x - s1.s.relativePosition, .s = s1.s.radius};
    Spin<double, D, Radius> s21 = {.x = s2.x + s2.s.relativePosition, .s = s2.s.radius};
    Spin<double, D, Radius> s22 = {.x = s2.x - s2.s.relativePosition, .s = s2.s.radius};

    return zSingle(s11, s21) + zSingle(s12, s21) + zSingle(s11, s22) + zSingle(s12, s22);
  };

  std::function<double(Spin<double, D, Dimer<double, D>>)> B = [L, H](Spin<double, D, Dimer<double, D>> s) -> double {
    return H * s.x(1);
  };

  auto g1 = uniformGenTorus<D, Dimer<double, D>>(L);

  auto tag = std::chrono::high_resolution_clock::now();

  Animation<double, D, TorusGroup<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait, true);
  model sphere(L, Z, B);

  Rng rng;

  sphere.s.resize(n);

  unsigned nx = ceil(sqrt(n));
  for (unsigned i = 0; i < sphere.s.size(); i++) {
    Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>();
    ss->x = {(i / nx) * L / nx, (i % nx) * L / nx};
    ss->s = Dimer<double, D>{.relativePosition = {0.2, 0}, .radius = 0.25};
    sphere.s[i] = ss;
    sphere.dict.insert(ss);
  }

  sphere.wolff(T, {g1}, A, N);

  /*
  std::ofstream snapfile;
  snapfile.open("sphere_snap.dat");

  for (Spin<double, D, double>* s : sphere.s) {
    Spin<double, D, double> rs = sphere.s0.inverse().act(*s);
    snapfile << rs.s << " " << rs.x.transpose() << "\n";
    delete s;
  }

  snapfile.close();
  */

  return 0;
}