!
[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 #ifdef USE_DFB_API
382  char x_txt[MAX_TXT];
383  char y_txt[MAX_TXT];
384  char z_txt[MAX_TXT];
385  char status_txt[MAX_TXT];
386  char conc_txt[MAX_TXT];
387  char steps_txt[MAX_TXT];
388  char cc_txt[MAX_TXT];
389  char a_txt[MAX_TXT];
390  char s_txt[MAX_TXT];
391  char ap_txt[MAX_TXT];
392  char el_txt[MAX_TXT];
393  char cd_txt[MAX_TXT];
394  char r_txt[MAX_TXT];
395  char ap2_txt[MAX_TXT];
396  char cd2_txt[MAX_TXT];
397  char el2_txt[MAX_TXT];
398  char dr_txt[MAX_TXT];
399  char *arg_v[MAX_ARGV];
400 #endif
401  d3_lattice d3_l;
402  info my_info;
403
404  d3_l.max_x=X;
405  d3_l.max_y=Y;
406  d3_l.max_z=Z;
407  my_info.steps=STEPS;
408  my_info.range=RANGE;
409  refresh=REFRESH;
410  resave=RESAVE;
411  my_info.a_el=A_EL;
412  my_info.b_el=B_EL;
413  my_info.a_cd=A_CD;
414  my_info.b_cd=B_CD;
415  my_info.a_ap=A_AP;
416  my_info.b_ap=B_AP;
417  my_info.cc=CC;
418  my_info.d_r=D_R;
419  nowait=0;
420  quit=0;
421  escape=0;
422  strcpy(s_file,"");
423  strcpy(l_file,"");
424  strcpy(c_file,"");
425  convert=0;
426  strcpy(r_file,"");
427
428  for(i=1;i<argc;i++)
429  {
430   if(argv[i][0]=='-')
431   {
432    switch(argv[i][1])
433    {
434     case 'h':
435      usage();
436      return -1;
437     case 'n':
438      nowait=1;
439      break;
440     case 'a':
441      my_info.a_el=atof(argv[++i]);
442      break;
443     case 'b':
444      my_info.b_el=atof(argv[++i]);
445      break;
446     case 'x':
447      d3_l.max_x=atoi(argv[++i]);
448      break;
449     case 'y':
450      d3_l.max_y=atoi(argv[++i]);
451      break;
452     case 'z':
453      d3_l.max_z=atoi(argv[++i]);
454      break;
455     /*
456     case 'X':
457      x=atoi(argv[++i]);
458      break;
459     case 'Y':
460      y=atoi(argv[++i]);
461      break;
462     case 'Z':
463      z=atoi(argv[++i]);
464      break;
465     */
466     case 's':
467      my_info.steps=atoi(argv[++i]);
468      break;
469     case 'd':
470      refresh=atoi(argv[++i]);
471      break;
472     case 'r':
473      my_info.range=atoi(argv[++i]);
474      break;
475     case 'f':
476      my_info.a_ap=atof(argv[++i]);
477      break;
478     case 'p':
479      my_info.b_ap=atof(argv[++i]);
480      break;
481     case 'A':
482      my_info.a_cd=atof(argv[++i]);
483      break;
484     case 'B':
485      my_info.b_cd=atof(argv[++i]);
486      break;
487     /*
488     case 'C':
489      my_info.cc=atoi(argv[++i]);
490      break;
491     */
492     case 'W':
493      resave=atoi(argv[++i]);
494      break;
495     case 'C':
496      strcpy(l_file,argv[++i]);
497      if(i<argc-1) if(argv[i+1][0]!='-') strcpy(c_file,argv[++i]);
498      convert=1;
499      break;
500     case 'D':
501      my_info.d_r=atof(argv[++i]);
502      break;
503     case 'L':
504      strcpy(l_file,argv[++i]);
505      break;
506     case 'S':
507      strcpy(s_file,argv[++i]);
508      break;
509     case 'R':
510      strcpy(r_file,argv[++i]);
511      break;
512     default:
513      usage();
514      return -1;
515    }
516   } else usage();
517  }
518
519  x=d3_l.max_x/2-1;
520  y=d3_l.max_y/2-1;
521  z=d3_l.max_z/2-1;
522
523 #ifdef NODFB
524  if(!strcmp(s_file,""))
525  {
526   puts("NODFB defined, run with -S option");
527   return -1;
528  }
529 #endif
530
531  if(!strcmp(r_file,"")) rand_init(NULL);
532  else rand_init(r_file);
533
534  if(!strcmp(l_file,""))
535  {
536   i=d3_l.max_x*d3_l.max_y*d3_l.max_z;
537 #ifdef USE_DFB_API
538   d3_lattice_init(&argc,argv,&d3_l);
539 #endif
540   if((d3_l.status=(unsigned char *)malloc(i*sizeof(unsigned char)))==NULL)
541   {
542    puts("failed allocating status buffer");
543    return -1;
544   }
545   memset(d3_l.status,0,i*sizeof(unsigned char));
546   if((d3_l.extra=(int *)malloc(i*sizeof(int)))==NULL)
547   {
548    puts("failed allocating concentration buffer");
549    return -1;
550   }
551   memset(d3_l.extra,0,i*sizeof(int));
552  } else
553  {
554   load_from_file(l_file,&d3_l,&my_info);
555   if(convert) 
556   {   
557    if(!strcmp(c_file,"")) sprintf(c_file,"%s_gnuplot",l_file);
558    printf("converting file %s to %s\n",l_file,c_file);
559    convert_file(c_file,&d3_l);
560    puts("done");
561    return 1;
562   } 
563 #ifdef USE_DFB_API
564     else d3_lattice_init(&argc,argv,&d3_l);
565 #endif
566  }
567
568 #ifdef USE_DFB_API
569  d3_event_init(&d3_l);
570 #endif
571
572 #ifdef USE_DFB_API
573  strcpy(a_txt,"args:");
574  sprintf(s_txt,"steps: %d",my_info.steps);
575  sprintf(r_txt,"pressure range: %d",my_info.range);
576  sprintf(ap_txt,"pressure faktor: %.2f",my_info.a_ap);
577  sprintf(ap2_txt,"pressure offset: %.2f",my_info.b_ap);
578  sprintf(el_txt,"energy loss slope: %.2f",my_info.a_el);
579  sprintf(el2_txt,"energy loss offset: %.2f",my_info.b_el);
580  sprintf(cd_txt,"c distrib slope: %.2f",my_info.a_cd);
581  sprintf(cd2_txt,"c distrib offset: %.2f",my_info.b_cd);
582  sprintf(dr_txt,"diffusion rate: %.2f",my_info.d_r);
583  arg_v[1]=x_txt;
584  arg_v[2]=y_txt;
585  arg_v[3]=z_txt;
586  arg_v[4]=NULL;
587  arg_v[5]=status_txt;
588  arg_v[6]=conc_txt;
589  arg_v[7]=NULL;
590  arg_v[8]=NULL;
591  arg_v[9]=NULL;
592  arg_v[10]=steps_txt;;
593  arg_v[11]=cc_txt;
594  arg_v[12]=NULL;
595  arg_v[13]=NULL;
596  arg_v[14]=a_txt;
597  arg_v[15]=NULL;
598  arg_v[16]=s_txt;
599  arg_v[17]=r_txt;
600  arg_v[18]=ap_txt;
601  arg_v[19]=ap2_txt;
602  arg_v[20]=el_txt;
603  arg_v[21]=el2_txt;
604  arg_v[22]=cd_txt;
605  arg_v[23]=cd2_txt;
606  arg_v[24]=dr_txt;
607 #endif
608
609  if(!strcmp(l_file,""))
610  {
611   i=0;
612   while((i<my_info.steps) && (quit==0) && (escape==0))
613   {
614    x_c=get_rand(d3_l.max_x);
615    y_c=get_rand(d3_l.max_y);
616    z_c=get_rand_lgp(d3_l.max_z,my_info.a_el,my_info.b_el);
617    distrib_c(&d3_l,my_info.d_r,my_info.a_cd,my_info.b_cd);
618    process_cell(&d3_l,x_c,y_c,z_c,my_info.range,my_info.a_ap,my_info.b_ap,&(my_info.cc));
619 #ifdef USE_DFB_API
620    if(i%refresh==0)
621    {
622     sprintf(x_txt,"x: %d",x+1);
623     sprintf(y_txt,"y: %d",y+1);
624     sprintf(z_txt,"z: %d",z+1);
625     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');
626     sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
627     sprintf(steps_txt,"step: %d",i);
628     sprintf(cc_txt,"total c: %d",my_info.cc);
629     d3_lattice_draw(&d3_l,x,y,z,24,arg_v);
630    }
631 #endif
632    if(i%resave==0 && strcmp(s_file,"") && resave!=0)
633    {
634     sprintf(s_file_tmp,"%s_%d_of_%d",s_file,i,my_info.steps);
635     save_to_file(s_file_tmp,&d3_l,&my_info);
636 #ifdef NODFB
637     printf("saved %s\n",s_file_tmp);
638 #endif
639    }
640    i++;
641   }
642  }
643
644  if(strcmp(s_file,""))
645  {
646    printf("saved %s\n",s_file);
647    save_to_file(s_file,&d3_l,&my_info);
648  }
649
650 #ifdef USE_DFB_API
651  while((quit==0) && (escape==0) && (nowait==0))
652  {
653   sprintf(x_txt,"x: %d",x+1);
654   sprintf(y_txt,"y: %d",y+1);
655   sprintf(z_txt,"z: %d",z+1);
656   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');
657   sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
658   strcpy(steps_txt,"step: end!");
659   sprintf(cc_txt,"total c: %d",my_info.cc);
660   d3_lattice_draw(&d3_l,x,y,z,24,arg_v);
661   scan_event(&d3_l,&x,&y,&z,&quit,&escape);
662  }
663
664  d3_lattice_release(&d3_l);
665 #endif
666
667  return 1;
668 }