磁石のモデル ising

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

double minof(double e);

int main(void)
{
  FILE *MT,*Y;

  Y = fopen("spin.txt","w");
  MT = fopen("MT L20.txt","a");

#define L 100
int i,j,k,b,S[L][L],J,r;
double Ef[L][L],En[L][L],E[L][L],p,e,sum,zika,M;

J=1;
double kbT[20] = {0.2, 0.3, 0.5, 0.7, 1.0, 1.2, 1.5, 1.7, 2.0, 2.3, 2.5, 2.7, 3.0, 3.1, 3.3, 3.5, 3.8, 3.9, 4.2, 4.5};
for( b=0; b<20; b++){
sum = 0.0;
for (i=0;i<L;i++){
for (j=0;j<L;j++){

r =  rand()%100 ;

if( r % 2 == 0  ){
  S[i][j] = 1;
}
else{
  S[i][j] = -1;
}
}
}
fprintf(Y,"\nkbT=%f\n",kbT[b]);
for (k=0;k<1000;k++){
for (i=0;i<L;i++){
for (j=0;j<L;j++){

Ef[i][j] = -2*J*S[i][j] *(S[i+1][j] + S[i-1][j] + S[i][j+1] + S[i][j-1]);
En[i][j] = -Ef[i][j] ;

E[i][j] = En[i][j] - Ef[i][j];

e = exp(- E[i][j]/(kbT[b])) ;



        p = minof(e);


if ( p >= (double)rand()/RAND_MAX){
 S[i][j] = -S[i][j] ;
}
else{
  S[i][j] = S[i][j] ;
}


if (S[i][j]  == 1)
{
  fprintf (Y,"↑ ");
}
if (S[i][j] == -1)
{
  fprintf (Y,"↓ ");
}

}
fprintf(Y,"\n");
}
fprintf(Y,"\n");

}
for (k=0;k<150;k++){
for (i=0;i<L;i++){
for (j=0;j<L;j++){
sum = sum + S[i][j];

zika = sum/(L*L);
if(zika < 0){
  M = -zika;
}
else{
  M = zika;
}
}
}
}
fprintf(Y,"\nM = %f",M);
fprintf(MT,"%f %f\n",kbT[b],M);

fclose(Y);
}
fclose(MT);
return 0;
}

double minof(double e)
{
  double min;
  min = 1.0;

        if (min < e){
                min = min;}
        else{
                min = e;}
return min;

}