Suscetibilidade magnética e calor específico

#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 = Lado da matrix
 	int t = 0;
	double r, prob, M = 0, T = 2.2, ET, w[9], X, M2 = 0, M_2 = 0, C, E2, E_2; // 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("Suscetibilidade", "w");
	arq2 = fopen("Calor", "w");
	for(T = 2.1; T <= 2.6; 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];
					}
				}
			}
		}

		// Inicialização das variáveis de médias
		M2 = 0;
		M_2 = 0;
		E2 = 0;
		E_2 = 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];
					}
				}
			}

			// Somatórios de média e média quadráticas por spin para magnetização e energia
			M2 += (1.*M/(l*l))*(1.*M/(l*l));		// somatório da média quadrática de M
			E2 += (1.*ET/(l*l))*(1.*ET/(l*l));		// somatório da média quadrática de E
			E_2 += 1.*ET/(l*l);						// somatório da média de E
			M_2 += fabs(1.*M/(l*l));				// somatório da média de M
		}

		// Cálculo final das médias e médias quadráticas
		M_2 = (M_2/MCS);		// Média de M
		E_2 = (E_2/MCS);		// Média de E
		M_2 *= M_2;				// Quadrado da média de M
		E_2 *= E_2;				// Quadrado da média de E
		M2 /= MCS;				// Média quadrática de M
		E2 /= MCS;				// Média quadrática de E


		X = (l*l)*(M2-M_2)/T;			//	Calculo da suscetibilidade magnética
		C = (l*l)*(E2-E_2)/(T*T);		//	Calculo do calor específico


		// Impressão da suscetibilidade magnética e calor específico
		fprintf(arq, "%lf\t%lf\n", T, X);
		fprintf(arq2, "%lf\t%lf\n", T, C);
	}

	fclose(arq);
	fclose(arq2);
}


//Função de probabilidade de mudar de estado
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);
}