#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
// Parameters 
#define nspin   32
#define J       1.0
#define nmcs    1000
#define border  0
#define H_max   0.0
#define T_max   1.0
#define scale   0.4

double ransu()
{
  return rand()/(RAND_MAX+1.0);
  //return  genrand_real2();//Mersenne-Twister
}

void initialize1(int s[][nspin])
{
  int i,j;
  double r;
  for(i = 0; i < nspin ; i++){ 
    for(j = 0; j < nspin ; j++){
      r = ransu();
      //r = genrand_real2();//Mersenne-Twister
      if(r < 0.5) s[i][j] = -1;
      else        s[i][j] =  1;
    }
  }
}

double delta_E(int i,int j,int s[][nspin],double H)
{ 
  int iprev,inext,jprev,jnext;
  int topS, bottomS,leftS,rightS;
  //境界条件
  iprev = (i == 0        ? nspin-1 : i-1);
  inext = (i == nspin-1  ?  0      : i+1);
  jprev = (j == 0        ? nspin-1 : j-1);
  jnext = (j == nspin-1  ?  0      : j+1);
  topS    = s[iprev][j];
  bottomS = s[inext][j];
  leftS   = s[i][jprev];
  rightS  = s[i][jnext];
  return (2.0*J*s[i][j]*(leftS + rightS + topS + bottomS)+2.0*H*s[i][j]);
}


double total_magnetization(int s[][nspin])
{
  int i,j;
  double M;
  M = 0.0;
  for( i = 0 ; i < nspin ; i++){
    for( j = 0 ; j < nspin ; j++){
      M += s[i][j];
    }
  }
  return M;
}

double total_energy(int s[][nspin],double H)
{
  int i,j; 
  double E;
  int iprev,inext,jprev,jnext;
  int topS, bottomS,leftS,rightS;
  E = 0.0;
  for( i = 0 ; i < nspin ; i++){
    for( j = 0 ; j < nspin ; j++){
      //境界条件
      //三項演算子 
      iprev = (i == 0        ? nspin-1 : i-1); // if(i==0) {iprev=nspin-1} else  { iprev=i-1}
      inext = (i == nspin-1  ?  0      : i+1); //
      jprev = (j == 0        ? nspin-1 : j-1);
      jnext = (j == nspin-1  ?  0      : j+1);

      topS    = s[iprev][j];
      bottomS = s[inext][j];
      leftS   = s[i][jprev];
      rightS  = s[i][jnext];
      E += -1.0/2.0*J*s[i][j]*(leftS + rightS + topS + bottomS)-H*s[i][j];
    }
  }
  return E;
}

int main()
{
  FILE *fp1,*fp2;
  int s[nspin][nspin];
  int i,j,istep,mcs;
  int ispin,jspin;
  double hanten,dE;
  double H,T;
  double magnetization,magnetization_sum;
  double magnetization_2,magnetization_sum2;
  double energy,energy_sum;
  double energy_2,energy_sum2;
  double magnetic_susceptibility,specific_heat;

  initialize1(s);

  H = H_max;
  T = T_max;
  mcs = 1;
  magnetization_sum = 0.0;
  energy_sum = 0.0;
  magnetization = 0.0;
  energy = 0.0;
  magnetization_sum2 = 0.0;
  energy_sum2 = 0.0;
  magnetization_2 = 0.0;
  energy_2 = 0.0;
  magnetic_susceptibility = 0.0;
  specific_heat = 0.0;

  fp1 = fopen("step.dat","w");
  fp2 = fopen("step_spin.dat","w");
  if(fp1 == NULL || fp2 == NULL ){
    printf("failure\n");
    return 0;
  }

  fprintf(fp1,"# ising modelb step\n");
  fprintf(fp1,"# nspin=%d,J=%g,T=%g,H=%g,border=%d \n",nspin,J,T,H,border);
  fprintf(fp1,"#step   E[i]     M[i]     C[i]   x[i] \n");
  fprintf(fp2,"# ising modelb step\n");
  fprintf(fp2,"# nspin=%d,J=%g,T=%g,H=%g,border=%d \n",nspin,J,T,H,border);
  fprintf(fp2,"#   x    y     x + dx     y + dy \n");

  for( istep = 1 ; istep <= nmcs ; istep++ ){

    for( i = 0 ; i < nspin ; i++ ){
      for( j = 0 ; j < nspin ; j++ ){
	//ランダムにスピンを選ぶ
	ispin = (int)(nspin * ransu());
	jspin = (int)(nspin * ransu());
	hanten = ransu();
	dE = delta_E(ispin,jspin,s,H);
	if( dE < 0.0){
	  s[ispin][jspin] *= -1;
	}
	else{
	  if( hanten < exp(-dE/T))
	    s[ispin][jspin] *= -1;
	}
      }
    }

    if(mcs > border){
	  magnetization_sum += total_magnetization(s);// M = \sum_{i}^{nspin}M_{i}
      energy_sum += total_energy(s,H);
      magnetization = magnetization_sum/(mcs-border);//<|M|> or <M>
      energy = energy_sum/(mcs-border);//<E>
      magnetization /= (nspin*nspin);//<m>=<M>/spin
      energy /= (nspin*nspin);
      fprintf(fp1,"%5d %8.6f %8.6f \n",mcs,energy,magnetization);
    }
    mcs += 1;
    for( i = 0 ; i < nspin ; i++ ){
      for( j = 0 ; j < nspin ; j++ ){
	fprintf(fp2,"%5d %5.2f %4d %4.2f\n",i,j-scale*s[i][j]/2.0,0,s[i][j]*scale);
      }
    }
    fprintf(fp2,"\n\n");
    
  }
  
  fclose(fp1);
  fclose(fp2);
  return 0;
}