]> www.hackdaworld.org Git - physik/ising.git/blob - ising.c
0b342332c8cf0dde84c4ef50679e791f48753e2a
[physik/ising.git] / 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 }