/*
This program simulates time evolution of a 1D spin system with the
following properties:
1. Each spin is a Heisenberg spin (3 components).
2. There is a random magnetic field at each lattice site.
3. Only nearest neighbour interactions are taken into account and each
bond has its own coupling constant.
4. The initial orientations are random.
5. We ignore the kinetic energy.
*/

#include <iostream>
#include <fstream>
#include <math.h>
#include <ctime>
#include <boost/random.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;
using namespace boost;

// Define constants
const int length = 50; // number of lattice sites
const double energyTolerance = 0.01; // relative
const double spinLengthTolerance = 0.01; // relative
const double convergenceTolerance = 0.05; // relative
const double timeFrac = 0.1; // the relative length of the convergent region

// Define global variables
vector<double> Sinit[length]; // contains initial positions of all spins
vector<double> S[length][2]; // contains positions of all spins
vector<double> h[length]; // contains all magnetic fields
double J[length]; // contains all coupling constants
double conductivity[100000000]; // estimate of conductivity at different times
double energyCurrent; // sum of all currents until now
double deltat; // time step
double initialEnergy; // initial energy that MUST (!) be conserved throughout the simulation
bool converged; // indicates whether the calculation has converged
bool precision; // indicates whether numerical precision is ok
bool spinsConserved; // indicates whether spin lengths are conserved (check is somehow redundant but still performed)


int mod(int a, int b) // proper modulo function that does not return negative numbers!
{
  int c = a % b;

  if (c < 0)
  {
    c = c + b;
  }

  return c;
}


vector<double> planeRotation(vector<double> initVector, double angle, int axis) // perform a plane rotation about specified axis
{
  vector<double> finalVector(3);

  // apply the rotation
  finalVector(mod(axis, 3)) = initVector(mod(axis, 3));
  finalVector(mod(axis + 1, 3)) = initVector(mod(axis + 1, 3)) * cos(angle) + initVector(mod(axis + 2, 3)) * sin(angle);
  finalVector(mod(axis + 2, 3)) = initVector(mod(axis + 2, 3)) * cos(angle) - initVector(mod(axis + 1, 3)) * sin(angle);

  return finalVector;
}


vector<double> generateDirection(mt19937& gen) // generate random direction by rotating a unit vector
{
  vector<double> direction(3);
  vector<double> angles(3);

  std::tr1::uniform_real<double> rotDist(0, 2*M_PI); // uniform distribution for rotations

  // generate the angles
  for(int j = 0; j < 3; j++)
  {
    angles(j) = rotDist(gen);
  }

  // construct unit vector in the x direction
  direction(0) = 1;
  direction(1) = 0;
  direction(2) = 0;

  // apply 3 consecutive rotations
  for(int j = 0; j < 3; j++) 
  {
    direction = planeRotation(direction, angles(j), j);
  }

  return direction;
}


double calculateEnergyCurrent() // calculate the energy current
{
  double deltaEnergyCurrent = 0;

  // loop through all spins
  for(int i = 0; i < length; ++i)
  {
    if (mod(i, 2) == 0) // for even spins
    {
      deltaEnergyCurrent = deltaEnergyCurrent + prec_inner_prod( (S[mod(i, length)][1] - S[mod(i, length)][0]), (J[mod(i, length)]*S[mod(i+1, length)][0] - J[mod(i-1, length)]*S[mod(i-1, length)][0]));
    }
    else // for odd spins
    {
      deltaEnergyCurrent = deltaEnergyCurrent + prec_inner_prod( (S[mod(i, length)][1] - S[mod(i, length)][0]), (J[mod(i, length)]*S[mod(i+1, length)][1] - J[mod(i-1, length)]*S[mod(i-1, length)][1]));
    }
  }

  return deltaEnergyCurrent;
}


double calculateEnergy() // calculate the total energy of the system
{
  double energy = 0;

  // loop through all spins
  for(int i = 0; i < length; i++)
  {
    energy = energy + prec_inner_prod( S[i][1], (h[i] + J[i]*S[mod(i+1, length)][1]) );
  }

  return energy;
}


void printSpinLengths() // print length of each spin
{
  // loop through all spins
  for(int i = 0; i < length; i++)
  {
    std::cout << sqrt( prec_inner_prod( S[i][1], S[i][1] ) ) << "\n";
  }
}


void printVector(vector<double> a) // print given vector
{
  std::cout << a(0) << "," << a(1) << "," << a(2) << "\n";
}


void initialise(double J0, double epsilon, double delta) // generate random initial state
{
  // Create a Mersenne twister random number generator
  // that is seeded once with #seconds since 1970
  mt19937 gen(static_cast<unsigned> (std::time(0)));

  // uniform distribution for coupling constants
  std::tr1::uniform_real<double> coupDist(1 - epsilon, 1 + epsilon);

  // normal distribution for magnetic fields (strength)
  normal_distribution<double> hDist(0, delta);

  // bind random number generator to distribution, forming a function
  variate_generator<mt19937&, normal_distribution<double> >  hDistSampler(gen, hDist);

  for(int i = 0; i < length; i++)
  {
    J[i] = J0 * coupDist(gen); // generate coupling constants

    Sinit[i] = generateDirection(gen); // generate initial conditions

    h[i] = hDistSampler() * generateDirection(gen); // generate magnetic fields
  }
}


vector<double> getEffectiveField(int k) // calculate effective magnetic field at site k
{
  vector<double> effectiveField(3);
  effectiveField = h[mod(k, length)] + J[mod(k, length)]*S[mod(k+1, length)][1] + J[mod(k-1, length)]*S[mod(k-1, length)][1];

  return effectiveField;
}


vector<double> crossProduct(vector<double> a, vector<double> b) // calculate cross product of 2 vectors
{
  vector<double> c(3);
  c(0) = a(1)*b(2) - a(2)*b(1);
  c(1) = a(2)*b(0) - a(0)*b(2);
  c(2) = a(0)*b(1) - a(1)*b(0);
  return c;
}


void checkPrecision() // check whether total energy and spin lengths are within required tolerance
{
  // check energy conservation
  if (fabs(initialEnergy - calculateEnergy()) > (energyTolerance * fabs(initialEnergy)))
  {
    precision = false;
    std::cout << "Warning - energy is NOT conserved!\n";
    std::cout << "E0 = " << initialEnergy << ", E = " << calculateEnergy() <<"\n";
  }

  bool spinsConserved = true;

  // check the length of each spin
  for(int i = 0; i < length; i++)
  {
    if (fabs(sqrt( prec_inner_prod( S[i][1], S[i][1] ) ) - 1) > spinLengthTolerance)
    {
      precision = false;

      if (spinsConserved)
      {
        std::cout << "Warning - spin lengths are NOT conserved!\n";
        spinsConserved = false;
      }

      std::cout << "|S[" << i << "]| = " << sqrt( prec_inner_prod( S[i][1], S[i][1] ) ) << "\n";
    }
  }
}


void evolve(double dt) // evolve the system by dt
{
  vector<double> deltaS[length];  // contains the movements of all spins
  double theta;
  double phi;
  double theta2;
  double phi2;
  vector<double> H(3);

  for(int i = 0; i < length; i++)  // copy the old positions
  {
    S[i][0] = S[i][1];
  }

  for(int j = 0; j < 2; j++)  // loop through even spins first and then odd spins
  {
    for(int i = 0; i < length; i++)
    {
      if (mod(i, 2) == j)
      {
        // calculate theta from cartesian
        theta = atan(sqrt(S[i][1](0) * S[i][1](0) + S[i][1](1) * S[i][1](1)) / S[i][1](2));

        if (theta < 0)
        {
          theta = theta + M_PI;
        }

        // calculate phi from cartesian
        if (S[i][1](0) > 0)
        {
          phi = atan(S[i][1](1)/S[i][1](0));
        }
        else
        {
          phi = atan(S[i][1](1)/S[i][1](0)) + M_PI;
        }

        H = getEffectiveField(i);

        // calculate new coordinates in spherical polars
        theta2 = theta + (H(1) * cos(phi) - H(0) * sin(phi)) * deltat;
        phi2 = phi + (H(2) - (H(1) * sin(phi) + H(0) * cos(phi))/tan(theta)) * deltat;
      
        // update the cartesian coordinates
        S[i][1](0) = sin(theta2) * cos(phi2);
        S[i][1](1) = sin(theta2) * sin(phi2);
        S[i][1](2) = cos(theta2);
      }
    }
  }

  checkPrecision(); // check whether this evolution step conforms with the precision tolerance
}


bool hasConverged (int i, double frac) // check whether the conductivity has converged
{
  bool ret = true;
  int k = i;

  while (k > (1 - frac) * i)
  {
    k--;
    if (fabs(conductivity[i] - conductivity[k]) > (convergenceTolerance * fabs(conductivity[i])))
    {
      ret = false;
    }
  }

  return ret;
}


//int main (int argc, char *argv[]) // main program
int main()
{
  //deltat = strtod(argv[1], NULL);

  std::cout.precision(8); // sets precision of the standard output
  int i;

  // open output files
  std::ofstream file("points.txt");
  std::ofstream result("result.txt");

  while (true)
  {
    converged = false;
    deltat = 0.001; // set initial time step
    initialise(0.5, 0, 0.707); // generate an instance of the system (including randomness); takes J0, epsilon, delta

    while (!converged) // keep trying until converges
    {
      // set all spins to their initial conditions
      for(int j = 0; j < length; j++)
      {
        S[j][0] = Sinit[j];
        S[j][1] = Sinit[j];
      }

      initialEnergy = calculateEnergy(); // calculate initial energy

      i = 0;
      int k = 0;
      double t = 0;
      precision = true;
      energyCurrent = 0;

      while (precision & !converged) // keep evolving until either converges or loses precision
      {
        t = t + deltat;
        evolve(deltat); // evolve
        energyCurrent = energyCurrent + calculateEnergyCurrent(); // calculate new energy current
        conductivity[i] = energyCurrent*energyCurrent/(length*t); // calculate new estimate of conductivity

        if (i > 10000) // check whether the system has reached the minimum evolution time
        {
          converged = hasConverged(i, timeFrac); // check whether the system has converged
        }

        if (t > 0.1*k) // every 0.1 units of time output the conductivity
        {
          std::cout << conductivity[i] << "   time " << t << "   step " << i << "\n";
          file << t << "," << conductivity[i] << "\n";
          file.flush();
          k++;
        }

        i++;
      }

      if (!precision) // check for numerical precision
      {
        // output error and reduce time step
        std::cout << "Numerical errors. Simulation terminated.\n";
        deltat = 0.5 * deltat;
      }
      else
      {
        // save data to the output file
        result << conductivity[i - 1] << "," << t << "," << initialEnergy << "\n";
        result.flush();
      }
    }
  }

  // close both files
  file.close();
  result.close();
  return 0;
}