#include <iostream>
//include for random numbers
#include <cstdlib>
//include for pow function
#include <cmath>
//include to write to file
#include <fstream>
//include to get the time for seed
#include <time.h>
 
 
using namespace std;
 
//number of time steps N, total time T = Na
const unsigned int N=400;
 
//time step a
const float a=0.25;
 
const int m=1;
const int w=1;
const int hbar=1;
const float delta=1.5;
 
float du, k, g;
float accepted = 0;
 
//generates a random number between -k and k
float randomNumber(int k) {
  return (rand()/(float)RAND_MAX)*2*k -k;
}
 
 
 
//method to write to file
int writeOutputToFile(char * filename, float x[]) {
  ofstream f;
  f.open(filename);
  for(int i=0; i<N+1; i++) {
    f << x[i] << " " << i << endl;
  }
  f.close();
  return 0;
}
 
 
//this infact computes partial S/a. Leave x a till the end to improve speed slightly?
/*
float partialS(float xi1, float xi) {
  float s;
  s = (m/2)* pow( (xi1 - xi)/a, 2 ) + 0.5 * m * pow(w*(xi1+xi)/2, 2);
  return s;
}
 
//calculate the total action
float totalS(float x[]) {
  float action = 0;
 
  //loop through all sites to get the total action by summing each partial action
  for (int i=0; i<N+1; i++) {
    action = action + partialS(x[i+1], x[i]);
  }
  
  action = a * action;
  return action;
}
*/
 
//calculate the change in the action/hbar
//x1 = u_{i-1}, x2= u_{i}, x3 = u_{i+1}
float dAction(float x1, float x2, float x3) {
  return du*(du + 2*x2 - g*(x1 + x3));
}
 
//metropolis algorithm
int metropolis(float x[]) {
  float newx, newPartialS, oldPartialS, r;
  float dS;
  //I initially had i<N+1 but this changes x[0] = x[N+1] and I need end points fixed
 
  for(int i=1; i<N; i++) {
 
    //newx is the oldx +- randomNumber between delta and -delta
    newx = x[i] + randomNumber(delta);
 
    //calculate the change in action/hbar
    dS = dAction(x[i-1], x[i], x[i+1]);
 
 
    //if the action is reduced with the new x, then set x[i] = newx;
    if (dS<0) {
      x[i] = newx;
      accepted++;
      //otherwise...
    } else {
 
      //create a random number between 0..1 with uniform probability
      r = rand()/(float)RAND_MAX;
 
 
      //check to see if we accept newx based on the random number generated and the change in action
      if (exp(-1 * dS ) > r) {
	//if true, ACCEPT!
	x[i] = newx;
	accepted++;
      } else {
	//cout << "here" << endl;
 
    }
    //cout << (accepted/(float)N) << endl;
    accepted = 0;
  }
 
  return 0;
 
}
 
float monteCarlo(float x[]) {
  float integral =0;
 
  //number of monte carlo interations
  int iter = 300;
  writeOutputToFile("output1.txt", x);
  //300 monte carlo iterations
  for(int i=0; i<iter; i++) {
 
 
 
 
    //calculate new path/configuration, a few times
    for(int j=0; j<50; j++) {
      metropolis(x);
    }
 
    //compute correlation function
    integral = integral + x[N-1] * x[1];
 
  }
 
  writeOutputToFile("output.txt", x);
  return integral/(float)iter;
 
}
 
 
//the main method!
int main() {
 
  //set up the variables a little
 
  k = 0.25 * a*a * w*w;
  g = (1-k)/(1+k);
 
  //using float because RAND_MAX is only 10^9, no point in using double
  //create an array of size N, since j=0..N-1, the Markov chain
  float u[N+1];
 
  //set the seed for rand()
  srand(time(0));
 
  //inputs random variables into sites, being between -5 and 5
  u[0] = 0;
  u[N] = 0;
  for(int i=1; i<N; i++) {
    u[i] = randomNumber(5);
  }
 
  //monteCarlo(u);
  //cout <<  monteCarlo(u) << endl;
 
  return 0;
 
}