X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fnlsop.git;a=blobdiff_plain;f=nlsop.c;h=ff89e7cfbd5795c62666a0e6d31a385454c54e89;hp=2d093069f1c2ee5c24b421e3746b01d18f8efd69;hb=807dbc1ca4779971778796ffe0ccc97477166aba;hpb=126057c72d3dcf4cea9f5e46f36606289703bd17 diff --git a/nlsop.c b/nlsop.c index 2d09306..ff89e7c 100644 --- a/nlsop.c +++ b/nlsop.c @@ -29,6 +29,10 @@ #define MAKE_AMORPH(N) *(N)|=AMORPH #define MAKE_CRYST(N) *(N)&=~AMORPH +/* test globals - get removed or included in my_info struct later */ +int amorph_count; +int cryst_count; + int usage(void) { puts("usage:"); @@ -40,23 +44,17 @@ int usage(void) printf("-x \t # x cells (default %d)\n",X); printf("-y \t # y cells (default %d)\n",Y); printf("-z \t # z cells (default %d)\n",Z); - /* - printf("-X \t display x (default %d)\n",X/2-1); - printf("-Y \t display y (default %d)\n",Y/2-1); - printf("-Z \t display z (default %d)\n",Z/2-1); - */ printf("-s \t steps (default %d)\n",STEPS); printf("-d \t refresh display (default %d)\n",REFRESH); printf("-r \t amorphous influence range (default %d)\n",RANGE); printf("-f \t pressure = * 1/distance^2 (default %f)\n",A_AP); printf("-p \t pressure offset (default %f)\n",B_AP); + printf("-F \t proportionality constant between c conc and ability to get amorphous (default %f)\n",A_CP); printf("-A \t slope of linear c distribution (default %f)\n",A_CD); printf("-B \t linear c distribution offset (default %f)\n",B_CD); - /* - printf("-C \t initial c concentration (default %d)\n",CC); - */ printf("-D \t diffusion rate from cryst to amorph cells (default %f)\n",DR_AC); printf("-c \t diffusion rate in cryst cells (default %f)\n",DR_CC); + printf("-e \t do diffusion every steps (default %d)\n",DIFF_RATE); printf("-W \t write every steps to save file (default %d)\n",RESAVE); puts("-C \t convert file to gnuplot format"); puts("-L \t load from file"); @@ -90,22 +88,27 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info) } } p+=*conc*my_info->a_cp*URAND_MAX; - printf("debug: p = %f\n",p); if(!(*thiz&AMORPH)) { - if(get_rand(URAND_MAX)<=p) MAKE_AMORPH(thiz); + if(get_rand(URAND_MAX)<=p) + { + MAKE_AMORPH(thiz); + amorph_count++; + } } else { /* assume 1-p probability */ - if(get_rand(URAND_MAX)>p) MAKE_CRYST(thiz); + if(get_rand(URAND_MAX)>p) + { + MAKE_CRYST(thiz); + cryst_count++; + } } - my_info->cc+=1; - return 1; } -int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,char z_diff) */ +int distrib_c(d3_lattice *d3_l,info *my_info,int step,double c_ratio) { u32 x,y,z; int i,j,k,c; @@ -113,12 +116,18 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch int carry; /* put one c ion somewhere in the lattice */ - x=get_rand(d3_l->max_x); - y=get_rand(d3_l->max_y); - z=get_rand_lgp(d3_l->max_z,a,b); - *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1; + if(my_info->ccmax_x); + y=get_rand(d3_l->max_y); + z=get_rand_lgp(d3_l->max_z,my_info->a_cd,my_info->b_cd); + *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1; + (my_info->cc)++; + } + + if(step%my_info->diff_rate==0) + { - /* diffusion in layer */ for(i=0;imax_x;i++) { for(j=0;jmax_y;j++) @@ -126,24 +135,16 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch for(k=0;kmax_z;k++) { offset=i+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y; - /* case amorph */ + /* case amorph - amorph <-> cryst diffusion */ if(*(d3_l->status+offset)&AMORPH) { - /* look at neighbours and move c ions */ for(c=-1;c<=1;c++) { if(c!=0) { 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; carry=0; - /* case neighbour not amorph */ - if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off))); - /* case neighbour amorph */ - /* - * no diffusion between amorphous cells - * - else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; - */ + if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off))); if(carry!=0) { *(d3_l->extra+offset)+=carry; @@ -157,15 +158,7 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch { 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; carry=0; - /* case neighbour not amorph */ - if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off))); - /* case neighbour amorph */ - /* - * - * no diffusion between amorphous cells - * - else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; - */ + if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off))); if(carry!=0) { *(d3_l->extra+offset)+=carry; @@ -173,20 +166,35 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch } } } + if(my_info->z_diff) + { + for(c=-1;c<=1;c++) + { + if(c!=0) + { + 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; + carry=0; + if(!*(d3_l->status+off)&AMORPH) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off))); + if(carry!=0) + { + *(d3_l->extra+off)-=carry; + *(d3_l->extra+offset)+=carry; + } + } + } + } } else - /* case not amorph */ + /* case not amorph - cryst <-> cryst diffusion */ { - /* look at neighbours and move c ions */ for(c=-1;c<=1;c++) { if(c!=0) { 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; carry=0; - /* case neighbour not amorph */ if(!(*(d3_l->status+off)&AMORPH)) { - carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; + carry=(int)(my_info->dr_cc*(*(d3_l->extra+off)-*(d3_l->extra+offset))/2); if(carry!=0) { *(d3_l->extra+offset)+=carry; @@ -201,10 +209,9 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch { 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; carry=0; - /* case neighbour not amorph */ if(!(*(d3_l->status+off)&AMORPH)) { - carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; + carry=(int)(my_info->dr_cc*(*(d3_l->extra+off)-*(d3_l->extra+offset))/2); if(carry!=0) { *(d3_l->extra+offset)+=carry; @@ -213,33 +220,20 @@ int distrib_c(d3_lattice *d3_l,info *my_info) /* double d_r,double a,double b,ch } } } - /* cryst -> amorph diffusion in z direction */ - if(z_diff) - { - for(c=1;c>=-1;c--) - { - 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; - carry=0; - if(*(d3_l->status+off)&AMORPH) carry=(int)(d_r*(*(d3_l->extra+offset))); - if(carry!=0) - { - *(d3_l->extra+offset)-=carry; - *(d3_l->extra+off)+=carry; - } - } - } /* if z_diff */ } } /* for z */ } /* for y */ } /* for x */ + } /* if step modulo diff_rate == 0 */ + return 1; } int calc_pressure(d3_lattice *d3_l,int range) { int i,j,off; - double count; + double count,max=0; int x,y,z; for(x=0;xmax_x;x++) @@ -260,7 +254,30 @@ int calc_pressure(d3_lattice *d3_l,int range) } } } - *(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_VPTR); + if(count>max) max=count; + } + } + } + + for(x=0;xmax_x;x++) + { + for(y=0;ymax_y;y++) + { + for(z=0;zmax_z;z++) + { + count=0; + for(i=-range;i<=range;i++) + { + for(j=-range;j<=range;j++) + { + if(i!=0 && j!=0) + { + 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; + if(*(d3_l->status+off)&AMORPH) count+=((double)*(d3_l->extra+off))/(i*i+j*j); + } + } + } + *(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); } } } @@ -268,6 +285,151 @@ int calc_pressure(d3_lattice *d3_l,int range) return 1; } +int calc_max_extra(d3_lattice *d3_l) +{ + int x,y,z; + int off,max=0; + + for(x=0;xmax_x;x++) + { + for(y=0;ymax_y;y++) + { + for(z=0;zmax_z;z++) + { + off=x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y; + if(*(d3_l->extra+off)>max) max=*(d3_l->extra+off); + } + } + } + + return max; +} + +int write_bmp(d3_lattice *d3_l,int window,u32 x,u32 y,u32 z) +{ + int fd,i,j,size=0; + int width=0,height=0; + char bmpfile[MAX_CHARS]; + char buf[128]; + + if(window==1) + { + sprintf(bmpfile,"x-z_%d.bmp",y); + size=x*z; + width=x; + height=z; + } + if(window==2) + { + sprintf(bmpfile,"y-z_%d.bmp",x); + size=y*z; + width=y; + height=z; + } + if(window==3) + { + sprintf(bmpfile,"x-y_%d.bmp",z); + size=x*y; + width=x; + height=y; + } + + if((fd=open(bmpfile,O_WRONLY|O_CREAT))<0) + { + puts("cannot open bmp file"); + return -1; + } + + /* bmpheader */ + buf[0]='B'; /* std header start */ + buf[1]='M'; + buf[2]=(size+0x36)&0xff; /* file size */ + buf[3]=(size+0x36)>>8; + memset(buf+4,0,6); + buf[10]=0x36; /* offset to data */ + memset(buf+11,0,3); + buf[14]=0x28; /* length of bmp info header */ + memset(buf+15,0,3); + buf[18]=width&0xff; /* width and height */ + buf[19]=width>>8; + memset(buf+20,0,2); + buf[22]=height&0xff; + buf[23]=height>>8; + memset(buf+24,0,2); + buf[26]=1; /* # planes -> 1 */ + buf[27]=0; + buf[28]=24; /* bits per pixel -> 2^24 (true color) */ + buf[29]=0; + memset(buf+30,0,4); /* compression -> none */ + buf[34]=size&0xff; /* data size */ + buf[35]=size>>8; + memset(buf+36,0,2); + buf[38]=0x12; /* res: pixel/meter */ + buf[39]=0x0b; + memset(buf+40,0,2); + buf[42]=0x12; + buf[43]=0x0b; + memset(buf+44,0,2); + memset(buf+46,0,8); /* no colors, no important colors */ + + if(write(fd,buf,54)<54) + { + puts("failed writing bmp header"); + return -1; + } + if(window==1) + { + for(i=0;imax_x;i++) + { + for(j=0;jmax_z;j++) + { + if(*(d3_l->status+i+y*d3_l->max_x+j*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3); + else memset(buf,0,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + } + } + if(window==2) + { + for(i=0;imax_y;i++) + { + for(j=0;jmax_z;j++) + { + if(*(d3_l->status+x+i*d3_l->max_x+j*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3); + else memset(buf,0,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + } + } + if(window==3) + { + for(i=0;imax_x;i++) + { + for(j=0;jmax_y;j++) + { + if(*(d3_l->status+i+j*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&RED) memset(buf,0xff,3); + else memset(buf,0,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return -1; + } + } + } + } + close(fd); + + return 1; +} + int save_to_file(char *sf,d3_lattice *d3_l,info *my_inf) { int sf_fd,c; @@ -374,10 +536,51 @@ int convert_file(char *cf,d3_lattice *d3_l) return 1; } +int get_c_ratio(double *c_ratio,char *pfile,info *my_info,d3_lattice *d3_l) +{ + double all,a,b,d; + int i,k; + int p_fd; + unsigned char buf[32]; + char *p; + unsigned char c; + + if((p_fd=open(pfile,O_RDONLY))<0) + { + puts("cannot open profile file"); + return -1; + } + k=1; + d=0; + all=0; + while(k) + { + for(i=0;i<32;i++) + { + k=read(p_fd,&c,1); + buf[i]=c; + if(c=='\n') break; + } + if(k) + { + p=strtok(buf," "); + a=atof(p)/10; /* nm */ + p=strtok(NULL," "); + b=atof(p); + if(a>my_info->b_cd*CELL_LENGTH && a<(my_info->b_cd+d3_l->max_z)*CELL_LENGTH) d+=b; + all+=b; + } + } + *c_ratio=d/all; + close(p_fd); + + return 1; +} + int main(int argc,char **argv) { u32 x,y,z,x_c,y_c,z_c; - int i,quit,escape,switchmode,nowait; + int i,quit,escape,switchmode,nowait,bmp; int refresh,resave; char s_file[MAX_CHARS]; char s_file_tmp[MAX_CHARS]; @@ -398,17 +601,22 @@ int main(int argc,char **argv) char el_txt[MAX_TXT]; char cd_txt[MAX_TXT]; char r_txt[MAX_TXT]; - char ap2_txt[MAX_TXT]; - char cd2_txt[MAX_TXT]; - char el2_txt[MAX_TXT]; + char cp_txt[MAX_TXT]; + char zdiff_txt[MAX_TXT]; + char diff_txt[MAX_TXT]; char dr_ac_txt[MAX_TXT]; char dr_cc_txt[MAX_TXT]; + char dose_txt[MAX_TXT]; char mode_txt[MAX_TXT]; char *arg_v[MAX_ARGV]; #endif d3_lattice d3_l; info my_info; unsigned char mode; + double c_ratio; +#ifdef USE_DFB_API + int max_extra; +#endif d3_l.max_x=X; d3_l.max_y=Y; @@ -424,9 +632,11 @@ int main(int argc,char **argv) my_info.b_cd=B_CD; my_info.a_ap=A_AP; my_info.b_ap=B_AP; + my_info.a_cp=A_CP; my_info.cc=CC; my_info.dr_ac=DR_AC; my_info.dr_cc=DR_CC; + my_info.diff_rate=DIFF_RATE; nowait=0; quit=0; escape=0; @@ -439,6 +649,9 @@ int main(int argc,char **argv) strcpy(r_file,""); mode=0; + amorph_count=0; + cryst_count=0; + for(i=1;i