/***************************************************/
/*          Algorithm for 2-spins problem           */
/*                        J.Inoue                  */
/***************************************************/
/***************************************************/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
/*************************************************/
/*--------Definitions of Parameters--------------*/
/*************************************************/
#define N 10000    /* アンサンブルのシステム総数 */ 
#define p 0.7    /* s1の初期状態を選ぶための確率 */
#define pp 0.7   /* s2の初期状態を選ぶための確率 */ 
#define tmax 500  /* アルゴリズムをまわす回数の上限 */ 
#define T 1.01   /* ノイズレベル. アニーリングする場合にはこれをコメントアウト */
#define Z (exp((J+h1+h2)/T)+exp((-J+h1-h2)/T)+exp((-J-h1+h2)/T)+exp((J-h1-h2)/T))
#define uu0 (exp((J+h1+h2)/T)/Z) /* e1をとる定常状態確率 */
#define ud0 (exp((-J+h1-h2)/T)/Z) /* e2をとる定常状態確率 */
#define du0 (exp((-J+h2-h1)/T)/Z) /* e3をとる定常状態確率 */
#define dd0 (exp((J-h1-h2)/T)/Z)  /* e4をとる定常状態確率 */ 
#define J 1.0
#define h1 0.5
#define h2 0.1
/*****************************************************************/
/*------ "Seeds" of rundomnumbers appearing in this program -----*/
/*****************************************************************/
#define SEEDNOISE 121
#define SEEDNOISE2 142
#define SEEDNOISE3 132
#define SEEDNOISE4 45
#define SEEDNOISE5 92
/*****************************************************************/
int s1[N];
int s2[N];               
double uu;
double ud;
double du;
double dd;
//double T;  /* アニ―ルする際にはこのコメントをはずす */
/*******************************************************/
/*     Function of generating randomnumber             */
/*       by Numerical recepies in C  ran3             */
/*******************************************************/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

float ran3(long *idum)
{
	static int inext,inextp;
	static long ma[56];
	static int iff=0;
	long mj,mk;
	int i,ii,k;

	if (*idum < 0 || iff == 0) {
		iff=1;
		mj=MSEED-(*idum < 0 ? -*idum : *idum);
		mj %= MBIG;
		ma[55]=mj;
		mk=1;
		for (i=1;i<=54;i++) {
			ii=(21*i) % 55;
			ma[ii]=mk;
			mk=mj-mk;
			if (mk < MZ) mk += MBIG;
			mj=ma[ii];
		}
		for (k=1;k<=4;k++)
			for (i=1;i<=55;i++) {
				ma[i] -= ma[1+(i+30) % 55];
				if (ma[i] < MZ) ma[i] += MBIG;
			}
		inext=0;
		inextp=31;
		*idum= 1;
	}
	if (++inext == 56) inext=1;
	if (++inextp == 56) inextp=1;
	mj=ma[inext]-ma[inextp];
	if (mj < MZ) mj += MBIG;
	ma[inext]=mj;
	return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
/* (C) Copr. 1986-92 Numerical Recipes Software 1+5-5i. */

/* s1の初期値を確率的に選ぶ */ 
void set_initial()
{
        long idum=(-SEEDNOISE2);
        int x;
	for(x = 0; x < N-1;  x++){
	  if(p > ran3(&idum)){
	    s1[x] = 1;
	  }else{
	    s1[x] = -1;
	  }
	}
}

/* s2の初期値を確率的に選ぶ */ 
void set_initial2()
{
        long idum=(-SEEDNOISE3);
        int x;
	for(x = 0; x < N-1;  x++){
	  if(p > ran3(&idum)){
	    s2[x] = 1;
	  }else{
	    s2[x] = -1;
	  }
	}
}

/* 状態更新部分 */ 
void metro(int k)
{
        long idum=(-SEEDNOISE+k);
        long idum2=(-SEEDNOISE5+k);
        int x;
        double r;
           for(x = 0; x < N-1; x++){
              r = ran3(&idum2);
                if((r <= 0.5) && (ran3(&idum) 
                  < exp(-(2.0*J*s1[x]*s2[x]+2.0*h1*s1[x])/T))){
	       s1[x] = -s1[x];}else{
		 s1[x]= s1[x];}	
	      if((r > 0.5) && (ran3(&idum) 
                  < exp(-(2.0*J*s1[x]*s2[x]+2.0*h2*s2[x])/T))){
	       s2[x] = -s2[x];}else{
		 s2[x]= s2[x];}	
	   }
}
/*****************************************************/
/*                Main Program                       */
/*****************************************************/
main()
{   
	FILE *pt;
	int i,t;
	  set_initial();
	  set_initial2();
          if((pt = fopen("result.dat", "wt")) != NULL){ 
	    for(i=1,uu=0.0,ud=0.0,du=0.0,dd=0.0;i < N-1; i++){
	      if((s1[i]==1) && (s2[i]==1)){uu = uu + 1.0/N;}
	      if((s1[i]==1) && (s2[i]==-1)){ud = ud + 1.0/N;}
	      if((s1[i]==-1) && (s2[i]==1)){du = du + 1.0/N;}
	      if((s1[i]==-1) && (s2[i]==-1)){dd = dd + 1.0/N;}
	    }
	    fprintf(pt,"%d %lf %lf %lf %lf %lf %lf %lf %lf\n",0,uu,ud,du,dd,uu0,ud0,du0,dd0);
	    for(t=0; t < tmax; t=t+1){
      /* アニールする場合には下記コメントをはずす */ 
	      // T = 2.0/(log(1.0+2.0*t));
	      // T = 2.0/(1.0+2.0*t);
	      metro(t); 
	      for(i=1,uu=0.0,ud=0.0,du=0.0,dd=0.0; i < N-1; i++){
	      if((s1[i]==1) && (s2[i]==1)){uu = uu + 1.0/N;}
	      if((s1[i]==1) && (s2[i]==-1)){ud = ud + 1.0/N;}
	      if((s1[i]==-1) && (s2[i]==1)){du = du + 1.0/N;}
	      if((s1[i]==-1) && (s2[i]==-1)){dd = dd + 1.0/N;}
	      }
	       fprintf(pt,"%d %lf %lf %lf %lf %lf %lf %lf %lf\n",t,uu,ud,du,dd,uu0,ud0,du0,dd0);
	    }
	  }
             fclose(pt); 
}