added NODFB mode
[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 implantation at temperatures
6  * below 400 degree celsius.
7  * hopefully the program will simulate the stabilization of the
8  * selforganizing lamella structure in the observed behaviour.
9  *
10  * refs: 
11  *  - J. K. N. Lindner. Habilationsschrift, Universitaet Augsburg.
12  *  - Maik Haeberlen. Diplomarbeit, Universitaet Augsburg.
13  */
14
15 #define _GNU_SOURCE
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <sys/types.h>
20 #include <sys/stat.h>
21 #include <fcntl.h>
22 #include <unistd.h>
23
24 #include "nlsop.h"
25
26 #include "dfbapi.h"
27 #include "random.h"
28
29 #define MAKE_AMORPH(N) *(N)|=AMORPH
30 #define MAKE_CRYST(N) *(N)&=~AMORPH
31
32 int usage(void)
33 {
34  puts("usage:");
35  puts("-h \t\t help");
36  puts("-n \t\t no user interaction");
37  printf("-a <value> \t slope of nuclear energy loss (default %f)\n",A_EL);
38  printf("-b <value> \t nuclear energy loss offset (default %f)\n",B_EL);
39  printf("-x <value> \t # x cells (default %d)\n",X);
40  printf("-y <value> \t # y cells (default %d)\n",Y);
41  printf("-z <value> \t # z cells (default %d)\n",Z);
42  /*
43  printf("-X <value> \t display x (default %d)\n",X/2-1);
44  printf("-Y <value> \t display y (default %d)\n",Y/2-1);
45  printf("-Z <value> \t display z (default %d)\n",Z/2-1);
46  */
47  printf("-s <value> \t steps (default %d)\n",STEPS);
48  printf("-d <value> \t refresh display (default %d)\n",REFRESH);
49  printf("-r <value> \t amorphous influence range (default %d)\n",RANGE);
50  printf("-f <value> \t pressure = <value> * 1/distance^2 (default %f)\n",A_AP);
51  printf("-p <value> \t pressure offset (default %f)\n",B_AP);
52  printf("-A <value> \t slope of linear c distribution (default %f)\n",A_CD);
53  printf("-B <value> \t linear c distribution offset (default %f)\n",B_CD);
54  /*
55  printf("-C <value> \t initial c concentration (default %d)\n",CC);
56  */
57  printf("-D <value> \t diffusion rate from cryst to amorph cells (default %f)\n",D_R);
58  printf("-W <value> \t write every <value> steps to save file (default %d)\n",RESAVE);
59  puts("-C <file> \t convert file to gnuplot format");
60  puts("-L <file> \t load from file");
61  puts("-S <file> \t save to file");
62  puts("-R <file> \t read from random file");
63  
64  return 1;
65 }
66
67 int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,int r,double a,double b,int *t_c)
68 {
69  unsigned char *thiz;
70  int *conc;
71  int i,j;
72  int off;
73  double p;
74
75  thiz=d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
76  conc=d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
77  p=b*URAND_MAX;
78  for(i=-r;i<=r;i++)
79  {
80   for(j=-r;j<=r;j++)
81   {
82    if(!(i==0 && j==0))
83    {
84     off=((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;
85     if(*(d3_l->status+off)&AMORPH) p+=a*(*(d3_l->extra+off))*URAND_MAX/(i*i+j*j);
86    } 
87   }
88  }
89  // printf("debug: p = %f\n",p);
90  if(!(*thiz&AMORPH))
91  {
92   if(get_rand(URAND_MAX)<=p) MAKE_AMORPH(thiz);
93  } else
94  {
95   /* assume 1-p probability */
96   if(get_rand(URAND_MAX)>p) MAKE_CRYST(thiz);
97  }
98  
99  *t_c+=1;
100
101  return 1;
102 }
103
104 int distrib_c(d3_lattice *d3_l,double d_r,double a,double b)
105 {
106  u32 x,y,z;
107  int i,j,k,c;
108  int offset,off;
109  int carry;
110
111  /* put one c ion somewhere in the lattice */
112  x=get_rand(d3_l->max_x);
113  y=get_rand(d3_l->max_y);
114  z=get_rand_lgp(d3_l->max_z,a,b);
115  *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1;
116
117  /* diffusion in layer */
118  for(i=0;i<d3_l->max_x;i++)
119  {
120   for(j=0;j<d3_l->max_y;j++)
121   {
122    for(k=0;k<d3_l->max_z;k++)
123    {
124     offset=i+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
125     /* case amorph */
126     if(*(d3_l->status+offset)&AMORPH)
127     {
128      /* look at neighbours and move c ions */
129      for(c=-1;c<=1;c++)
130      {
131       if(c!=0)
132       {
133        off=((i+d3_l->max_x+c)%d3_l->max_x)+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
134        carry=0;
135        /* case neighbour not amorph */
136        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off)));
137        /* case neighbour amorph */
138        /*
139         * no diffusion between amorphous cells
140         *
141        else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2;
142         */
143        if(carry!=0)
144        {
145         *(d3_l->extra+offset)+=carry;
146         *(d3_l->extra+off)-=carry;
147        }
148       }
149      }
150      for(c=-1;c<=1;c++)
151      {
152       if(c!=0)
153       {
154        off=i+((j+c+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
155        carry=0;
156        /* case neighbour not amorph */  
157        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off)));                             
158        /* case neighbour amorph */
159        /*
160         *
161         * no diffusion between amorphous cells
162         *
163        else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2;
164         */
165        if(carry!=0)
166        {
167         *(d3_l->extra+offset)+=carry; 
168         *(d3_l->extra+off)-=carry; 
169        }
170       }
171      }
172     } else
173     /* case not amorph */
174     {
175      /* look at neighbours and move c ions */     
176      for(c=-1;c<=1;c++) 
177      {
178       if(c!=0)
179       {
180        off=i+((j+c+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
181        carry=0;
182        /* case neighbour not amorph */
183        if(!(*(d3_l->status+off)&AMORPH))
184        {
185         carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2;
186         if(carry!=0)
187         {
188          *(d3_l->extra+offset)+=carry;
189          *(d3_l->extra+off)-=carry;
190         }
191        }
192       }
193      }
194     }
195    } /* for z */
196   } /* for y */
197  } /* for x */
198
199  return 1;
200 }
201
202 int distrib_c_old(d3_lattice *d3_l,int t_c,double a,double b)
203 {
204  int i,j,k,total,area;
205  double sum;
206  int temp,left;
207  int *area_h;
208  u32 x,y,z;
209
210  area=d3_l->max_x*d3_l->max_y;
211  area_h=(int *)malloc(d3_l->max_z*sizeof(int));
212
213  total=0;
214  sum=b*d3_l->max_z+a*d3_l->max_z*(d3_l->max_z+1)/2;
215  for(i=0;i<d3_l->max_z;i++)
216  {
217   area_h[i]=0;
218   for(j=0;j<area;j++)
219   {
220    if(!(*(d3_l->status+(i*area)+j)&AMORPH))
221    {
222     area_h[i]+=1;
223    }
224   }
225   temp=(int)((i+1)*a+b)*t_c/(sum*area_h[i]);
226   for(j=0;j<area;j++)
227   {
228    if(!(*(d3_l->status+(i*area)+j)&AMORPH))
229    {
230     *(d3_l->extra+(i*area)+j)=temp;
231     total+=temp;
232    }
233   }
234   left=(int)(((i+1)*a+b)*t_c/sum)%area_h[i];
235   while(left)
236   {
237    x=get_rand(d3_l->max_x);
238    y=get_rand(d3_l->max_y);
239    if(!(*(d3_l->status+(i*area)+x+y*d3_l->max_x)&AMORPH))
240    {
241     *(d3_l->extra+(i*area)+x+y*d3_l->max_x)+=1;
242     total+=1;
243     left-=1;
244    }
245   }
246  }
247  left=t_c-total;
248  while(left)
249  {
250   x=get_rand(d3_l->max_x);
251   y=get_rand(d3_l->max_y);
252   z=get_rand_lgp(d3_l->max_z,a,b);
253   if(!(*(d3_l->status+x+y*d3_l->max_x+z*area)&AMORPH))
254   {
255    *(d3_l->extra+x+y*d3_l->max_x+z*area)+=1;
256    left-=1;
257   }
258  }
259  free(area_h);
260
261  return 1;
262 }
263
264 int save_to_file(char *sf,d3_lattice *d3_l,info *my_inf)
265 {
266  int sf_fd,c;
267
268  if((sf_fd=open(sf,O_WRONLY|O_CREAT))<0)
269  {
270   puts("cannot open save file");
271   return -1;
272  }
273  if(write(sf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
274  {
275   puts("failed saving d3 lattice struct");
276   return -1;
277  }
278  if(write(sf_fd,my_inf,sizeof(info))<sizeof(info))
279  {
280   puts("failed saving info struct");
281   return-1;
282  }
283  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
284  if(write(sf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
285  {
286   puts("failed saving status of d3 lattice sites");
287   return -1;
288  }
289  if(write(sf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
290  {
291   puts("failed saving sites concentration");
292   return -1;
293  }
294  close(sf_fd);
295
296  return 1;
297 }
298
299 int load_from_file(char *lf,d3_lattice *d3_l,info *my_inf)
300 {
301  int lf_fd,c;
302
303  if((lf_fd=open(lf,O_RDONLY))<0)
304  {
305   puts("cannot open load file");
306   return -1;
307  }
308  if(read(lf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
309  {
310   puts("failed reading d3 lattice struct");
311   return -1;
312  }
313  if(read(lf_fd,my_inf,sizeof(info))<sizeof(info))
314  {
315   puts("failed reading info struct");
316   return-1;
317  }
318  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
319  if((d3_l->status=(unsigned char*)malloc(c*sizeof(unsigned char)))==NULL)
320  {
321   puts("cannot allocate status buffer");
322   return -1;
323  }
324  if((d3_l->extra=(int *)malloc(c*sizeof(int)))==NULL)
325  {
326   puts("cannot allocate concentration buffer");
327   return -1;
328  }
329  if(read(lf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
330  {
331   puts("failed reading status of d3 lattice sites");
332   return -1;
333  }
334  if(read(lf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
335  {
336   puts("failed reading sites concentration");
337   return -1;
338  }
339  close(lf_fd);
340
341  return 1;
342 }
343
344 int convert_file(char *cf,d3_lattice *d3_l)
345 {
346  int x,y,z;
347  int c_fd;
348
349  if((c_fd=open(cf,O_WRONLY|O_CREAT))<0)
350  {
351   puts("cannot open convert file");
352   return -1;
353  }
354  dprintf(c_fd,"# created by nlsop (gnuplot format)\n");
355  for(x=0;x<d3_l->max_x;x++)
356  {
357   for(y=0;y<d3_l->max_y;y++)
358   {
359    for(z=0;z<d3_l->max_z;z++)
360    {
361     if(*(d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH) dprintf(c_fd,"%d %d %d\n",x,y,z);
362    }
363   }
364  }
365  close(c_fd);
366
367  return 1;
368 }
369
370 int main(int argc,char **argv)
371 {
372  u32 x,y,z,x_c,y_c,z_c;
373  int i,quit,escape,nowait;
374  int refresh,resave;
375  char s_file[MAX_CHARS];
376  char s_file_tmp[MAX_CHARS];
377  char l_file[MAX_CHARS];
378  char c_file[MAX_CHARS];
379  char convert;
380  char r_file[MAX_CHARS];
381  char x_txt[MAX_TXT];
382  char y_txt[MAX_TXT];
383  char z_txt[MAX_TXT];
384  char status_txt[MAX_TXT];
385  char conc_txt[MAX_TXT];
386  char steps_txt[MAX_TXT];
387  char cc_txt[MAX_TXT];
388  char a_txt[MAX_TXT];
389  char s_txt[MAX_TXT];
390  char ap_txt[MAX_TXT];
391  char el_txt[MAX_TXT];
392  char cd_txt[MAX_TXT];
393  char r_txt[MAX_TXT];
394  char ap2_txt[MAX_TXT];
395  char cd2_txt[MAX_TXT];
396  char el2_txt[MAX_TXT];
397  char dr_txt[MAX_TXT];
398  char *arg_v[MAX_ARGV];
399  d3_lattice d3_l;
400  info my_info;
401
402  d3_l.max_x=X;
403  d3_l.max_y=Y;
404  d3_l.max_z=Z;
405  my_info.steps=STEPS;
406  my_info.range=RANGE;
407  refresh=REFRESH;
408  resave=RESAVE;
409  my_info.a_el=A_EL;
410  my_info.b_el=B_EL;
411  my_info.a_cd=A_CD;
412  my_info.b_cd=B_CD;
413  my_info.a_ap=A_AP;
414  my_info.b_ap=B_AP;
415  my_info.cc=CC;
416  my_info.d_r=D_R;
417  nowait=0;
418  quit=0;
419  escape=0;
420  strcpy(s_file,"");
421  strcpy(l_file,"");
422  strcpy(c_file,"");
423  convert=0;
424  strcpy(r_file,"");
425
426  for(i=1;i<argc;i++)
427  {
428   if(argv[i][0]=='-')
429   {
430    switch(argv[i][1])
431    {
432     case 'h':
433      usage();
434      return -1;
435     case 'n':
436      nowait=1;
437      break;
438     case 'a':
439      my_info.a_el=atof(argv[++i]);
440      break;
441     case 'b':
442      my_info.b_el=atof(argv[++i]);
443      break;
444     case 'x':
445      d3_l.max_x=atoi(argv[++i]);
446      break;
447     case 'y':
448      d3_l.max_y=atoi(argv[++i]);
449      break;
450     case 'z':
451      d3_l.max_z=atoi(argv[++i]);
452      break;
453     /*
454     case 'X':
455      x=atoi(argv[++i]);
456      break;
457     case 'Y':
458      y=atoi(argv[++i]);
459      break;
460     case 'Z':
461      z=atoi(argv[++i]);
462      break;
463     */
464     case 's':
465      my_info.steps=atoi(argv[++i]);
466      break;
467     case 'd':
468      refresh=atoi(argv[++i]);
469      break;
470     case 'r':
471      my_info.range=atoi(argv[++i]);
472      break;
473     case 'f':
474      my_info.a_ap=atof(argv[++i]);
475      break;
476     case 'p':
477      my_info.b_ap=atof(argv[++i]);
478      break;
479     case 'A':
480      my_info.a_cd=atof(argv[++i]);
481      break;
482     case 'B':
483      my_info.b_cd=atof(argv[++i]);
484      break;
485     /*
486     case 'C':
487      my_info.cc=atoi(argv[++i]);
488      break;
489     */
490     case 'W':
491      resave=atoi(argv[++i]);
492      break;
493     case 'C':
494      strcpy(l_file,argv[++i]);
495      if(i<argc-1) if(argv[i+1][0]!='-') strcpy(c_file,argv[++i]);
496      convert=1;
497      break;
498     case 'D':
499      my_info.d_r=atof(argv[++i]);
500      break;
501     case 'L':
502      strcpy(l_file,argv[++i]);
503      break;
504     case 'S':
505      strcpy(s_file,argv[++i]);
506      break;
507     case 'R':
508      strcpy(r_file,argv[++i]);
509      break;
510     default:
511      usage();
512      return -1;
513    }
514   } else usage();
515  }
516
517  x=d3_l.max_x/2-1;
518  y=d3_l.max_y/2-1;
519  z=d3_l.max_z/2-1;
520
521 #ifdef NODFB
522  if(!strcmp(s_file,"")
523  {
524   puts("NODFB defined, run with -S option");
525   return -1;
526  }
527
528  if(!strcmp(r_file,"")) rand_init(NULL);
529  else rand_init(r_file);
530
531  if(!strcmp(l_file,""))
532  {
533   i=d3_l.max_x*d3_l.max_y*d3_l.max_z;
534 #ifdef USE_DFB_API
535   d3_lattice_init(&argc,argv,&d3_l);
536 #endif
537   if((d3_l.status=(unsigned char *)malloc(i*sizeof(unsigned char)))==NULL)
538   {
539    puts("failed allocating status buffer");
540    return -1;
541   }
542   memset(d3_l.status,0,i*sizeof(unsigned char));
543   if((d3_l.extra=(int *)malloc(i*sizeof(int)))==NULL)
544   {
545    puts("failed allocating concentration buffer");
546    return -1;
547   }
548   memset(d3_l.extra,0,i*sizeof(int));
549  } else
550  {
551   load_from_file(l_file,&d3_l,&my_info);
552   if(convert) 
553   {   
554    if(!strcmp(c_file,"")) sprintf(c_file,"%s_gnuplot",l_file);
555    printf("converting file %s to %s\n",l_file,c_file);
556    convert_file(c_file,&d3_l);
557    puts("done");
558    return 1;
559   } 
560 #ifdef USE_DFB_API
561     else d3_lattice_init(&argc,argv,&d3_l);
562 #endif
563  }
564
565 #ifedef USE_DFB_API
566  d3_event_init(&d3_l);
567 #endif
568
569  strcpy(a_txt,"args:");
570  sprintf(s_txt,"steps: %d",my_info.steps);
571  sprintf(r_txt,"pressure range: %d",my_info.range);
572  sprintf(ap_txt,"pressure faktor: %.2f",my_info.a_ap);
573  sprintf(ap2_txt,"pressure offset: %.2f",my_info.b_ap);
574  sprintf(el_txt,"energy loss slope: %.2f",my_info.a_el);
575  sprintf(el2_txt,"energy loss offset: %.2f",my_info.b_el);
576  sprintf(cd_txt,"c distrib slope: %.2f",my_info.a_cd);
577  sprintf(cd2_txt,"c distrib offset: %.2f",my_info.b_cd);
578  sprintf(dr_txt,"diffusion rate: %.2f",my_info.d_r);
579  arg_v[1]=x_txt;
580  arg_v[2]=y_txt;
581  arg_v[3]=z_txt;
582  arg_v[4]=NULL;
583  arg_v[5]=status_txt;
584  arg_v[6]=conc_txt;
585  arg_v[7]=NULL;
586  arg_v[8]=NULL;
587  arg_v[9]=NULL;
588  arg_v[10]=steps_txt;;
589  arg_v[11]=cc_txt;
590  arg_v[12]=NULL;
591  arg_v[13]=NULL;
592  arg_v[14]=a_txt;
593  arg_v[15]=NULL;
594  arg_v[16]=s_txt;
595  arg_v[17]=r_txt;
596  arg_v[18]=ap_txt;
597  arg_v[19]=ap2_txt;
598  arg_v[20]=el_txt;
599  arg_v[21]=el2_txt;
600  arg_v[22]=cd_txt;
601  arg_v[23]=cd2_txt;
602  arg_v[24]=dr_txt;
603
604  if(!strcmp(l_file,""))
605  {
606   i=0;
607   while((i<my_info.steps) && (quit==0) && (escape==0))
608   {
609    x_c=get_rand(d3_l.max_x);
610    y_c=get_rand(d3_l.max_y);
611    z_c=get_rand_lgp(d3_l.max_z,my_info.a_el,my_info.b_el);
612    distrib_c(&d3_l,my_info.d_r,my_info.a_cd,my_info.b_cd);
613    process_cell(&d3_l,x_c,y_c,z_c,my_info.range,my_info.a_ap,my_info.b_ap,&(my_info.cc));
614    if(i%refresh==0)
615    {
616     sprintf(x_txt,"x: %d",x+1);
617     sprintf(y_txt,"y: %d",y+1);
618     sprintf(z_txt,"z: %d",z+1);
619     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');
620     sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
621     sprintf(steps_txt,"step: %d",i);
622     sprintf(cc_txt,"total c: %d",my_info.cc);
623 #ifdef USE_DFB_API
624     d3_lattice_draw(&d3_l,x,y,z,24,arg_v);
625 #endif
626     // scan_event(&d3_l,&x,&y,&z,&quit,&escape);
627    }
628    if(i%resave==0 && strcmp(s_file,"") && resave!=0)
629    {
630     sprintf(s_file_tmp,"%s_%d_of_%d",s_file,i,my_info.steps);
631     save_to_file(s_file_tmp,&d3_l,&my_info);
632    }
633    i++;
634   }
635  }
636
637  if(strcmp(s_file,"")) save_to_file(s_file,&d3_l,&my_info);
638
639  while((quit==0) && (escape==0) && (nowait==0))
640  {
641   sprintf(x_txt,"x: %d",x+1);
642   sprintf(y_txt,"y: %d",y+1);
643   sprintf(z_txt,"z: %d",z+1);
644   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');
645   sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
646   strcpy(steps_txt,"step: end!");
647   sprintf(cc_txt,"total c: %d",my_info.cc);
648 #ifdef USE_DFB_API
649   d3_lattice_draw(&d3_l,x,y,z,24,arg_v);
650   scan_event(&d3_l,&x,&y,&z,&quit,&escape);
651 #endif
652  }
653
654 #ifdef USE_DFB_API
655  d3_lattice_release(&d3_l);
656 #endif
657
658  return 1;
659 }