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
  40  return 1;
  41 }
  42
  43 int main(int argc, char **argv)
  44 {
  45  unsigned char *atom;
  46  char *arg_v[7];
  47  char m_text[64];
  48  char t_text[64];
  49  char b_text[64];
  50  char s_text[64];
  51  char output_file[64];
  52  int of_fd=0;
  53  int max_x,x_c,max_y,y_c;
  54  int i;
  55  int itt;
  56  int count_p;
  57  double T;
  58  double s;
  59  int M=0;
  60  double beta;
  61  double delta_e;
  62  int runfree=0;
  63  double max_T;
  64
  65  /* random stuff*/
  66  srand(time(0));
  67
  68  /* display stuff */
  69  d2_lattice d2_l;
  70
  71  max_x=X;
  72  max_y=Y;
  73  itt=I;
  74  s=S;
  75  strcpy(output_file,"");
  76  /* parse argv */
  77  for(i=1;i<argc;i++)
  78  {
  79   if(argv[i][0]=='-')
  80   {
  81    switch(argv[i][1])
  82    {
  83     case 'x':
  84      max_x=atoi(argv[++i]);
  85      break;
  86     case 'y':
  87      max_y=atoi(argv[++i]);
  88      break;
  89     case 'i':
  90      itt=atoi(argv[++i]);
  91      break;
  92     case 's':
  93      s=atof(argv[++i]);
  94      break;
  95     case 'o':
  96      strcpy(output_file,argv[++i]);
  97      break;
  98     case 'r':
  99      runfree=1;
 100      break;
 101     default:
 102      usage();
 103      return -1;
 104    }
 105   } else usage();
 106  }
 107
 108  /* prepare lattice */
 109  d2_lattice_init(&argc,argv,&d2_l,max_x,max_y);
 110  atom=(unsigned char *)(malloc(max_x*max_y*sizeof(unsigned char)));
 111  d2_l.status=atom;
 112
 113  if(strcmp(output_file,""))
 114  {
 115   of_fd=open(output_file,O_WRONLY|O_CREAT);
 116   if(of_fd<1)
 117   {
 118    puts("can't open output file ...");
 119    return -1;
 120   }
 121  }
 122
 123  /* begin at T=0 M=1 situation */
 124  memset(atom,0,max_x*max_y*sizeof(unsigned char));
 125
 126  max_T=1.3*s;
 127
 128  for(T=.05;T<max_T;T+=.05)
 129  {
 130   beta=1.0/T; /* k_B = 1 */
 131   for(i=0;i<itt;i++)
 132   {
 133
 134  M=0;
 135  for(x_c=0;x_c<max_x;x_c++)
 136  {
 137   for(y_c=0;y_c<max_y;y_c++)
 138   {
 139    count_p=0;
 140    if((*(atom+x_c+((y_c+max_y+1)%max_y)*max_x))&1) ++count_p;
 141    if((*(atom+x_c+((y_c+max_y-1)%max_y)*max_x))&1) ++count_p;
 142    if((*(atom+((max_x+x_c+1)%max_x)+y_c*max_x))&1) ++count_p;
 143    if((*(atom+((max_x+x_c-1)%max_x)+y_c*max_x))&1) ++count_p;
 144    if(((*(atom+x_c+y_c*max_x))&1)==0) count_p=4-count_p;
 145    delta_e=(2*count_p-4)*s;
 146    if(delta_e<0) *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 147    else
 148    {
 149     if(1.0*rand()/RAND_MAX<exp(-1.0*delta_e*beta))
 150      *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 151    }
 152    if((*(atom+x_c+y_c*max_x))&1) ++M;
 153   }
 154  }
 155
 156   sprintf(t_text," temp: %.3f",T);
 157   arg_v[1]=t_text;
 158   sprintf(b_text," beta: %.3f",beta);
 159   arg_v[2]=b_text;
 160   arg_v[3]=NULL;
 161   sprintf(s_text," interaction strength: %.3f",s);
 162   arg_v[4]=s_text;
 163   arg_v[5]=NULL;
 164   sprintf(m_text," magnetization: %.3f",1.0-2.0*M/(max_x*max_y));
 165   arg_v[6]=m_text;
 166   d2_lattice_draw(&d2_l,0,0,6,arg_v);
 167
 168   }
 169   if(of_fd) dprintf(of_fd,"%f %f\n",T,1.0-2.0*M/(max_x*max_y));
 170  }
 171
 172  for(T=max_T;T>0;T-=.05)
 173  {
 174   beta=1.0/T; /* k_B = 1 */
 175   for(i=0;i<itt;i++)
 176   {
 177
 178  M=0;
 179  for(x_c=0;x_c<max_x;x_c++)
 180  {
 181   for(y_c=0;y_c<max_y;y_c++)
 182   {
 183    count_p=0;
 184    if((*(atom+x_c+((y_c+max_y+1)%max_y)*max_x))&1) ++count_p;
 185    if((*(atom+x_c+((y_c+max_y-1)%max_y)*max_x))&1) ++count_p;
 186    if((*(atom+((max_x+x_c+1)%max_x)+y_c*max_x))&1) ++count_p;
 187    if((*(atom+((max_x+x_c-1)%max_x)+y_c*max_x))&1) ++count_p;
 188    if(((*(atom+x_c+y_c*max_x))&1)==0) count_p=4-count_p;
 189    delta_e=(2*count_p-4)*s;
 190    if(delta_e<0) *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 191    else
 192    {
 193     if(1.0*rand()/RAND_MAX<exp(-1.0*delta_e*beta))
 194      *(atom+x_c+y_c*max_x)=(*(atom+x_c+y_c*max_x)+1)&1;
 195    }
 196    if((*(atom+x_c+y_c*max_x))&1) ++M;
 197   }
 198  }
 199
 200   sprintf(t_text," temp = %.3f",T);
 201   arg_v[1]=t_text;
 202   sprintf(b_text," beta = %.3f",beta);
 203   arg_v[2]=b_text;
 204   arg_v[3]=NULL;
 205   sprintf(s_text," interaction strength: %.3f",s);
 206   arg_v[4]=s_text;
 207   arg_v[5]=NULL;
 208   sprintf(m_text," magnetization: %.3f",1.0-2.0*M/(max_x*max_y));
 209   arg_v[6]=m_text;
 210   d2_lattice_draw(&d2_l,0,0,6,arg_v);
 211
 212   }
 213   if(of_fd) dprintf(of_fd,"%f %f\n",T,1.0-2.0*M/(max_x*max_y));
 214  }
 215
 216  if(of_fd) close(of_fd);
 217
 218  getchar();
 219
 220  return 1;
 221 }