// finalVersionCurieTemp.cpp - Alibek Medetbekov
#include <iostream>
#include <fstream>
#include <cmath>
#include <string>
#include <stdio.h>
#include <string.h>
#include <iomanip>

using namespace std;
//declaration of the size of the lattice
const unsigned int NROWS = 10, NCOLS = 10;
const unsigned int N = NROWS * NCOLS;

void save (char* pFile, int domain[][NCOLS]) //saves cofiguration of lattice electron spins
{
	int i, j;
	ofstream fout;
	fout.open(pFile, fstream::out);

	for(i = 0; i < NROWS; i++)
	{
		for(j = 0; j < NCOLS; j++)
		{
			  fout << domain[i][j];
			  if (j != NCOLS-1)
				 fout << ",";
		}          
		fout << "\n";
	}     
	fout.close();
}

double randomN() //generates a random number between 0 and 1
{
       return (double)rand()/ RAND_MAX;
}

void initialize(int domain[][NCOLS]) //initializes the array
{
	int i, j;
	for(i = 0; i < NROWS; i++)
	{
		for(j = 0; j < NCOLS; j++)
		{
			if ((double)rand()/RAND_MAX < 0.5)
				domain[i][j] = -1;
			else
				domain[i][j] = 1;   
		}
	}
}

double calcMagnetization(int domain[][NCOLS]) //calculates magnetization
{
	int i, j;
	int magnetization = 0;
	for (i = 0; i < NROWS; i++)
	{
	   for (j = 0; j < NCOLS; j++)
	   {
		   magnetization += domain[i][j];
	   }    
	}
	if (magnetization < 0) magnetization *= -1;
	return (double)magnetization / (NROWS * NCOLS);
}

double calcEnergy(int domain[][NCOLS]) // calculates energy
{
	int i, j, i_new, j_new;
    double energy = 0.0;
    for(i = 0; i < NROWS; i++)
    {
		for(j = 0; j < NCOLS; j++)
		{
			if (i == NROWS-1) i_new = 0;
			else i_new = i + 1;
			if (j == NCOLS-1) j_new = 0;
			else j_new = j + 1;

			energy += domain[i][j] * (domain[i][j_new] + domain[i_new][j]);
		}       
    }                     
    return -1.0 * energy;
}

int main()
{
	int steps = 0;
	int ddArray[NROWS][NCOLS];
	double newEnergy, oldEnergy, p, pOld, T = 2.5961, dE, totalMag = 0.0;
	int rRow, rColumn, randNum;
    double e_Gr  = -2.0 * (double)N;

    // starts streams that will save the data
	ofstream en_fout("2DEnergy.txt", fstream::out);
	ofstream mag_fout("Magnetizations.txt", fstream::out);
    ofstream fout("probability.txt", fstream::out);
    ofstream avgMag_fout("avgMag.txt", fstream::out);
	initialize(ddArray);
	save("initial2DConfiguration.txt", ddArray);
	
	
	oldEnergy = calcEnergy(ddArray);
	
	en_fout << steps << "\t" << oldEnergy << "\n";
    //random seed inorder to repeat experiements with the same random numbers
	srand(876);
    //begining of loop
	while  (steps < 60000) 
    {
	       //selects a random atom
           randNum = rand() % N;
	       rRow = randNum % NROWS;
	       rColumn = randNum / NROWS;
	       //switches the spin value at the random location to the opposite value
	       ddArray[rRow][rColumn] = -1 * ddArray[rRow][rColumn];
           //calculates energy
           newEnergy = calcEnergy(ddArray);
           //calculates the difference in energy after a spin is changed
           dE = newEnergy - oldEnergy;
        
           if(dE <= 0) // tests if the difference in energy is less than or equal to zero
           {	
    		        oldEnergy = newEnergy;      //sets new energy to old energy
    	 	        en_fout << steps << "\t" << oldEnergy << "\n"; //saves energy and magnetization
    		        mag_fout << steps << "\t" << calcMagnetization(ddArray) << "\n";
           }         
           else 
           {
                     p = exp(-dE/T); //defines probability
                     if (p > randomN())
    	             {
                  		   oldEnergy = newEnergy;
    			           steps++;       
    			           en_fout << steps << "\t" << oldEnergy << "\n";
    			           mag_fout << steps << "\t" << calcMagnetization(ddArray) << "\n";
                     }
                     else
                     {
	             	     ddArray[rRow][rColumn] = -1 * ddArray[rRow][rColumn]; 	
    			         en_fout << steps << "\t" << oldEnergy << "\n";
    			         mag_fout << steps << "\t" << calcMagnetization(ddArray) << "\n";
                     }
                     //saves probability values for trouble shooting is needed
		             fout << steps << "\t" << p << "\n"; 
		             steps++; //increments steps
           }
    }
    save("final2DConfiguration.txt", ddArray);
    avgMag_fout << totalMag/(double)steps;
	en_fout.close();
	mag_fout.close();
    avgMag_fout.close();
    fout.close();

	system("PAUSE");
	return 0;
}