Série temporal e histograma da magnetização com campo magnético

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


// Função para atribuir as energias do sistema
double W(double, 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 na mudança de estado com campo magnético externo não-nulo
double Eflip(int l, int** matrix, int i, int j, int **vizinhoi, int **vizinhoj, double H)
{
	int k, J = 1;
	double E = 0, dE;

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

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


// Função responsável por calcular a energia total do estado com campo magnético externo não-nulo
double ETotal(int** matrix, int **vizinhoi, int **vizinhoj, int l, double H)
{
	int i, j, k, J = 1;
	double E = 0, E1 = 0, E2 = 0;

	for(i = 0; i < l; i++)
	{
		for (j = 0; j < l; j++)
		{
			for (k = 1; k < 5; k++)
			{
				E1 += matrix[i][j]*matrix[vizinhoi[i][k]][vizinhoj[j][k]];
			}

		}
		E2 += matrix[i][j];
	}
	E1 = E1*J/2.0;	//Divisão por 2 por somamos duas vezes cada interação entre vizinhos
	E2 = E2*H;
	E = -E1-E2;
	return E;
}




main()
{
	int i, j, k, l = 20, MCS = 0;	// l = tamanho do lado da matriz
 	int t = 0, h[2*l*l];
	char m[50], e[50], MH[50], R[50];
	double r, prob, M = 0, T = 2.2, ET, H = -0.005, E; 	// T = temperatura; H = Campo magnético
	FILE *arq, *arq2, *arq3;

	sprintf(m, "MagnH%lfT%lf", H, T);
	sprintf(MH, "MagnhistH%lfT%lf", H, T);
	sprintf(e, "EnerH%lfT%lf", H, T);
	sprintf(R, "ResultsH%lfT%lf", H, T);



// 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);


	// Incialização do vetor do histograma da magnetização
	for(i = 0; i < 2*l*l; i++)
		h[i] = 0;



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

	// Atribuição aletória de estratégias //

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

	srand(time(0));

	for (i = 0; i < l; i++)
	{
		for (j = 0; j < l; j++)
		{
			r = (double)rand()/RAND_MAX;
			if (r < 0.5)
				// Spins negativos recebem valor -1
				matrix[i][j] = (int)-1;
			else
				// Spins positivos recebem valor +1
				matrix[i][j] = (int)1;
		}
	}

	// Chamada para calcular a energia total do sistema
	ET = ETotal(matrix, vizinhoi, vizinhoj, l, H);
	
	// Calculo da magnetização inicial do sistema
	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 Monte Carlo Steps para o equilíbrio do sistema
	{
		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, H);
			
			// 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,T);
				r = 1.*rand()/RAND_MAX;
				if(r < prob)
				{
					matrix[i][j] = -matrix[i][j];
					ET += E;
					M += 2*matrix[i][j];
				}
			}
		}
	}
	

	arq = fopen(m, "w");
	arq2 = fopen(MH, "w");
	arq3 = fopen(e, "w");

	for (MCS = 0; MCS < 1e6; MCS++)		// Laço de Monte Carlo Steps
	{
		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, H);
			
			// 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,T);
				r = 1.*rand()/RAND_MAX;
				if(r < prob)
				{
					matrix[i][j] = -matrix[i][j];
					ET += E;
					M += 2*matrix[i][j];
				}
			}

		}
	    
		// Impressão da série temporal da magnetização
		fprintf(arq, "%d\t%lf\n", MCS, M/(l*l));
		
		// Impressão da série temporal da energia
		fprintf(arq3, "%d\t%lf\n", MCS, ET/(l*l));
		
		// Soma no histograma da magnetização
		h[(int)M+(l*l)]++;

	}


	// Impressão do histograma
	for (i = 0; i < 2*l*l; ++i)
	{
		if(i%2 == 0) // Impressão apenas dos pares pois a magnetização sempre será um número par
			fprintf(arq2, "%lf\t%lf\n", (double)(i-(l*l))/(l*l), (double)h[i]/MCS);
	}
	
	fclose(arq);
	fclose(arq2);
	fclose(arq3);


	// Impressão do estado final do sistema
	arq = fopen(R, "w");
	for (i = l-1; i >= 0; i--)
	{
		for (j = 0; j < l; j++)
		{
			fprintf(arq,"%d\t", matrix[i][j]);
		}
		fprintf(arq,"\n");
	}
	fclose(arq);

}


//Função de probabilidade de troca de estado
double W(double E, double T)
{
	return exp(-1.*E/(T));
	
}