#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <math.h>
#include <TRandom3.h>
#include <TH1.h>
#include <TFile.h>
#include <TCanvas.h>

//Ising Model Program
//By Jake Connors 11 Feb 2011

using namespace std;

void Ising(){

//	TFile outfile("Isingoutput.root","RECREATE");

//	if( outfile->IsWritable()){ cout << "File is Writable!" << endl;}

///////////////////PARAMETER INPUT////////////////////////////////////

	bool debug = false; 
	
	bool corrcalc = false;//set to false if you want no correlation function calc

	static const unsigned size = 5; //length of one side

	float temperature = 4.0;//really (kT/epsilon)

	float epsilon = 1; //energy splitting between levels

	static const unsigned timesteps = 1000;

//////////////////////////////////////////////////////////////////////

	float energy = 0;

	float aveE = 0;

	int spins[size][size];

	int spinhistory[size][size][size*size*timesteps+1];

	int numiters = size*size;
	
	static const int numcorrbins = floor(size/2);

	float corr_function[numcorrbins];
	
	TH1F* MagHist = new TH1F("MagHist","Net Magnetization",100,-numiters,numiters);

	TH1F* CorrFunc = new TH1F("CorrFunc","Correlation Function",numcorrbins,0.5,numcorrbins+0.5);

	float minE = -2.0*epsilon*size*size;

	TH1F* EnergyHist = new TH1F("EnergyHist","Average Energy",100,minE,-minE);
	
	cout << "Program Initialized" << endl << endl;

	//Initialize Spins

	TRandom *r0 = new TRandom3();
	
	r0->SetSeed(0);

	int tempspin;
	
	//initializing Spins
	for(int i =0; i<size; i++){
		for(int j = 0; j<size; j++){

			float rand = r0->Rndm();

			if(rand>0.5){
				spins[i][j] = 1;
			}
			if(rand<0.5){
				spins[i][j] = -1;
			}
		
		}
	}

	system("cls");

	// To print initial spin configuration and calculate energy
	float tempE = 0;
	for(int i =0; i<size; i++){
		for(int j = 0; j<size; j++){

			spinhistory[i][j][0] = spins[i][j];
			if(spins[i][j] == 1){
				cout << "X";
			}
			if(spins[i][j] == -1){
				cout << "O";
			}

			//Calculate Energy
			int top;
			int left;
			int right;
			int bottom;
		
			if(i==0){
				top = spins[size-1][j];
			}
			if(i>0){
				top = spins[i-1][j];
			}
			if(i==(size-1)){
				bottom = spins[0][j];
			}
			if(i<(size-1)){
				bottom = spins[i+1][j];
			}
			if(j==0){
				left = spins[i][size-1];
			}
			if(j>0){
				left = spins[i][j-1];
			}
			if(j==(size-1)){
				right = spins[i][0];
			}
			if(j<(size-1)){
				right = spins[i][j+1];
			}
			tempE += (-1)*epsilon*spins[i][j]*(top+bottom+left+right);
		
		}
		cout << endl;
	}
	cout << endl << endl;
	
	energy = tempE/2.0;

	int magnetization = 0; 
	bool changed = false;

	//Looping over time from t=0 (intial config) to t = timesteps
	for(int m = 1; m < size*size*timesteps+1; m++){

		//Beginning of Ising Model Simulation

		int x = floor((r0->Rndm())*size);
		int y = floor((r0->Rndm())*size);

		if (debug){ cout << "At position :" << x << ", " << y << endl;}

		//Calculate Energy Diff
		int top;
		int left;
		int right;
		int bottom;

		if(y==0){
			top = spins[size-1][x];
		}
		if(y>0){
			top = spins[y-1][x];
		}
		if(y==(size-1)){
			bottom = spins[0][x];
		}
		if(y<(size-1)){
			bottom = spins[y+1][x];
		}
		if(x==0){
			left = spins[y][size-1];
		}
		if(x>0){
			left = spins[y][x-1];
		}
		if(x==(size-1)){
			right = spins[y][0];
		}
		if(x<(size-1)){
			right = spins[y][x+1];
		}
		float Ediff = 2*epsilon*spins[y][x]*(top+bottom+left+right);
		
		if(debug){ cout << "Ediff = " << Ediff << endl; }
		
		changed = false;

		if(Ediff<=0){
			spins[y][x] = (-1)*spins[y][x];	
			changed = true;	
		}

		if(Ediff>0){
			float chance = r0->Rndm();
			if(chance < exp((-1)*Ediff/temperature) ){
				spins[y][x] = (-1)*spins[y][x];
				changed = true;
			}
		}

		if(changed){ energy += Ediff;}

		if(m > size*size*timesteps/2.0){
			EnergyHist->Fill(energy);
			aveE += 2.0*energy/(size*size*timesteps);
		}

		//End of Ediff Calulation

		magnetization = 0;

		//Saving spin config and calculate magnetization
		for(int i =0; i<size; i++){
			for(int j = 0; j<size; j++){

				spinhistory[i][j][m] = spins[i][j];
				magnetization += spins[i][j];
			}
		}

		if(m >= size*size*timesteps/2.0){ MagHist->Fill(magnetization);}

		if(corrcalc){
		
			float tempcorr = 0.0;

			//Beginning of Correlation Function Calculation
			for(int rad = 1; rad<numcorrbins+1; rad++){
				
				tempcorr = 0.0;
				
				for(int i=0; i<size; i++){
					for(int j=0; j<size; j++){
						
						//go rad spaces left/right/up/down
						int left = j-rad;
						if(left <0){left = left+size;}
						tempcorr += 1.0*spins[i][j]*spins[i][left]/(4.0*size*size);
						int right = j+rad;
						if(right >=size){right = right-size;}
						tempcorr += 1.0*spins[i][j]*spins[i][right]/(4.0*size*size);
						int up = i-rad;
						if(up<0){up = up+size;}
						tempcorr += 1.0*spins[i][j]*spins[up][j]/(4.0*size*size);
						int down = i+rad;
						if(down>=size){down = down -size;}
						tempcorr += 1.0*spins[i][j]*spins[down][j]/(4.0*size*size);	
						
					}
				}//end loop over sites
				tempcorr = tempcorr - (1.0*magnetization*magnetization)/(1.0*size*size*size*size);
				corr_function[rad-1] = tempcorr;
				//cout << "CorrFunc[" << rad << "] = " << corr_function[rad-1] << endl;
				if(m > size*size*timesteps/2.0){
					CorrFunc->Fill(rad,corr_function[rad-1]/(0.5*size*size*timesteps));
				}
			}//end loop over radii			

		}//end if corrcalc
	
		//cout << "m is " << m << endl;

		//Printing out spin config every 100 steps
		if((floor((1.0*m)/100.0)-(1.0*m)/100.0)>-0.0001){ 
			system("cls");
			//cout << floor((1.0*m)/100.0)-(1.0*m)/100.0) << endl;
			for(int i =0; i<size; i++){
				for(int j = 0; j<size; j++){

					if(spins[i][j] == 1){
						printf("X");
					}
					if(spins[i][j] == -1){
						printf("O");
					}
				}
				cout << endl;
			}
			//gSystem->Sleep(10);
			cout << endl << m << " / " << size*size*timesteps << endl;
		}

			
	}
	
	TCanvas *c1 = new TCanvas("c1","c1");
	c1->Divide(1,3);
	c1->cd(1);
	MagHist->Draw();
	c1->cd(2);
	CorrFunc->Draw();
	c1->cd(3);
	EnergyHist->Draw();
	
/*
	outfile.cd();
	outfile.Write();
	outfile.Close();
*/
	cout << "Average Energy: " << aveE << endl;

	cout << "Program Finished" << endl << endl;
}