Magnetização e energia em função da temperatura

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>


// Função para atribuir as energias do sistema
double W(double w[9], double);
///////////////////////////////////////////////////////////////////////////////////////////////////



// Função de atribuição de vizinhos em horizontais
int vizinhancaj(int l, int **vizinhoj)
{
	int i;
	for (i = 0; i < l; i++)
	{
		vizinhoj[i][1] = i+1;
		vizinhoj[i][3] = i-1;
		vizinhoj[i][2] = i;
		vizinhoj[i][4] = i;
	}
	vizinhoj[0][3] = l-1;
	vizinhoj[l-1][1] = 0;
}


// Função de atribuição de vizinhos verticais
int vizinhancai(int l, int **vizinhoi)
{
	int i;
	for (i = 0; i < l; i++)
	{
		vizinhoi[i][2] = i+1;
		vizinhoi[i][4] = i-1;
		vizinhoi[i][1] = i;
		vizinhoi[i][3] = i;
	}
	vizinhoi[0][4] = l-1;
	vizinhoi[l-1][2] = 0;
}

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



// Função responsável por calcular a variação de energia entre dois estados
double Eflip(int l, int** matrix, int i, int j, int **vizinhoi, int **vizinhoj)
{
	int k, J = 1;
	int E = 0;

	for (k = 1; k < 5; k++)
	{
		E += matrix[vizinhoi[i][k]][vizinhoj[j][k]];
	}
	E = 2*J*E*matrix[i][j];
	return E;
}

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


// Função responsável por calcular a energia total do sistema
double ETotal(int** matrix, int **vizinhoi, int **vizinhoj, int l)
{
	int i, j, k, J = 1;
	double E = 0;

	for(i = 0; i < l; i++)
		for (j = 0; j < l; j++)
		{
			for (k = 1; k < 5; k++)
			{
				E += matrix[i][j]*matrix[vizinhoi[i][k]][vizinhoj[j][k]];
			}
		}
	E *= -J;
	return E/2.0; //Divisão por 2 por somamos duas vezes cada interação entre vizinhos
}






main()
{
	int i, j, k, l = 20, MCS = 0, E;  // l = tamanho do lado da matriz
 	int t = 0;
	double r, prob, M = 0, T = 2.2, ET, w[9], X, M_t = 0, E_t = 0;	//T = temperatura
	FILE *arq, *arq2;


	

	// Incialização das matrizes utilizadas no sistema
	int ** matrix = (int **)malloc(l*sizeof(int*));
	int **vizinhoi = (int **)malloc(l*sizeof(int*));
	int **vizinhoj = (int **)malloc(l*sizeof(int*));

	for(i=0;i<l;i++)
	{
		vizinhoi[i] = (int*)malloc(5*sizeof(int));
		vizinhoj[i] = (int*)malloc(5*sizeof(int));
		matrix[i] = (int *) malloc(l*sizeof(int*));
	}

	vizinhancai(l,vizinhoi);
	vizinhancaj(l,vizinhoj);




	srand(time(0));


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

	arq = fopen("MagnXTl20", "w");
	arq2 = fopen("EnerXTl20", "w");

	for(T = 2.1; T <= 2.8; T+=0.01)		// Laço de variação da temperatura
	{
		W(w,T);		// Inicializar probabilidades de mudar de estado


		// Incialização da matriz de spins aleatóriamente para cada temperatura
		for (i = 0; i < l; i++)
		{
			for (j = 0; j < l; j++)
			{
				r = (double)rand()/RAND_MAX;
				if (r < 0.5)
					//Cooperadores recebem valor 0
					matrix[i][j] = (int)-1;
				else
					//Desertores ganham valor 1
					matrix[i][j] = (int)1;
			}
		}


		ET = ETotal(matrix, vizinhoi, vizinhoj, l);		// Calculo da energia total
		M = 0;


		// Calculo da magnetização total
		for (i = 0; i < l; i++)
			for (j = 0; j < l; j++)
				M += matrix[i][j];


		for (MCS = 0; MCS < 1e5; MCS++)			// Laço de equilíbrio do sistema
		{
			for (t = 0; t < l*l; t++)
			{
				i = rand()%l;
				j = rand()%l;
				E = Eflip(l, matrix, i, j, vizinhoi, vizinhoj);
				if(E <= 0)
				{
					matrix[i][j] = -matrix[i][j];
					ET += E;
					M += 2*matrix[i][j];
				}
				else
				{
					prob = w[E];
					r = 1.*rand()/RAND_MAX;
					if(r < prob)
					{
						matrix[i][j] = -matrix[i][j];
						ET += E;
						M += 2*matrix[i][j];
					}
				}
			}
		}


		// Incialização das variáveis de temperatura e magnetização em função da temperatura
		M_t = 0;
		E_t = 0;

		for (MCS = 0; MCS < 1e6; MCS++)		// Laço de MCS
		{
			for (t = 0; t < l*l; t++)		// Laço de tempo para 1 MCS
			{

				// Escolha aleatória de um spin
				i = rand()%l;
				j = rand()%l;

				// Calculo da variação de energia
				E = Eflip(l, matrix, i, j, vizinhoi, vizinhoj);
			
				// Condicional do método de Metropolis
				if(E <= 0)
				{
					matrix[i][j] = -matrix[i][j];
					ET += E;
					M += 2*matrix[i][j];
				}
				else
				{
					prob = w[E];
					r = 1.*rand()/RAND_MAX;
					if(r < prob)
					{
						matrix[i][j] = -matrix[i][j];
						ET += E;
						M += 2*matrix[i][j];
					}
				}
			}
			M_t += abs(M);
			E_t += ET/(l*l);
		}
	fprintf(arq, "%lf\t%lf\n", T, M_t/(l*l*MCS));
	fprintf(arq2, "%lf\t%lf\n", T, E_t/(MCS));
	}

	fclose(arq);
	fclose(arq2);
}


//Função de probabilidade de mudar para a estratégia do vizinho
double W(double w[9], double T)
{
	//return exp(-1.*E/(T));
	w[2] = exp(-2.0/T);
	w[4] = exp(-4.0/T);
	w[6] = exp(-6.0/T);
	w[8] = exp(-8.0/T);
}