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 #define _GNU_SOURCE
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #include <time.h>
  16 #include <sys/types.h>
  17 #include <sys/stat.h>
  18 #include <fcntl.h>
  19 #include <unistd.h>
  20
  21 #include "dfbapi.h"
  22
  23 #define X 200
  24 #define Y 200
  25 #define I 20
  26 #define S 1
  27
  28 int usage(void)
  29 {
  30  puts("usage:");
  31  puts("");
  32  puts("-h \t show this help");
  33  puts("-o <file> \t save T - M values to file");
  34  puts("-i <value> \t specify itteration value");
  35  puts("-x <value> \t # x lattice sites");
  36  puts("-y <value> \t # y lattice sites");
  37  puts("-s <value> \t spin interaction strength");
  38  puts("-r \t run in interactive mode (still in work)");
  39  puts("-d <value> \t refresh display rate");
  40  puts("-t <value> maximal temperature");
  41
  42  return 1;
  43 }
  44
  45 int main(int argc, char **argv)
  46 {
  47  unsigned char *atom;
  48  char *arg_v[7];
  49  char m_text[64];
  50  char t_text[64];
  51  char b_text[64];
  52  char s_text[64];
  53  char output_file[64];
  54  int of_fd=0;
  55  int max_x,x_c,max_y,y_c;
  56  int i;
  57  int itt;
  58  int count_p;
  59  double T;
  60  double s;
  61  int M=0;
  62  double beta;
  63  double delta_e;
  64  int runfree=0;
  65  double max_T;
  66  int dr=1;
  67
  68  /* random stuff*/
  69  srand(time(0));
  70
  71  /* display stuff */
  72  d2_lattice d2_l;
  73
  74  max_x=X;
  75  max_y=Y;
  76  itt=I;
  77  s=S;
  78  strcpy(output_file,"");
  79  max_T=0;
  80  /* parse argv */
  81  for(i=1;i<argc;i++)
  82  {
  83   if(argv[i][0]=='-')
  84   {
  85    switch(argv[i][1])
  86    {
  87     case 'x':
  88      max_x=atoi(argv[++i]);
  89      break;
  90     case 'y':
  91      max_y=atoi(argv[++i]);
  92      break;
  93     case 'i':
  94      itt=atoi(argv[++i]);
  95      break;
  96     case 's':
  97      s=atof(argv[++i]);
  98      break;
  99     case 'o':
 100      strcpy(output_file,argv[++i]);
 101      break;
 102     case 'r':
 103      runfree=1;
 104      break;
 105     case 'd':
 106      dr=atoi(argv[++i]);
 107      break;
 108     case 't':
 109      max_T=atof(argv[++i]);
 110      break;
 111     default:
 112      usage();
 113      return -1;
 114    }
 115   } else usage();
 116  }
 117
 118  /* prepare lattice */
 119  d2_lattice_init(&argc,argv,&d2_l,max_x,max_y);
 120  atom=(unsigned char *)(malloc(max_x*max_y*sizeof(unsigned char)));
 121  d2_l.status=atom;
 122
 123  if(strcmp(output_file,""))
 124  {
 125   of_fd=open(output_file,O_WRONLY|O_CREAT);
 126   if(of_fd<1)
 127   {
 128    puts("can't open output file ...");
 129    return -1;
 130   }
 131  }
 132
 133  /* begin at T=0 M=1 situation */
 134  memset(atom,0,max_x*max_y*sizeof(unsigned char));
 135
 136  if(max_T==0) max_T=3.0*s;
 137
 138  for(T=.05;T<max_T;T+=.05)
 139  {
 140   beta=1.0/T; /* k_B = 1 */
 141   for(i=0;i<itt;i++)
 142   {
 143
 144  M=0;
 145  for(x_c=0;x_c<max_x;x_c++)
 146  {
 147   for(y_c=0;y_c<max_y;y_c++)
 148   {
 149    count_p=0;
 150    if((*(atom+x_c+((y_c+max_y+1)%max_y)*max_x))&1) ++count_p;
 151    if((*(atom+x_c+((y_c+max_y-1)%max_y)*max_x))&1) ++count_p;
 152    if((*(atom+((max_x+x_c+1)%max_x)+y_c*max_x))&1) ++count_p;
 153    if((*(atom+((max_x+x_c-1)%max_x)+y_c*max_x))&1) ++count_p;
 154    if(((*(atom+x_c+y_c*max_x))&1)==0) count_p=4-count_p;
 155    delta_e=(2*count_p-4)*s;
 156    if(delta_e<=0) *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 157    else
 158    {
 159     if(1.0*rand()/RAND_MAX<exp(-1.0*delta_e*beta))
 160      *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 161    }
 162    if((*(atom+x_c+y_c*max_x))&1) ++M;
 163   }
 164  }
 165
 166  if(i%dr==0)
 167  {
 168   sprintf(t_text," temp: %.3f",T);
 169   arg_v[1]=t_text;
 170   sprintf(b_text," beta: %.3f",beta);
 171   arg_v[2]=b_text;
 172   arg_v[3]=NULL;
 173   sprintf(s_text," interaction strength: %.3f",s);
 174   arg_v[4]=s_text;
 175   arg_v[5]=NULL;
 176   sprintf(m_text," magnetization: %.3f",1.0-2.0*M/(max_x*max_y));
 177   arg_v[6]=m_text;
 178   d2_lattice_draw(&d2_l,0,0,6,arg_v);
 179  }
 180
 181   }
 182   if(of_fd) dprintf(of_fd,"%f %f\n",T,1.0-2.0*M/(max_x*max_y));
 183  }
 184
 185  for(T=max_T;T>0;T-=.05)
 186  {
 187   beta=1.0/T; /* k_B = 1 */
 188   for(i=0;i<itt;i++)
 189   {
 190
 191  M=0;
 192  for(x_c=0;x_c<max_x;x_c++)
 193  {
 194   for(y_c=0;y_c<max_y;y_c++)
 195   {
 196    count_p=0;
 197    if((*(atom+x_c+((y_c+max_y+1)%max_y)*max_x))&1) ++count_p;
 198    if((*(atom+x_c+((y_c+max_y-1)%max_y)*max_x))&1) ++count_p;
 199    if((*(atom+((max_x+x_c+1)%max_x)+y_c*max_x))&1) ++count_p;
 200    if((*(atom+((max_x+x_c-1)%max_x)+y_c*max_x))&1) ++count_p;
 201    if(((*(atom+x_c+y_c*max_x))&1)==0) count_p=4-count_p;
 202    delta_e=(2*count_p-4)*s;
 203    if(delta_e<=0) *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 204    else
 205    {
 206     if(1.0*rand()/RAND_MAX<exp(-1.0*delta_e*beta))
 207      *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 208    }
 209    if((*(atom+x_c+y_c*max_x))&1) ++M;
 210   }
 211  }
 212
 213  if(i%dr==0)
 214  {
 215   sprintf(t_text," temp = %.3f",T);
 216   arg_v[1]=t_text;
 217   sprintf(b_text," beta = %.3f",beta);
 218   arg_v[2]=b_text;
 219   arg_v[3]=NULL;
 220   sprintf(s_text," interaction strength: %.3f",s);
 221   arg_v[4]=s_text;
 222   arg_v[5]=NULL;
 223   sprintf(m_text," magnetization: %.3f",1.0-2.0*M/(max_x*max_y));
 224   arg_v[6]=m_text;
 225   d2_lattice_draw(&d2_l,0,0,6,arg_v);
 226  }
 227
 228   }
 229   if(of_fd) dprintf(of_fd,"%f %f\n",T,1.0-2.0*M/(max_x*max_y));
 230  }
 231
 232  if(of_fd) close(of_fd);
 233
 234  getchar();
 235
 236  return 1;
 237 }