/* kazoeage de <E>,<E>**2,<M>,<M>**2 no keisan */

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

/* parameter */
#define J        1.0
#define T_max    10.0
#define Tbunkatu 100.0
#define dT       T_max/Tbunkatu
#define H_max    10.0
#define Hbunkatu 100
#define dH       H_max/Hbunkatu
#define kb       1.0
#define Tc       z*J/kb
#define N        100
#define nspin    2

double Z(double T)
{
  // Z=\sum_{i}^{nspin}exp(-E_{i}(S_{1},S_{2},...,S_{nspin})/(k_{B}*T))
  return (2.0*exp(1.0*J/(kb*T))+2.0*exp(-1.0*J/(kb*T)));
}

double energy(double T)
{
  // <E> = \sum_{i}^{nspin}p_{i}*E_{i}
  return ((-2.0*1.0*exp(1.0*J/(kb*T))+2.0*1.0*exp(-1.0*J/(kb*T)))/Z(T));
}

//<M>=? or <|M|>=?

void initialize1(double x[])
{
  int i;
  for(i=0;i <= N ;i++) 
    x[i] = 0.0;
}

int main()
{
  FILE *fp1;
  int i;
  double H,T;
  double E[N+1],M[N+1];

  H = 0.0;
  initialize1(E);
  initialize1(M);

	
  fp1 = fopen("k_step.dat","w");
  if(fp1 == NULL){
		printf("failure\n");
		return 0;
  }
  printf("#ising 2*2 kazoeage\n");
  fprintf(fp1,"#ising 2*2 kazoeage\n");
  //printf("#step temper mag_field energy  magnetiz spe_heat mag_susc \n");
  printf("#  i    T     H      E[i]     M[i]    \n");
  fprintf(fp1,"#  i    T     H      E[i]     M[i]    \n");
	
  for(i = Tbunkatu;i >= 0 ;i--){
    T = i * dT;
    if(i==0){
      T = 0.1;
    }
    
    E[i] = energy(T); 
	//M[i] = ?
    
    //printf("%5.2f  %8.6f  %8.6f\n",T,energy2(T),energy(T)*energy(T));
    printf("%4d %5.2f %5.2f %8.6f %8.6f \n",i,T,H,E[i],M[i]);
	fprintf(fp1,"%4d %5.2f %5.2f %8.6f %8.6f \n",i,T,H,E[i],M[i]);
  }
  printf("\n\n");

  return 0;
}