added test: no c diff (-i option)
[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 /* test globals - get removed or included in my_info struct later */
33 int amorph_count;
34 int cryst_count;
35 unsigned char c_diff;
36
37 int usage(void)
38 {
39  puts("usage:");
40  puts("-h \t\t help");
41  puts("-n \t\t no user interaction");
42  puts("-Z \t\t cryst -> amorph c diffusion in z direction");
43  puts("-i \t\t no cryst to cryst diffusion");
44  printf("-a <value> \t slope of nuclear energy loss (default %f)\n",A_EL);
45  printf("-b <value> \t nuclear energy loss offset (default %f)\n",B_EL);
46  printf("-x <value> \t # x cells (default %d)\n",X);
47  printf("-y <value> \t # y cells (default %d)\n",Y);
48  printf("-z <value> \t # z cells (default %d)\n",Z);
49  printf("-s <value> \t steps (default %d)\n",STEPS);
50  printf("-d <value> \t refresh display (default %d)\n",REFRESH);
51  printf("-r <value> \t amorphous influence range (default %d)\n",RANGE);
52  printf("-f <value> \t pressure = <value> * 1/distance^2 (default %f)\n",A_AP);
53  printf("-p <value> \t pressure offset (default %f)\n",B_AP);
54  printf("-F <value> \t proportionality constant between c conc and ability to get amorphous (default %f)\n",A_CP);
55  printf("-A <value> \t slope of linear c distribution (default %f)\n",A_CD);
56  printf("-B <value> \t linear c distribution offset (default %f)\n",B_CD);
57  printf("-D <value> \t diffusion rate from cryst to amorph cells (default %f)\n",DR_AC);
58  printf("-c <value> \t diffusion rate in cryst cells (default %f)\n",DR_CC);
59  printf("-e <value> \t do diffusion every <value> steps (default %d)\n",DIFF_RATE);
60  puts("-g <file> <step> continue simulation from file and step (step > 0)!");
61  printf("-W <value> \t write every <value> steps to save file (default %d)\n",RESAVE);
62  puts("-C <file> \t convert file to gnuplot format");
63  puts("-L <file> \t load from file");
64  puts("-S <file> \t save to file");
65  puts("-R <file> \t read from random file");
66  puts("-P <file> \t specify implantatin profile file");
67  
68  return 1;
69 }
70
71 int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info)
72 {
73  unsigned char *thiz;
74  int *conc;
75  int i,j;
76  int off;
77  double p;
78
79  thiz=d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
80  conc=d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
81  p=my_info->b_ap*URAND_MAX;
82  for(i=-(my_info->range);i<=my_info->range;i++)
83  {
84   for(j=-(my_info->range);j<=my_info->range;j++)
85   {
86    if(!(i==0 && j==0))
87    {
88     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;
89     if(*(d3_l->status+off)&AMORPH) p+=my_info->a_ap*(*(d3_l->extra+off))*URAND_MAX/(i*i+j*j);
90    } 
91   }
92  }
93  p+=*conc*my_info->a_cp*URAND_MAX;
94  if(!(*thiz&AMORPH))
95  {
96   if(get_rand(URAND_MAX)<=p)
97   {
98    MAKE_AMORPH(thiz);
99    amorph_count++;
100   }
101  } else
102  {
103   /* assume 1-p probability */
104   if(get_rand(URAND_MAX)>p)
105   {
106    MAKE_CRYST(thiz);
107    cryst_count++;
108   }
109  }
110  
111  return 1;
112 }
113
114 int distrib_c(d3_lattice *d3_l,info *my_info,int step,double c_ratio)
115 {
116  u32 x,y,z;
117  int i,j,k,c;
118  int offset,off;
119  int carry;
120
121  /* put one c ion somewhere in the lattice */
122  if(my_info->cc<c_ratio*step)
123  {
124   x=get_rand(d3_l->max_x);
125   y=get_rand(d3_l->max_y);
126   z=get_rand_lgp(d3_l->max_z,my_info->a_cd,my_info->b_cd);
127   *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1;
128   (my_info->cc)++;
129  }
130
131  if(step%my_info->diff_rate==0)
132  {
133
134  for(i=0;i<d3_l->max_x;i++)
135  {
136   for(j=0;j<d3_l->max_y;j++)
137   {
138    for(k=0;k<d3_l->max_z;k++)
139    {
140     offset=i+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
141     /* case amorph: amorph <- cryst diffusion */
142     if(*(d3_l->status+offset)&AMORPH)
143     {
144      for(c=-1;c<=1;c++)
145      {
146       if(c!=0)
147       {
148        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;
149        carry=0;
150        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
151        if(carry!=0)
152        {
153         *(d3_l->extra+offset)+=carry;
154         *(d3_l->extra+off)-=carry;
155        }
156       }
157      }
158      for(c=-1;c<=1;c++)
159      {
160       if(c!=0)
161       {
162        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;
163        carry=0;
164        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));                             
165        if(carry!=0)
166        {
167         *(d3_l->extra+offset)+=carry; 
168         *(d3_l->extra+off)-=carry; 
169        }
170       }
171      }
172      if(my_info->z_diff)
173      {
174       for(c=-1;c<=1;c++)
175       {
176        if(c!=0)
177        {
178         off=i+j*d3_l->max_x+((k+c+d3_l->max_z)%d3_l->max_z)*d3_l->max_x*d3_l->max_y;
179         carry=0;
180         if(!*(d3_l->status+off)&AMORPH) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
181         if(carry!=0)
182         {
183          *(d3_l->extra+off)-=carry;
184          *(d3_l->extra+offset)+=carry;
185         }
186        }
187       }
188      }  
189     } else
190     /* case not amorph: cryst <-> cryst diffusion */
191
192     /* test ! */
193     if(c_diff) {
194     /* */
195     {
196      for(c=-1;c<=1;c++) 
197      {
198       if(c!=0)
199       {
200        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;
201        carry=0;
202        if(!(*(d3_l->status+off)&AMORPH))
203        {
204         carry=(int)(my_info->dr_cc*(*(d3_l->extra+off)-*(d3_l->extra+offset))/2);
205         if(carry!=0)
206         {
207          *(d3_l->extra+offset)+=carry;
208          *(d3_l->extra+off)-=carry;
209         }
210        }
211       }
212      }
213      for(c=-1;c<=1;c++)
214      {
215       if(c!=0)
216       {
217        off=((i+c+d3_l->max_x)%d3_l->max_x)+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
218        carry=0;
219        if(!(*(d3_l->status+off)&AMORPH))
220        {
221         carry=(int)(my_info->dr_cc*(*(d3_l->extra+off)-*(d3_l->extra+offset))/2);
222         if(carry!=0)
223         {
224          *(d3_l->extra+offset)+=carry;
225          *(d3_l->extra+off)-=carry;
226         }
227        }
228       }
229      }
230     }
231     /* end test */
232     }
233     /* */
234    } /* for z */
235   } /* for y */
236  } /* for x */
237
238  } /* if step modulo diff_rate == 0 */
239
240  return 1;
241 }
242
243 int calc_pressure(d3_lattice *d3_l,int range)
244 {
245  int i,j,off;
246  double count,max=0;
247  int x,y,z;
248
249  for(x=0;x<d3_l->max_x;x++)
250  {
251   for(y=0;y<d3_l->max_y;y++)
252   {
253    for(z=0;z<d3_l->max_z;z++)
254    {
255     count=0;
256     for(i=-range;i<=range;i++)
257     {
258      for(j=-range;j<=range;j++)
259      {
260       if(i!=0 && j!=0)
261       {
262        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;
263        if(*(d3_l->status+off)&AMORPH) count+=((double)*(d3_l->extra+off))/(i*i+j*j);
264       }
265      }
266     }
267     if(count>max) max=count;
268    }
269   }
270  }
271
272  for(x=0;x<d3_l->max_x;x++)
273  {
274   for(y=0;y<d3_l->max_y;y++)
275   {
276    for(z=0;z<d3_l->max_z;z++)
277    {
278     count=0;
279     for(i=-range;i<=range;i++)
280     {
281      for(j=-range;j<=range;j++)
282      {
283       if(i!=0 && j!=0)
284       {
285        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;
286        if(*(d3_l->status+off)&AMORPH) count+=((double)*(d3_l->extra+off))/(i*i+j*j);
287       }
288      }
289     }
290     *(unsigned char *)(d3_l->v_ptr+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)=(unsigned char)(count*255/max);
291    }
292   }
293  }
294
295  return 1;
296 }
297
298 int calc_max_extra(d3_lattice *d3_l)
299 {
300  int x,y,z;
301  int off,max=0;
302
303  for(x=0;x<d3_l->max_x;x++)
304  {
305   for(y=0;y<d3_l->max_y;y++)
306   {
307    for(z=0;z<d3_l->max_z;z++)
308    {
309     off=x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
310     if(*(d3_l->extra+off)>max) max=*(d3_l->extra+off);
311    }
312   }
313  }
314
315  return max;
316 }
317
318 int write_bmp(d3_lattice *d3_l,int window,u32 x,u32 y,u32 z)
319 {
320  int fd,i,j,size=0,foo=0;
321  int width=0,height=0;
322  char bmpfile[MAX_CHARS];
323  char buf[128];
324
325  if(window==1)
326  {
327   sprintf(bmpfile,"x-z_%d.bmp",y);
328   foo=3*d3_l->max_x;
329   size=(foo+(4-foo%4))*d3_l->max_z;
330   width=d3_l->max_x;
331   height=d3_l->max_z;
332  }
333  if(window==2)
334  {
335   sprintf(bmpfile,"y-z_%d.bmp",x);
336   foo=3*d3_l->max_y;
337   size=(foo+(4-foo%4))*d3_l->max_z;
338   width=d3_l->max_y;
339   height=d3_l->max_z;
340  }
341  if(window==3)
342  {
343   sprintf(bmpfile,"x-y_%d.bmp",z);
344   foo=3*d3_l->max_x;
345   size=(foo+(4-foo%4))*d3_l->max_y;
346   width=d3_l->max_x;
347   height=d3_l->max_y;
348  }
349
350  if((fd=open(bmpfile,O_WRONLY|O_CREAT))<0)
351  {
352   puts("cannot open bmp file");
353   return -1;
354  }
355
356  /* bmpheader */
357  buf[0]='B'; /* std header start */
358  buf[1]='M';
359  buf[2]=(size+0x36)&0xff; /* file size */
360  buf[3]=(size+0x36)>>8;
361  memset(buf+4,0,6);
362  buf[10]=0x36; /* offset to data */
363  memset(buf+11,0,3);
364  buf[14]=0x28; /* length of bmp info header */
365  memset(buf+15,0,3);
366  buf[18]=width&0xff; /* width and height */
367  buf[19]=width>>8;
368  memset(buf+20,0,2);
369  buf[22]=height&0xff;
370  buf[23]=height>>8;
371  memset(buf+24,0,2);
372  buf[26]=1; /* # planes -> 1 */
373  buf[27]=0;
374  buf[28]=24; /* bits per pixel -> 2^24 (true color) */
375  buf[29]=0;
376  memset(buf+30,0,4); /* compression -> none */
377  buf[34]=size&0xff; /* data size */
378  buf[35]=size>>8;
379  memset(buf+36,0,2);
380  buf[38]=0x12; /* res: pixel/meter */
381  buf[39]=0x0b;
382  memset(buf+40,0,2);
383  buf[42]=0x12;
384  buf[43]=0x0b;
385  memset(buf+44,0,2); 
386  memset(buf+46,0,8); /* no colors, no important colors */
387
388  if(write(fd,buf,54)<54)
389  {
390   puts("failed writing bmp header");
391   return -1;
392  }
393  if(window==1)
394  {
395   for(j=0;j<d3_l->max_z;j++)
396   {
397    for(i=0;i<d3_l->max_x;i++)
398    {
399     if(*(d3_l->status+i+y*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3);
400     else memset(buf,0,3);
401     if(write(fd,buf,3)<3)
402     {
403      puts("failed writing rgb values to bmp file");
404      return-1;
405     }
406    }
407    if(foo%4)
408    {
409     memset(buf,0,4-foo%4);
410     if(write(fd,buf,4-foo%4)<4-foo%4)
411     {
412      puts("failed writing 4 byte ending");
413      return -1;
414     }
415    } 
416   }
417  }
418  if(window==2)
419  {
420   for(j=0;j<d3_l->max_z;j++)
421   {
422    for(i=0;i<d3_l->max_y;i++)
423    {
424     if(*(d3_l->status+x+i*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3);
425     else memset(buf,0,3);
426     if(write(fd,buf,3)<3)
427     {
428      puts("failed writing rgb values to bmp file");
429      return-1;
430     }
431    }
432    if(foo%4)
433    {
434     memset(buf,0,4-foo%4);
435     if(write(fd,buf,4-foo%4)<4-foo%4)
436     {
437      puts("failed writing 4 byte ending");
438      return -1;
439     }
440    }
441   }
442  }
443  if(window==3)
444  {
445   for(j=0;j<d3_l->max_y;j++)
446   {
447    for(i=0;i<d3_l->max_x;i++)
448    {
449     if(*(d3_l->status+i+(d3_l->max_y-j-1)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3);
450     else memset(buf,0,3);
451     if(write(fd,buf,3)<3)
452     {
453      puts("failed writing rgb values to bmp file");
454      return -1;
455     }
456    }
457    if(foo%4)
458    {
459     memset(buf,0,4-foo%4);
460     if(write(fd,buf,4-foo%4)<4-foo%4)
461     {
462      puts("failed writing 4 byte ending");
463      return -1;
464     }
465    }
466   }
467  }
468  close(fd);
469
470  return 1;
471 }
472
473 int save_to_file(char *sf,d3_lattice *d3_l,info *my_inf)
474 {
475  int sf_fd,c;
476
477  if((sf_fd=open(sf,O_WRONLY|O_CREAT))<0)
478  {
479   puts("cannot open save file");
480   return -1;
481  }
482  if(write(sf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
483  {
484   puts("failed saving d3 lattice struct");
485   return -1;
486  }
487  if(write(sf_fd,my_inf,sizeof(info))<sizeof(info))
488  {
489   puts("failed saving info struct");
490   return-1;
491  }
492  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
493  if(write(sf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
494  {
495   puts("failed saving status of d3 lattice sites");
496   return -1;
497  }
498  if(write(sf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
499  {
500   puts("failed saving sites concentration");
501   return -1;
502  }
503  close(sf_fd);
504
505  return 1;
506 }
507
508 int load_from_file(char *lf,d3_lattice *d3_l,info *my_inf)
509 {
510  int lf_fd,c;
511
512  if((lf_fd=open(lf,O_RDONLY))<0)
513  {
514   puts("cannot open load file");
515   return -1;
516  }
517  if(read(lf_fd,d3_l,sizeof(d3_lattice))<sizeof(d3_lattice))
518  {
519   puts("failed reading d3 lattice struct");
520   return -1;
521  }
522  if(read(lf_fd,my_inf,sizeof(info))<sizeof(info))
523  {
524   puts("failed reading info struct");
525   return-1;
526  }
527  c=d3_l->max_x*d3_l->max_y*d3_l->max_z;
528  if((d3_l->status=(unsigned char*)malloc(c*sizeof(unsigned char)))==NULL)
529  {
530   puts("cannot allocate status buffer");
531   return -1;
532  }
533  if((d3_l->extra=(int *)malloc(c*sizeof(int)))==NULL)
534  {
535   puts("cannot allocate concentration buffer");
536   return -1;
537  }
538  if(read(lf_fd,d3_l->status,c*sizeof(unsigned char))<c*sizeof(unsigned char))
539  {
540   puts("failed reading status of d3 lattice sites");
541   return -1;
542  }
543  if(read(lf_fd,d3_l->extra,c*sizeof(int))<c*sizeof(int))
544  {
545   puts("failed reading sites concentration");
546   return -1;
547  }
548  close(lf_fd);
549
550  return 1;
551 }
552
553 int convert_file(char *cf,d3_lattice *d3_l)
554 {
555  int x,y,z;
556  int c_fd;
557
558  if((c_fd=open(cf,O_WRONLY|O_CREAT))<0)
559  {
560   puts("cannot open convert file");
561   return -1;
562  }
563  dprintf(c_fd,"# created by nlsop (gnuplot format)\n");
564  for(x=0;x<d3_l->max_x;x++)
565  {
566   for(y=0;y<d3_l->max_y;y++)
567   {
568    for(z=0;z<d3_l->max_z;z++)
569    {
570     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);
571    }
572   }
573  }
574  close(c_fd);
575
576  return 1;
577 }
578
579 int get_c_ratio(double *c_ratio,char *pfile,info *my_info,d3_lattice *d3_l)
580 {
581  double all,a,b,d;
582  int i,k;
583  int p_fd;
584  unsigned char buf[32];
585  char *p;
586  unsigned char c;
587
588  if((p_fd=open(pfile,O_RDONLY))<0)
589  {
590   puts("cannot open profile file");
591   return -1;
592  }
593  k=1;
594  d=0;
595  all=0;
596  while(k)
597  {
598   for(i=0;i<32;i++)
599   {
600    k=read(p_fd,&c,1);
601    buf[i]=c;
602    if(c=='\n') break;
603   }
604   if(k)
605   {
606    p=strtok(buf," ");
607    a=atof(p)/10; /* nm */
608    p=strtok(NULL," ");
609    b=atof(p);
610    if(a>my_info->b_cd*CELL_LENGTH && a<(my_info->b_cd+d3_l->max_z)*CELL_LENGTH) d+=b;
611    all+=b;
612   }
613  }
614  *c_ratio=d/all;
615  close(p_fd);
616
617  return 1;
618 }
619
620 int main(int argc,char **argv)
621 {
622  u32 x,y,z,x_c,y_c,z_c;
623  int i,quit,escape,switchmode,nowait,bmp;
624  int refresh,resave;
625  int c_step;
626  char s_file[MAX_CHARS];
627  char s_file_tmp[MAX_CHARS];
628  char l_file[MAX_CHARS];
629  char c_file[MAX_CHARS];
630  char p_file[MAX_CHARS];
631  char convert;
632  char r_file[MAX_CHARS];
633 #ifdef USE_DFB_API
634  char xyz_txt[MAX_TXT];
635  char status_txt[MAX_TXT];
636  char conc_txt[MAX_TXT];
637  char steps_txt[MAX_TXT];
638  char cc_txt[MAX_TXT];
639  char a_txt[MAX_TXT];
640  char s_txt[MAX_TXT];
641  char ap_txt[MAX_TXT];
642  char el_txt[MAX_TXT];
643  char cd_txt[MAX_TXT];
644  char r_txt[MAX_TXT];
645  char cp_txt[MAX_TXT];
646  char zdiff_txt[MAX_TXT];
647  char diff_txt[MAX_TXT];
648  char dr_ac_txt[MAX_TXT];
649  char dr_cc_txt[MAX_TXT];
650  char dose_txt[MAX_TXT];
651  char mode_txt[MAX_TXT];
652  char *arg_v[MAX_ARGV];
653 #endif
654  d3_lattice d3_l;
655  info my_info;
656  unsigned char mode;
657  double c_ratio;
658 #ifdef USE_DFB_API
659  int max_extra;
660 #endif
661
662  d3_l.max_x=X;
663  d3_l.max_y=Y;
664  d3_l.max_z=Z;
665  my_info.steps=STEPS;
666  my_info.range=RANGE;
667  refresh=REFRESH;
668  resave=RESAVE;
669  my_info.z_diff=0;
670  c_diff=1;
671  my_info.a_el=A_EL;
672  my_info.b_el=B_EL;
673  my_info.a_cd=A_CD;
674  my_info.b_cd=B_CD;
675  my_info.a_ap=A_AP;
676  my_info.b_ap=B_AP;
677  my_info.a_cp=A_CP;
678  my_info.cc=CC;
679  my_info.dr_ac=DR_AC;
680  my_info.dr_cc=DR_CC;
681  my_info.diff_rate=DIFF_RATE;
682  nowait=0;
683  quit=0;
684  escape=0;
685  switchmode=0;
686  c_step=0;
687  strcpy(s_file,"");
688  strcpy(l_file,"");
689  strcpy(c_file,"");
690  strcpy(p_file,IMP_PROFILE);
691  convert=0;
692  strcpy(r_file,"");
693  mode=0;
694
695  amorph_count=0;
696  cryst_count=0;
697
698  for(i=1;i<argc;i++)
699  {
700   if(argv[i][0]=='-')
701   {
702    switch(argv[i][1])
703    {
704     case 'h':
705      usage();
706      return -1;
707     case 'n':
708      nowait=1;
709      break;
710     case 'a':
711      my_info.a_el=atof(argv[++i]);
712      break;
713     case 'b':
714      my_info.b_el=atof(argv[++i]);
715      break;
716     case 'x':
717      d3_l.max_x=atoi(argv[++i]);
718      break;
719     case 'y':
720      d3_l.max_y=atoi(argv[++i]);
721      break;
722     case 'z':
723      d3_l.max_z=atoi(argv[++i]);
724      break;
725     case 'Z':
726      my_info.z_diff=1;
727      break;
728     case 'i':
729      c_diff=0;
730      break;
731     case 's':
732      my_info.steps=atoi(argv[++i]);
733      break;
734     case 'd':
735      refresh=atoi(argv[++i]);
736      break;
737     case 'r':
738      my_info.range=atoi(argv[++i]);
739      break;
740     case 'f':
741      my_info.a_ap=atof(argv[++i]);
742      break;
743     case 'p':
744      my_info.b_ap=atof(argv[++i]);
745      break;
746     case 'F':
747      my_info.a_cp=atof(argv[++i]);
748      break;
749     case 'A':
750      my_info.a_cd=atof(argv[++i]);
751      break;
752     case 'B':
753      my_info.b_cd=atof(argv[++i]);
754      break;
755     case 'W':
756      resave=atoi(argv[++i]);
757      break;
758     case 'C':
759      strcpy(l_file,argv[++i]);
760      if(i<argc-1) if(argv[i+1][0]!='-') strcpy(c_file,argv[++i]);
761      convert=1;
762      break;
763     case 'D':
764      my_info.dr_ac=atof(argv[++i]);
765      break;
766     case 'c':
767      my_info.dr_cc=atof(argv[++i]);
768      break;
769     case 'e':
770      my_info.diff_rate=atoi(argv[++i]);
771      break;
772     case 'L':
773      strcpy(l_file,argv[++i]);
774      break;
775     case 'S':
776      strcpy(s_file,argv[++i]);
777      break;
778     case 'R':
779      strcpy(r_file,argv[++i]);
780      break;
781     case 'P':
782      strcpy(p_file,argv[++i]);
783      break;
784     case 'g':
785      strcpy(l_file,argv[++i]);
786      if(i<argc-1) if(argv[i+1][0]!='-') c_step=atoi(argv[++i]);
787      break;
788     default:
789      usage();
790      return -1;
791    }
792   } else usage();
793  }
794
795  x=d3_l.max_x/2-1;
796  y=d3_l.max_y/2-1;
797  z=d3_l.max_z/2-1;
798
799 #ifdef NODFB
800  if(!strcmp(s_file,""))
801  {
802   puts("NODFB defined, run with -S option");
803   return -1;
804  }
805 #endif
806
807  if(!strcmp(r_file,"")) rand_init(NULL);
808  else rand_init(r_file);
809
810  if(!strcmp(l_file,""))
811  {
812   i=d3_l.max_x*d3_l.max_y*d3_l.max_z;
813 #ifdef USE_DFB_API
814   d3_lattice_init(&argc,argv,&d3_l);
815 #endif
816   if((d3_l.status=(unsigned char *)malloc(i*sizeof(unsigned char)))==NULL)
817   {
818    puts("failed allocating status buffer");
819    return -1;
820   }
821   memset(d3_l.status,0,i*sizeof(unsigned char));
822   if((d3_l.extra=(int *)malloc(i*sizeof(int)))==NULL)
823   {
824    puts("failed allocating concentration buffer");
825    return -1;
826   }
827   memset(d3_l.extra,0,i*sizeof(int));
828  } else
829  {
830   load_from_file(l_file,&d3_l,&my_info);
831   if(convert) 
832   {   
833    if(!strcmp(c_file,"")) sprintf(c_file,"%s_gnuplot",l_file);
834    printf("converting file %s to %s\n",l_file,c_file);
835    convert_file(c_file,&d3_l);
836    puts("done");
837    return 1;
838   } 
839 #ifdef USE_DFB_API
840   else d3_lattice_init(&argc,argv,&d3_l);
841 #endif
842  }
843
844 #ifdef USE_DFB_API
845  d3_event_init(&d3_l);
846 #endif
847
848 #ifdef USE_DFB_API
849  strcpy(a_txt,"args:");
850  sprintf(s_txt,"steps: %d",my_info.steps);
851  sprintf(dose_txt,"dose: %.2fe+17 C/cm²",my_info.steps*1.0/(d3_l.max_x*d3_l.max_y*CELL_LENGTH*CELL_LENGTH*1000));
852  sprintf(r_txt,"pressure range: %d",my_info.range);
853  sprintf(ap_txt,"a_ap: %.3f  b_ap: %.3f",my_info.a_ap,my_info.b_ap);
854  sprintf(el_txt,"a_el: %.3f  b_el: %.3f",my_info.a_el,my_info.b_el);
855  sprintf(cd_txt,"a_cd: %.3f  b_cd: %.3f",my_info.a_cd,my_info.b_cd);
856  sprintf(cp_txt,"a_cp: %.4f",my_info.a_cp);
857  sprintf(dr_ac_txt,"a/c diffusion rate: %.4f",my_info.dr_ac);
858  sprintf(dr_cc_txt,"c/c diffusion rate: %.4f",my_info.dr_cc);
859  sprintf(zdiff_txt,"diffusion in z direction: %c",my_info.z_diff?'y':'n');
860  sprintf(diff_txt,"diffusion every %d steps",my_info.diff_rate);
861  strcpy(mode_txt,"view: a/c mode");
862  arg_v[1]=xyz_txt;
863  arg_v[2]=NULL;
864  arg_v[3]=status_txt;
865  arg_v[4]=conc_txt;
866  arg_v[5]=NULL;
867  arg_v[6]=mode_txt;
868  arg_v[7]=NULL;
869  arg_v[8]=steps_txt;
870  arg_v[9]=cc_txt;
871  arg_v[10]=NULL;
872  arg_v[11]=diff_txt;
873  arg_v[12]=zdiff_txt;
874  arg_v[13]=NULL;
875  arg_v[14]=a_txt;
876  arg_v[15]=s_txt;
877  arg_v[16]=dose_txt;
878  arg_v[17]=NULL;
879  arg_v[18]=r_txt;
880  arg_v[19]=ap_txt;
881  arg_v[20]=el_txt;
882  arg_v[21]=cd_txt;
883  arg_v[22]=cp_txt;
884  arg_v[23]=NULL;
885  arg_v[24]=dr_ac_txt;
886  arg_v[25]=dr_cc_txt;
887 #endif
888
889  if((!strcmp(l_file,""))||(c_step))
890  {
891   if(get_c_ratio(&c_ratio,p_file,&my_info,&d3_l)!=1)
892   {
893    puts("failed calculating ratio");
894    return -1;
895   }
896   i=(c_step?c_step:0);
897   while((i<my_info.steps) && (quit==0) && (escape==0))
898   {
899    x_c=get_rand(d3_l.max_x);
900    y_c=get_rand(d3_l.max_y);
901    z_c=get_rand_lgp(d3_l.max_z,my_info.a_el,my_info.b_el);
902    distrib_c(&d3_l,&my_info,i,c_ratio);
903    process_cell(&d3_l,x_c,y_c,z_c,&my_info);
904 #ifdef USE_DFB_API
905    if(i%refresh==0)
906    {
907     sprintf(xyz_txt,"x: %d  y: %d  z: %d",x+1,y+1,z+1);
908     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');
909     sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
910     sprintf(steps_txt,"step: %d  a_count: %d",i,amorph_count);
911     sprintf(cc_txt,"total c: %d  c_count: %d",my_info.cc,cryst_count);
912     d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,0);
913    }
914 #endif
915    if(i%resave==0 && strcmp(s_file,"") && resave!=0 && i!=0)
916    {
917     sprintf(s_file_tmp,"%s_%d_of_%d",s_file,i,my_info.steps);
918     save_to_file(s_file_tmp,&d3_l,&my_info);
919 #ifdef NODFB
920     printf("saved %s\n",s_file_tmp);
921 #endif
922    }
923    i++;
924   }
925  }
926
927  if(strcmp(s_file,""))
928  {
929    printf("saved %s\n",s_file);
930    save_to_file(s_file,&d3_l,&my_info);
931  }
932
933 #ifdef USE_DFB_API
934  /* allocating buffer for pressure values */
935  if((d3_l.v_ptr=malloc(d3_l.max_x*d3_l.max_y*d3_l.max_z*sizeof(unsigned char)))==NULL)
936  {
937   puts("cannot allocate buffer for pressure values");
938   return -1;
939  }
940  /* calc values */
941  calc_pressure(&d3_l,my_info.range);
942  max_extra=calc_max_extra(&d3_l);
943
944  while((quit==0) && (escape==0) && (nowait==0))
945  {
946   /* bahh! */
947   if(switchmode==0) mode=0;
948   if(switchmode==1) mode=1;
949   if(switchmode==2) mode=2;
950   /* end of bahh! */
951   sprintf(xyz_txt,"x: %d  y: %d  z: %d",x+1,y+1,z+1);
952   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');
953   sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y));
954   sprintf(steps_txt,"step: end");
955   sprintf(cc_txt,"total c: %d",my_info.cc);
956   if(switchmode==0) strcpy(mode_txt,"view: a/c mode");
957   if(switchmode==1) strcpy(mode_txt,"view: c conc mode");
958   if(switchmode==2) strcpy(mode_txt,"view: a pressure mode");
959   d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,max_extra);
960   bmp=0;
961   scan_event(&d3_l,&x,&y,&z,&quit,&escape,&switchmode,&bmp);
962   if(bmp) write_bmp(&d3_l,bmp,x,y,z);
963
964  }
965
966  d3_lattice_release(&d3_l);
967 #endif
968
969  return 1;
970 }