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