#include "animation.hpp"
#include <fstream>
#include "ising.hpp"

int main(int argc, char* argv[]) {
  unsigned L = 32;
  unsigned N = 1000;
  unsigned mod = 0;
  double mag = 0.5;
  double pop = 1.0;
  double T = 2.0 / log(1.0 + sqrt(2.0));
  double H = 1.0;
  bool color = false;
  unsigned wait = N;

  int opt;

  while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:cw:")) != -1) {
    switch (opt) {
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 'L':
      L = atoi(optarg);
      break;
    case 'T':
      T = atof(optarg);
      break;
    case 'H':
      H = atof(optarg);
      break;
    case 'm':
      mod = atoi(optarg);
      break;
    case 'r':
      mag = atof(optarg);
      break;
    case 'p':
      pop = atof(optarg);
      break;
    case 'c':
      color = true;
      break;
    case 'w':
      wait = atoi(optarg);
      break;
    default:
      exit(1);
    }
  }

  std::function<double(Spin<signed, D, signed>)> B;

  if (mod == 0) {
    B = isingBFace(L, H);
  } else {
    B = isingBMod(L, mod, H);
  }

  isingModel ising(L, isingZ(L), B);
  isingPopulate(ising, L, pop, mag);

  auto g = isingGen(L);

  Animation<signed, D, TorusGroup<signed, D>, IsingSpin> A(L, 750, argc, argv, wait, true);

  ising.wolff(T, {g}, A, N);

  std::vector<std::vector<int>> finalState(L);

  for (std::vector<int>& v : finalState) {
    v.resize(L);
  }

  for (const isingSpin* s: ising.s) {
    Vector<int, 2> v = ising.s0.inverse().act(s->x);
    finalState[v[0]][v[1]] = s->s;
  }

  std::ofstream outfile("ising_snap.dat");

  for (const std::vector<int>& v : finalState) {
    for (int s : v) {
      outfile << s << " ";
    }
    outfile << "\n";
  }

  outfile.close();

  return 0;
}