ising.c

   1 /*
   2  * ising.c - visualization of ising spins in an N x N lattice
   3  *
   4  * ... actually N x M
   5  *
   6  * author: hackbard@hackdaworld.dyndns.org
   7  *
   8  */
   9
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #include <string.h>
  13 #include <math.h>
  14 #include <time.h>
  15
  16 #include "dfbapi.h"
  17
  18 #define X 200
  19 #define Y 200
  20 #define MAX_T 10
  21
  22 int main(int argc, char **argv)
  23 {
  24  unsigned char *atom;
  25  char *arg_v[4];
  26  char m_text[64];
  27  char t_text[64];
  28  char b_text[64];
  29  int max_x,x_c,max_y,y_c;
  30  int i;
  31  int count_p;
  32  unsigned int T;
  33  double S;
  34  int M;
  35  double beta;
  36  double delta_e;
  37
  38  /* random stuff*/
  39  srand(time(0));
  40
  41  /* display stuff */
  42  d2_lattice d2_l;
  43
  44  /* we will parse argv later ... */
  45  max_x=X;
  46  max_y=Y;
  47
  48  d2_lattice_init(&argc,argv,&d2_l,max_x,max_y);
  49
  50  atom=(unsigned char *)(malloc(max_x*max_y*sizeof(unsigned char)));
  51
  52  d2_l.status=atom;
  53
  54  /* begin at T=0 M=1 situation */
  55  memset(atom,0,max_x*max_y*sizeof(unsigned char));
  56
  57  S=1; /* i have no idea! */
  58
  59  for(T=1;T<MAX_T;T++)
  60  {
  61   beta=5.0/T; /* k_B = 1 */
  62   /* do 100 itterations, we will need more */
  63   for(i=0;i<100;i++)
  64   {
  65
  66  M=0;
  67  for(x_c=0;x_c<max_x;x_c++)
  68  {
  69   for(y_c=0;y_c<max_y;y_c++)
  70   {
  71    count_p=0;
  72    if((*(atom+x_c+((y_c+max_y+1)%max_y)*max_x))&1) ++count_p;
  73    if((*(atom+x_c+((y_c+max_y-1)%max_y)*max_x))&1) ++count_p;
  74    if((*(atom+((max_x+x_c+1)%max_x)+y_c*max_x))&1) ++count_p;
  75    if((*(atom+((max_x+x_c-1)%max_x)+y_c*max_x))&1) ++count_p;
  76    if(((*(atom+x_c+y_c*max_x))&1)==0) count_p=4-count_p;
  77    delta_e=(2*count_p-4)*S;
  78    if(delta_e<0) *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
  79    else
  80    {
  81     if(1.0*rand()/RAND_MAX<exp(-1.0*delta_e*beta))
  82      *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
  83    }
  84    if((*(atom+x_c+y_c*max_x))&1) ++M;
  85   }
  86  }
  87   sprintf(t_text,"temp = %d",T);
  88   arg_v[1]=t_text;
  89   sprintf(b_text,"beta = %f",beta);
  90   arg_v[2]=b_text;
  91   sprintf(m_text,"magn = %f",1.0-2.0*M/(max_x*max_y));
  92   arg_v[3]=m_text;
  93   d2_lattice_draw(&d2_l,0,0,3,arg_v);
  94   }
  95  }
  96
  97  getchar();
  98
  99  return 1;
 100 }