//2D Ising Model
//Mayrolin Garcia ~ August 27, 2009


#include <iostream>
#include <cstdlib>
#include <fstream>
#include <cmath>
#include <stdio.h>
using namespace std;

#define NX 14
#define NY 14
#define N (NX*NY)

/////////////////////////////////////////////////////////////////////
// Function to Calculate Energy
//
double EnergyCalc(int spins[][NX])
{     
	double energy = 0.0;
	int ir, ic, row, col;
	for (row = 0; row < NY; row++) {
		for (col = 0; col < NX; col++){
			// periodic... boundary...
			if (row == NY-1) ir = 0;
			else ir = row + 1;
			if (col == NX-1) ic = 0;
			else ic = col + 1;

			energy = energy + spins[row][col] * (spins[row][ic] + spins[ir][col]);
		}
	}
	return -1.0 * energy;
}

/////////////////////////////////////////////////////////////////////
// Function to Calculate Magnetization
//
double Magnetization(int spins[][NX]) 
{
	int MG = 0;
	int row, col;
	for (row = 0; row < NY; row++){
		for (col=0; col < NX; col++){
			MG += spins[row][col];
		}
	}
	if (MG < 0) MG = MG * -1;
	return (double)MG / N;
}
               
/////////////////////////////////////////////////////////////////////
// Function to Export Spin-Map
//
void SaveSpins(int spins[][NX], char filename[])
{
	ofstream fout(filename, fstream::out);
	int row, col;
	for (row = 0; row < NY; row++){
		for (col = 0; col < NX; col++){
			fout << spins[row][col];
			if (col < NX-1) fout << ",";
		}
		fout << "\n";
	}
	fout.close();
}
/////////////////////////////////////////////////////////////////////
// Function to Get Spins
//
int GetRandomSpin()
{ 
	if ( ((double)rand() / RAND_MAX) < 0.5)
		return 1;
	else
		return  -1;
}

/////////////////////////////////////////////////////////////////////
// Function to Initialize Array
//
void Initialize(int spins [][NX])
{
	int row, col;
	for (row = 0; row < NY; row++) {
		for (col = 0; col < NX; col++) {
		   spins[row][col] = GetRandomSpin();
		}
	}
}


////////////////////////////////////////////////////////////////////
// Random Function Generator
//

double UnitRandom ()
{
    return (double)rand()/RAND_MAX;           
}


/////////////////////////////////////////////////////////////////////
// Code entry point
//
int main()
{
	double en_old, en_new, T;
    int p,q;
	int steps = 0;
	int spins[NY][NX];
	//energy gound state
    int en_grst = -2 * N;
    int counter = 0;
    
    srand (77);

	Initialize(spins);
	SaveSpins(spins, "init-config .8.txt");
	
	ofstream fout("ground-state .8.txt", fstream::out);
	ofstream fout2("magentization .8.txt", fstream::out);

    en_old = EnergyCalc(spins);   
   
    while (steps < 50000)//en_old>en_grst && 
	{
		// pick random spin site
        p = rand() % NY; 
        q = rand() % NX;
        spins[p][q] = -1 * spins[p][q];

        // calculate the new energy of the new configuration
        en_new = EnergyCalc(spins);
        
        // calculate delta Energy
        double dEnergy = en_new-en_old;
        
        T = .8;
        
       if(dEnergy <= 0)
       {
                       en_old = en_new;
                       fout << steps << "\t" << en_old << "\n";
                       fout2 << steps << "\t" << Magnetization(spins) << "\n";
                       //totalMag += calcMagnetization(ddArray);
               }

               else
       {
           double prob = exp(-dEnergy/T);
            if (prob > UnitRandom() )//(newEnergy <= oldEnergy)
                {
                  en_old = en_new;
                          fout << steps << "\t" << en_old << "\n";
                          fout2 << steps << "\t" << Magnetization(spins) << "\n";
                          //totalMag += calcMagnetization(ddArray);
            }
            else
                            spins[p][q] = -1 * spins[p][q];
                            fout << steps << "\t" << en_old << "\n";
                            fout2 << steps << "\t" << Magnetization(spins) << "\n";
               }
		steps++;
    } 


    SaveSpins(spins, "final-config .8.txt");

    fout.close();
    fout2.close();

	// pausing the screen, comment this out when don't need it
	system("PAUSE");
	return 0;
}