fixed distrib_c
[physik/nlsop.git] / nlsop.c
1 /*
2  * nlsop.c 
3  *
4  * this program tries helping to understand the amorphous depuration
5  * and recrystallization of SiCx while ion implanation. hopefully the program 
6  * will simulate the stabilization of the selforganizing structure in the
7  * observed behaviour.
8  *
9  * refs: 
10  *  - J. K. N. Lindner. Habilationsschrift, Universitaet Augsburg.
11  *  - Maik Haeberlen. Diplomarbeit, Universitaet Augsburg.
12  */
13
14 #define _GNU_SOURCE
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <sys/types.h>
19 #include <sys/stat.h>
20 #include <fcntl.h>
21 #include <unistd.h>
22
23 #include "nlsop.h"
24
25 #include "dfbapi.h"
26 #include "random.h"
27
28 #define MAKE_AMORPH(N) *(N)|=AMORPH
29 #define MAKE_CRYST(N) *(N)&=~AMORPH
30
31 int usage(void)
32 {
33  puts("usage:");
34  puts("-h \t\t help");
35  puts("-n \t\t no user interaction");
36  printf("-a <value> \t slope of nuclear energy loss (default %f)\n",A_EL);
37  printf("-b <value> \t nuclear energy loss offset (default %f)\n",B_EL);
38  printf("-x <value> \t # x cells (default %d)\n",X);
39  printf("-y <value> \t # x cells (default %d)\n",Y);
40  printf("-z <value> \t # x cells (default %d)\n",Z);
41  /*
42  printf("-X <value> \t display x (default %d)\n",X/2-1);
43  printf("-Y <value> \t display y (default %d)\n",Y/2-1);
44  printf("-Z <value> \t display z (default %d)\n",Z/2-1);
45  */
46  printf("-s <value> \t steps (default %d)\n",STEPS);
47  printf("-d <value> \t refresh display (default %d)\n",REFRESH);
48  printf("-r <value> \t amorphous influence range (default %d)\n",RANGE);
49  printf("-f <value> \t pressure = <value> * 1/distance^2 (default %f)\n",A_AP);
50  printf("-p <value> \t pressure offset (default %f)\n",B_AP);
51  printf("-A <value> \t slope of linear c distribution (default %f)\n",A_CD);
52  printf("-B <value> \t linear c distribution offset (default %f)\n",B_CD);
53  printf("-C <value> \t initial c concentration (default %d)\n",CC);
54  puts("-L <file> \t load from file");
55  puts("-S <file> \t save to file");
56  puts("-R <file> \t read from random file");
57  
58  return 1;
59 }
60
61 int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,int r,double a,double b,int *t_c)
62 {
63  unsigned char *thiz;
64  int *conc;
65  int i,j;
66  double p;
67
68  thiz=d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
69  conc=d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
70  p=b*URAND_MAX;
71  for(i=-r;i<=r;i++)
72  {
73   for(j=-r;j<=r;j++)
74   {
75    if(!(i==0 && j==0))
76    {
77     if(*(d3_l->status+((x+d3_l->max_x+i)%d3_l->max_x)+((y+d3_l->max_y+j)%d3_l->max_x)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH) p+=a*URAND_MAX/(i*i+j*j);
78    } 
79   }
80  }
81  // p*=*conc;
82  // printf("debug: p * conc = %f\n",p);
83  if(!(*thiz&AMORPH))
84  {
85   if(get_rand(URAND_MAX)<=p)
86   {
87    MAKE_AMORPH(thiz);
88    printf("debug: c->a %d - %d\n",*t_c,*conc);
89    *t_c=*t_c+1-*conc;
90   } else *t_c+=1;
91  } else
92  {
93   /* assume 1-p probability */
94   if(get_rand(URAND_MAX)>p)
95   {
96    MAKE_CRYST(thiz);
97    printf("debug: a->c %d - %d\n",*t_c,*conc);
98    *t_c=*t_c+1+*conc;
99   } else *t_c+=1;
100  }
101
102  return 1;
103 }
104
105 int distrib_c(d3_lattice *d3_l,int t_c,double a,double b)
106 {
107  int i,j,k,total,area;
108  double sum;
109  int temp,left;
110  int *area_h;
111  u32 x,y,z;
112
113  area=d3_l->max_x*d3_l->max_y;
114  area_h=(int *)malloc(d3_l->max_z*sizeof(int));
115
116  total=0;
117  sum=b*d3_l->max_z+a*d3_l->max_z*(d3_l->max_z+1)/2;
118  for(i=0;i<d3_l->max_z;i++)
119  {
120   area_h[i]=0;
121   for(j=0;j<area;j++)
122   {
123    if(!(*(d3_l->status+(i*area)+j)&AMORPH))
124    {
125     area_h[i]+=1;
126    }
127   }
128   temp=(int)((i+1)*a+b)*t_c/(sum*area_h[i]);
129   // if(temp)
130   // {
131    for(j=0;j<area;j++)
132    {
133     if(!(*(d3_l->status+(i*area)+j)&AMORPH))
134     {
135      *(d3_l->extra+(i*area)+j)=temp;
136      total+=temp;
137     }
138    }
139   // }
140   left=(int)(((i+1)*a+b)*t_c/sum)%area_h[i];
141   while(left)
142   {
143    x=get_rand(d3_l->max_x);
144    y=get_rand(d3_l->max_y);
145    if(!(*(d3_l->status+(i*area)+x+y*d3_l->max_x)&AMORPH))
146    {
147     *(d3_l->extra+(i*area)+x+y*d3_l->max_x)+=1;
148     total+=1;
149     left-=1;
150    }
151   }
152  }
153  left=t_c-total;
154  while(left)
155  {
156   x=get_rand(d3_l->max_x);
157   y=get_rand(d3_l->max_y);
158   z=get_rand_lgp(d3_l->max_z,a,b);
159   if(!(*(d3_l->status+x+y*d3_l->max_x+z*area)&AMORPH))
160   {
161    *(d3_l->extra+x+y*d3_l->max_x+z*area)+=1;
162    left-=1;
163   }
164  }
165  free(area_h);
166
167  return 1;
168 }
169
170 int save_to_file(char *sf,d3_lattice *d3_l)
171 {
172  int sf_fd,c;
173
174  if((sf_fd=open(sf,O_WRONLY|O_CREAT))<0)
175  {
176   puts("cannot open save file");
177   return -1;
178  }
179  if(write(sf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
180  {
181   puts("failed saving d3 lattice struct");
182   return -1;
183  }
184  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
185  if(write(sf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
186  {
187   puts("failed saving status of d3 lattice sites");
188   return -1;
189  }
190  if(write(sf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
191  {
192   puts("failed saving sites concentration");
193   return -1;
194  }
195  close(sf_fd);
196
197  return 1;
198 }
199
200 int load_from_file(char *lf,d3_lattice *d3_l)
201 {
202  int lf_fd,c;
203
204  if((lf_fd=open(lf,O_RDONLY))<0)
205  {
206   puts("cannot open load file");
207   return -1;
208  }
209  if(read(lf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
210  {
211   puts("failed reading d3 lattice struct");
212   return -1;
213  }
214  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
215  if((d3_l->status=(unsigned char*)malloc(c*sizeof(unsigned char)))==NULL)
216  {
217   puts("cannot allocate status buffer");
218   return -1;
219  }
220  if((d3_l->extra=(int *)malloc(c*sizeof(int)))==NULL)
221  {
222   puts("cannot allocate concentration buffer");
223   return -1;
224  }
225  if(read(lf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
226  {
227   puts("failed reading status of d3 lattice sites");
228   return -1;
229  }
230  if(read(lf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
231  {
232   puts("failed reading sites concentration");
233   return -1;
234  }
235  close(lf_fd);
236
237  return 1;
238 }
239
240 int main(int argc,char **argv)
241 {
242  u32 max_x,max_y,max_z,x,y,z,x_c,y_c,z_c;
243  int i,quit,escape,nowait;
244  double a_el,b_el,a_cd,b_cd,a_ap,b_ap;
245  int cc,steps,range,refresh;
246  char s_file[MAX_CHARS];
247  char l_file[MAX_CHARS];
248  char r_file[MAX_CHARS];
249  char x_txt[MAX_TXT];
250  char y_txt[MAX_TXT];
251  char z_txt[MAX_TXT];
252  char status_txt[MAX_TXT];
253  char conc_txt[MAX_TXT];
254  char steps_txt[MAX_TXT];
255  char cc_txt[MAX_TXT];
256  char *arg_v[MAX_ARGV];
257  d3_lattice d3_l;
258
259  max_x=X;
260  max_y=Y;
261  max_z=Z;
262  steps=STEPS;
263  range=RANGE;
264  refresh=REFRESH;
265  a_el=A_EL;
266  b_el=B_EL;
267  a_cd=A_CD;
268  b_cd=B_CD;
269  a_ap=A_AP;
270  b_ap=B_AP;
271  cc=CC;
272  nowait=0;
273  quit=0;
274  escape=0;
275  strcpy(s_file,"");
276  strcpy(l_file,"");
277  strcpy(r_file,"");
278
279  for(i=1;i<argc;i++)
280  {
281   if(argv[i][0]=='-')
282   {
283    switch(argv[i][1])
284    {
285     case 'h':
286      usage();
287      return -1;
288     case 'n':
289      nowait=1;
290      break;
291     case 'a':
292      a_el=atof(argv[++i]);
293      break;
294     case 'b':
295      b_el=atof(argv[++i]);
296      break;
297     case 'x':
298      max_x=atoi(argv[++i]);
299      break;
300     case 'y':
301      max_y=atoi(argv[++i]);
302      break;
303     case 'z':
304      max_z=atoi(argv[++i]);
305      break;
306     /*
307     case 'X':
308      x=atoi(argv[++i]);
309      break;
310     case 'Y':
311      y=atoi(argv[++i]);
312      break;
313     case 'Z':
314      z=atoi(argv[++i]);
315      break;
316     */
317     case 's':
318      steps=atoi(argv[++i]);
319      break;
320     case 'd':
321      refresh=atoi(argv[++i]);
322      break;
323     case 'r':
324      range=atoi(argv[++i]);
325      break;
326     case 'f':
327      a_ap=atof(argv[++i]);
328      break;
329     case 'p':
330      b_ap=atof(argv[++i]);
331      break;
332     case 'A':
333      a_cd=atof(argv[++i]);
334      break;
335     case 'B':
336      b_cd=atof(argv[++i]);
337      break;
338     case 'C':
339      cc=atoi(argv[++i]);
340      break;
341     case 'L':
342      strcpy(l_file,argv[++i]);
343      break;
344     case 'S':
345      strcpy(s_file,argv[++i]);
346      break;
347     case 'R':
348      strcpy(r_file,argv[++i]);
349      break;
350     default:
351      usage();
352      return -1;
353    }
354   } else usage();
355  }
356
357  x=max_x/2-1;
358  y=max_y/2-1;
359  z=max_z/2-1;
360
361  if(!strcmp(r_file,"")) rand_init(NULL);
362  else rand_init(r_file);
363
364  if(!strcmp(l_file,""))
365  {
366   i=max_x*max_y*max_z;
367   d3_lattice_init(&argc,argv,&d3_l,max_x,max_y,max_z);
368   if((d3_l.status=(unsigned char *)malloc(i*sizeof(unsigned char)))==NULL)
369   {
370    puts("failed allocating status buffer");
371    return -1;
372   }
373   memset(d3_l.status,0,i*sizeof(unsigned char));
374   if((d3_l.extra=(int *)malloc(i*sizeof(int)))==NULL)
375   {
376    puts("failed allocating concentration buffer");
377    return -1;
378   }
379   memset(d3_l.extra,0,i*sizeof(int));
380  } else
381  {
382   load_from_file(l_file,&d3_l);
383   d3_lattice_init(&argc,argv,&d3_l,d3_l.max_x,d3_l.max_y,d3_l.max_z);
384  }
385
386  d3_event_init(&d3_l);
387
388  if(!strcmp(l_file,""))
389  {
390   i=0;
391   while((i<steps) && (quit==0) && (escape==0))
392   {
393    x_c=get_rand(d3_l.max_x);
394    y_c=get_rand(d3_l.max_y);
395    z_c=get_rand_lgp(d3_l.max_z,a_el,b_el);
396    distrib_c(&d3_l,cc,a_cd,b_cd);
397    process_cell(&d3_l,x_c,y_c,z_c,range,a_ap,b_ap,&cc);
398    if(i%refresh==0)
399    {
400     sprintf(x_txt,"x: %d",x+1);
401     sprintf(y_txt,"y: %d",y+1);
402     sprintf(z_txt,"z: %d",z+1);
403     sprintf(status_txt,"status: %c",(*(d3_l.status+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y)&AMORPH)?'a':'c');
404     sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
405     sprintf(steps_txt,"step: %d",i);
406     sprintf(cc_txt,"total c: %d",cc);
407     arg_v[1]=x_txt;
408     arg_v[2]=y_txt;
409     arg_v[3]=z_txt;
410     arg_v[4]=NULL;
411     arg_v[5]=status_txt;
412     arg_v[6]=conc_txt;
413     arg_v[7]=NULL;
414     arg_v[8]=NULL;
415     arg_v[9]=steps_txt;
416     arg_v[10]=cc_txt;
417     d3_lattice_draw(&d3_l,x,y,z,10,arg_v);
418     // scan_event(&d3_l,&x,&y,&z,&quit,&escape);
419    }
420    i++;
421   }
422  }
423
424  if(strcmp(s_file,"")) save_to_file(s_file,&d3_l);
425
426  while((quit==0) && (escape==0) && (nowait==0))
427  {
428   sprintf(x_txt,"x: %d",x);
429   sprintf(y_txt,"y: %d",y);
430   sprintf(z_txt,"z: %d",z);
431   sprintf(status_txt,"status: %c",(*(d3_l.status+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y)&AMORPH)?'a':'c');
432   sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
433   strcpy(steps_txt,"step: end!");
434   sprintf(cc_txt,"total c: %d",cc);
435   arg_v[1]=x_txt;
436   arg_v[2]=y_txt;
437   arg_v[3]=z_txt;
438   arg_v[4]=NULL;
439   arg_v[5]=status_txt;
440   arg_v[6]=conc_txt;
441   arg_v[7]=NULL;
442   arg_v[8]=NULL;
443   arg_v[9]=steps_txt;
444   arg_v[10]=cc_txt;
445   d3_lattice_draw(&d3_l,x,y,z,10,arg_v);
446   scan_event(&d3_l,&x,&y,&z,&quit,&escape);
447  }
448
449  return 1;
450 }