#include "animation.hpp"

template <int D, int n> class SoftSphere : public Sphere<D> {
private:
  const double c0 = -pow(2, 2*n-3) * pow(5, -n) * (8 + 6*n + pow(n, 2));
  const double c2 =  pow(2, 2+2*n) * pow(5, -n-2) * n * (4+n);
  const double c4 = -pow(2, 5+2*n) * pow(5, -n-4) * n * (2+n);
public:
  using Sphere<D>::Sphere;

  double U(const Position<D>& y, double ry) const {
    double r2 = Sphere<D>::regularizedDistanceSquared(y, ry);
    if (r2 < pow(1.25, 2)) {
      return pow(1 / r2, n / 2) + c0 + c2 * r2 + c4 * pow(r2, 2);
    } else {
      return 0;
    }
  }
};

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

  unsigned N = 1200;
  double R = 0.025;
  double ε = 0.005;
  unsigned steps = 1e4;
  double β = 0.5;

  initializeAnimation(argc, argv);

  Rng r;

  Model<2, SoftSphere<2, n>> m(N, 2 * R * 1.25, gContinuousPolydisperse(R), r);

  for (unsigned i = 0; i < steps; i++) {
    for (unsigned j = 0; j < N; j++) {
      m.randomMove(β, ε, r);
      m.randomSwap(β, r);
    }

//    std::cout << m.randomClusterSwap(β, r) << std::endl;
    
    if (i % 3 == 0)
      draw(m);
  }

  for (unsigned i = 0; i < steps; i++) {
    for (unsigned j = 0; j < N; j++) {
      m.randomMove(β, ε, r);
    }

    unsigned k = m.randomClusterSwap(β, r);
    if (k != 0 && k != N) {
      std::cout << k << std::endl;
    }
    
    if (i % 3 == 0)
      draw(m);
  }

  return 0;
}