Tentei dormir. Fiquei pensando sobre um problema por resolver e vim aqui "codar".
Ó:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void simula_mc(int L, int Nsteps, double temp, int **S, int **resultados, double *magnetizacao, double *energia);
void inicializa_matriz(int **S, int L);
void troca(int **S, int L, double temp);
void calcula(int **S, int L, int *mag, int *ener);
int main(int argc, char** argv){
int L, Nsteps, counter, **S, **resultados;
double temp, magnetizacao, energia;
FILE *output;
srand (1);
/* Converte os parametros de entrada em variaveis utilizaveis */
if (argc != 5){
printf("Erro: Uso: %s L Nsteps temp outfile\n", argv[0]);
exit(1);
}
L = atoi(argv[1]);
Nsteps = atoi(argv[2]);
temp = atof(argv[3]);
/* Declara a matriz de Ising S*/
S = malloc(L*sizeof(int*));
if (S == NULL){
printf("Allocation failure\n");
exit (2);
}
for( counter=0; counter < L; counter++){
S[counter] = malloc(L * sizeof(int));
if (S[counter] == NULL){
printf("Allocation failure\n");
exit (2);
}
}
/* Declara a matriz de rezultados: primeira coluna = magnetizacao, segunda coluna = energia. uma linha por passo de MC */
resultados = malloc(Nsteps*sizeof(int *));
if (resultados == NULL){
printf("Allocation failure\n");
exit (2);
}
for( counter=0; counter < Nsteps; counter++){
resultados[counter] = malloc(2 * sizeof(int));
if (resultados[counter] == NULL){
printf("Allocation failure\n");
exit (2);
}
}
simula_mc(L, Nsteps, temp, S, resultados, &magnetizacao, &energia);
output = fopen(argv[4], "w");
if (output == NULL){
printf("Erro ao abrir arquivo %s para escrita\n", argv[6]);
exit(2);
}
fprintf(output, "#Saida montecarlo: T = %g; %d passos; L = %d\n",temp, Nsteps, L);
fprintf(output, "#Magnetizacao media: %g \t Energia media %g\n", magnetizacao, energia);
fprintf(output, "#Passo\tMagnetizacao\tEnergia\n");
for( counter=0; counter < Nsteps; counter++){
fprintf(output, "%d\t%d\t%d\n", counter, resultados[counter][0], resultados[counter][1]);
}
fclose(output);
for(counter = 0; counter < L; counter++) free(S[counter]);
free (S);
for (counter = 0; counter < Nsteps; counter++) free(resultados[counter]);
free(resultados);
return 0;
}
void simula_mc(int L, int Nsteps, double temp, int **S, int **resultados, double *magnetizacao, double *energia){
int passo, passoMC, L2 = L*L;
int mag_passo, ener_passo;
temp = 1/temp;
*magnetizacao= 0;
*energia = 0;
inicializa_matriz(S, L);
/* transiente */
for (passo = 0; passo < L2*Nsteps/10; passo++)
troca(S, L, temp);
for (passoMC = 0; passoMC < Nsteps; passoMC++)
{
for(passo = 0; passo < L2; passo++)
troca(S, L, temp);
calcula(S, L, &mag_passo, &ener_passo);
*magnetizacao += abs(mag_passo);
*energia += ener_passo;
resultados[passoMC][0]= mag_passo;
resultados[passoMC][1]= ener_passo;
}
*magnetizacao = *magnetizacao/Nsteps;
*energia = *energia/Nsteps;
}
void inicializa_matriz(int **S, int L){
int i, j;
for (i=0; i < L; i++)
for (j=0; j < L; j++){
S[i][j] = 1;
}
}
void troca(int **S, int L, double temp){
int i, j, imais, imenos, jmais, jmenos;
int delta_E;
double aleat;
i = rand() % L;
j = rand() % L;
imais = (i+L+1)%L;
imenos = (i+L-1)%L;
jmais = (j+L+1)%L;
jmenos = (j+L-1)%L;
delta_E = 2 * S[i][j] * (S[i][jmais]+ S[i][jmenos] + S[imais][j] + S[imenos][j]);
if(delta_E <= 0){
S[i][j] *= -1;
}
else{
aleat = (double) rand()/RAND_MAX;
if (aleat < exp(-delta_E*temp)){
S[i][j] *= -1;
}
}
}
void calcula(int **S, int L, int *mag, int *ener){
register int i, j;
*mag=0; *ener=0;
for (i=0; i < L; i++){
for (j=0; j < L; j++){
*mag += S[i][j];
*ener += S[i][j]*(S[i][(j+1)%L] + S[(i+1)%L][j]);
}
}
}
