X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fnlsop.git;a=blobdiff_plain;f=nlsop.c;h=350868e536751af130b839c9f4f87f9140c24b81;hp=26ba461032b81746f4cc46b7ea7b4d22fc7627f7;hb=f8eab7046ff397059600f9f25181b38c101e6f81;hpb=9160d1de6debb4b30ea196f6c05287eca7cc0366 diff --git a/nlsop.c b/nlsop.c index 26ba461..350868e 100644 --- a/nlsop.c +++ b/nlsop.c @@ -10,7 +10,7 @@ * selforganizing lamella structure in the observed behaviour. * * refs: - * - J. K. N. Lindner. Habilationsschrift, Universitaet Augsburg. + * - J. K. N. Lindner. Habil.Schrift, Universitaet Augsburg. * - Maik Haeberlen. Diplomarbeit, Universitaet Augsburg. * * Copyright (C) 2004 Frank Zirkelbach @@ -55,8 +55,6 @@ int usage(void) puts("-n \t\t no user interaction"); puts("-Z \t\t cryst -> amorph c diffusion in z direction"); puts("-i \t\t no cryst to cryst diffusion"); - printf("-a \t slope of nuclear energy loss (default %f)\n",A_EL); - printf("-b \t nuclear energy loss offset (default %f)\n",B_EL); 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); @@ -66,8 +64,6 @@ int usage(void) 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("-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); @@ -85,17 +81,17 @@ int usage(void) return 1; } -int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info) +int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info,u32 nel_z) { unsigned char *thiz; int *conc; int i,j; int off; - double p; + double p,q; thiz=d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y; conc=d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y; - p=my_info->b_ap*URAND_MAX; + p=my_info->b_ap*nel_z; for(i=-(my_info->range);i<=my_info->range;i++) { for(j=-(my_info->range);j<=my_info->range;j++) @@ -114,6 +110,17 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info) } else { /* assume 1-p probability */ + /* also look for neighbours ! */ + q=(URAND_MAX-p)>0?URAND_MAX-p:0; + j=0; + j+=(*(d3_l->status+((x+d3_l->max_x+1)%d3_l->max_x)+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + j+=(*(d3_l->status+((x+d3_l->max_x-1)%d3_l->max_x)+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + j+=(*(d3_l->status+x+((y+1+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + j+=(*(d3_l->status+x+((y-1+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + j+=(*(d3_l->status+x+y*d3_l->max_x+((z+1+d3_l->max_z)%d3_l->max_z)*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + j+=(*(d3_l->status+x+y*d3_l->max_x+((z-1+d3_l->max_z)%d3_l->max_z)*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0; + + p+=((q/6)*j); if(get_rand(URAND_MAX)>p) MAKE_CRYST(thiz); } @@ -334,11 +341,14 @@ int write_ac_distr(d3_lattice *d3_l,int ac_distr) { int fd,x,y,z; int count,offset; - char file[16]; - - if(ac_distr==1) strcpy(file,"a.plot"); - if(ac_distr==2) strcpy(file,"c.plot"); - if(ac_distr==3) strcpy(file,"b.plot"); + char file[32]; + int si_count; + + si_count=d3_l->max_x*d3_l->max_y*SI_PER_VOLUME; + if(ac_distr==1) strcpy(file,"carbon_in_av.plot"); + if(ac_distr==2) strcpy(file,"carbon_in_cv.plot"); + if(ac_distr==3) strcpy(file,"carbon.plot"); + if(ac_distr==4) strcpy(file,"amorphous_volumes.plot"); if((fd=open(file,O_WRONLY|O_CREAT))<0) { @@ -359,9 +369,12 @@ int write_ac_distr(d3_lattice *d3_l,int ac_distr) if(ac_distr==2) if(!(*(d3_l->status+offset)&AMORPH)) count+=*(d3_l->extra+offset); if(ac_distr==3) count+=*(d3_l->extra+offset); + if(ac_distr==4) + if(*(d3_l->status+offset)&AMORPH) count+=1; } } - dprintf(fd,"%d %d\n",z,count); + if(ac_distr==4) dprintf(fd,"%d %d\n",z*CELL_LENGTH,count); + else dprintf(fd,"%d %f\n",z*CELL_LENGTH,100.0*count/si_count); } close(fd); @@ -423,7 +436,54 @@ int write_bmp(d3_lattice *d3_l,int window,u32 x,u32 y,u32 z,int max) width=d3_l->max_x; height=d3_l->max_y; } - + if(window==7) + { + sprintf(bmpfile,"x-z_a_cdistr_%d.bmp",y); + foo=3*d3_l->max_x; + size=(foo+(4-foo%4))*d3_l->max_z; + height=d3_l->max_z; + width=d3_l->max_x; + } + if(window==8) + { + sprintf(bmpfile,"y-z_a_cdistr_%d.bmp",x); + foo=3*d3_l->max_y; + size=(foo+(4-foo%4))*d3_l->max_z; + height=d3_l->max_z; + width=d3_l->max_y; + } + if(window==9) + { + sprintf(bmpfile,"x-y_a_cdistr_%d.bmp",z); + foo=3*d3_l->max_x; + size=(foo+(4-foo%4))*d3_l->max_y; + height=d3_l->max_y; + width=d3_l->max_x; + } + if(window==10) + { + sprintf(bmpfile,"x-z_c_cdistr_%d.bmp",y); + foo=3*d3_l->max_x; + size=(foo+(4-foo%4))*d3_l->max_z; + height=d3_l->max_z; + width=d3_l->max_x; + } + if(window==11) + { + sprintf(bmpfile,"y-z_c_cdistr_%d.bmp",x); + foo=3*d3_l->max_y; + size=(foo+(4-foo%4))*d3_l->max_z; + height=d3_l->max_z; + width=d3_l->max_y; + } + if(window==12) + { + sprintf(bmpfile,"x-y_c_cdistr_%d.bmp",z); + foo=3*d3_l->max_x; + size=(foo+(4-foo%4))*d3_l->max_y; + height=d3_l->max_y; + width=d3_l->max_x; + } if((fd=open(bmpfile,O_WRONLY|O_CREAT))<0) { puts("cannot open bmp file"); @@ -602,7 +662,7 @@ int write_bmp(d3_lattice *d3_l,int window,u32 x,u32 y,u32 z,int max) } if(window==6) { - for(j=0;jmax_z;j++) + for(j=0;jmax_y;j++) { for(i=0;imax_x;i++) { @@ -626,6 +686,168 @@ int write_bmp(d3_lattice *d3_l,int window,u32 x,u32 y,u32 z,int max) } } } + if(window==7) + { + for(j=0;jmax_z;j++) + { + for(i=0;imax_x;i++) + { + 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) sum=*(d3_l->extra+i+y*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } + if(window==8) + { + for(j=0;jmax_z;j++) + { + for(i=0;imax_x;i++) + { + 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) sum=*(d3_l->extra+x+i*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } + if(window==9) + { + for(j=0;jmax_y;j++) + { + for(i=0;imax_x;i++) + { + 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) sum=*(d3_l->extra+i+(d3_l->max_y-j-1)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } + if(window==10) + { + for(j=0;jmax_z;j++) + { + for(i=0;imax_x;i++) + { + 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)) sum=*(d3_l->extra+i+y*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } + if(window==11) + { + for(j=0;jmax_z;j++) + { + for(i=0;imax_x;i++) + { + 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)) sum=*(d3_l->extra+x+i*d3_l->max_x+(d3_l->max_z-j-1)*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } + if(window==12) + { + for(j=0;jmax_y;j++) + { + for(i=0;imax_x;i++) + { + 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)) sum=*(d3_l->extra+i+(d3_l->max_y-j-1)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y); + else sum=0; + sum=sum*255/max; + memset(buf,(unsigned char)sum,3); + if(write(fd,buf,3)<3) + { + puts("failed writing rgb values to bmp file"); + return-1; + } + } + if(foo%4) + { + memset(buf,0,4-foo%4); + if(write(fd,buf,4-foo%4)<4-foo%4) + { + puts("failed writing 4 byte ending"); + return -1; + } + } + } + } close(fd); return 1; @@ -668,7 +890,8 @@ int save_to_file(char *sf,d3_lattice *d3_l,info *my_inf) int load_from_file(char *lf,d3_lattice *d3_l,info *my_inf) { - int lf_fd,c; + + int lf_fd,c,pos,end,data,data_len,strip; if((lf_fd=open(lf,O_RDONLY))<0) { @@ -680,11 +903,31 @@ int load_from_file(char *lf,d3_lattice *d3_l,info *my_inf) puts("failed reading d3 lattice struct"); return -1; } - if(read(lf_fd,my_inf,sizeof(info))max_x*d3_l->max_y*d3_l->max_z; + data_len=data*(sizeof(int)+sizeof(unsigned char)); + printf("there are %d volumes so we need %d of bytes\n",data,data_len); + end=lseek(lf_fd,0,SEEK_END); + c=end-pos-data_len; + printf("end: %d => length: %d => guessed info size: %d bytes\n",end,end-pos,c); + strip=sizeof(info)-c; + printf("as real programs info size is %d, we strip %d bytes\n",sizeof(info),strip); + lseek(lf_fd,pos,SEEK_SET); + c=sizeof(info); + if(strip>0) c-=strip; + if(c<0) + { + puts("info smaller then strip size"); + return -1; + } + if(read(lf_fd,my_inf,c)0) memset(my_inf+c,0,strip); + if(strip<0) lseek(lf_fd,(-1*strip),SEEK_CUR); c=d3_l->max_x*d3_l->max_y*d3_l->max_z; if((d3_l->status=(unsigned char*)malloc(c*sizeof(unsigned char)))==NULL) { @@ -768,7 +1011,7 @@ int get_c_ratio(double *c_ratio,char *pfile,info *my_info,d3_lattice *d3_l) 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; + if(a<=d3_l->max_z*CELL_LENGTH) d+=b; all+=b; } } @@ -817,7 +1060,7 @@ u32 get_reject_graph(info *my_info,d3_lattice *d3_l,char *file,u32 *graph) { if(a>d3_l->max_z*CELL_LENGTH) k=0; else { - graph[(int)(a/CELL_LENGTH)]=(int)(URAND_MAX*b); + graph[(int)(a/CELL_LENGTH)]=(int)(URAND_MAX/100*b); flag[(int)(a/CELL_LENGTH)]=1; } } @@ -838,7 +1081,7 @@ u32 get_reject_graph(info *my_info,d3_lattice *d3_l,char *file,u32 *graph) { free(flag); #ifdef DEBUG_INTERPOL_PROFILE - printf("debug: (interpolated c profile)\n"); + printf("debug: %s (interpolated profile)\n",file); for(i=0;imax_z;i++) printf("%d %d\n",i,graph[i]); #endif @@ -867,11 +1110,10 @@ int main(int argc,char **argv) char cc_txt[MAX_TXT]; char a_txt[MAX_TXT]; char s_txt[MAX_TXT]; - char ap_txt[MAX_TXT]; - char el_txt[MAX_TXT]; - char cd_txt[MAX_TXT]; + char ballistic_txt[MAX_TXT]; + char carbon_txt[MAX_TXT]; + char stress_txt[MAX_TXT]; char r_txt[MAX_TXT]; - char cp_txt[MAX_TXT]; char zdiff_txt[MAX_TXT]; char diff_txt[MAX_TXT]; char dr_ac_txt[MAX_TXT]; @@ -892,6 +1134,7 @@ int main(int argc,char **argv) u32 *c_profile; u32 *n_e_loss; u32 ne_max,ip_max; + u32 nel_z; d3_l.max_x=X; d3_l.max_y=Y; @@ -902,10 +1145,6 @@ int main(int argc,char **argv) resave=RESAVE; my_info.z_diff=0; my_info.c_diff=1; - my_info.a_el=A_EL; - my_info.b_el=B_EL; - my_info.a_cd=A_CD; - my_info.b_cd=B_CD; my_info.a_ap=A_AP; my_info.b_ap=B_AP; my_info.a_cp=A_CP; @@ -943,12 +1182,6 @@ int main(int argc,char **argv) case 'n': nowait=1; break; - case 'a': - my_info.a_el=atof(argv[++i]); - break; - case 'b': - my_info.b_el=atof(argv[++i]); - break; case 'x': d3_l.max_x=atoi(argv[++i]); break; @@ -983,12 +1216,6 @@ int main(int argc,char **argv) case 'F': my_info.a_cp=atof(argv[++i]); break; - case 'A': - my_info.a_cd=atof(argv[++i]); - break; - case 'B': - my_info.b_cd=atof(argv[++i]); - break; case 'W': resave=atoi(argv[++i]); break; @@ -1096,10 +1323,9 @@ int main(int argc,char **argv) sprintf(s_txt,"steps: %d",my_info.steps); 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)); sprintf(r_txt,"pressure range: %d",my_info.range); - sprintf(ap_txt,"a_ap: %f b_ap: %f",my_info.a_ap,my_info.b_ap); - sprintf(el_txt,"a_el: %f b_el: %f",my_info.a_el,my_info.b_el); - sprintf(cd_txt,"a_cd: %f b_cd: %f",my_info.a_cd,my_info.b_cd); - sprintf(cp_txt,"a_cp: %f",my_info.a_cp); + sprintf(stress_txt,"stress term: %f",my_info.a_ap); + sprintf(ballistic_txt,"ballistic term: %f",my_info.b_ap); + sprintf(carbon_txt,"carbon term: %f",my_info.a_cp); sprintf(dr_ac_txt,"a/c diffusion rate: %f",my_info.dr_ac); if(my_info.c_diff!=0) sprintf(dr_cc_txt,"c/c diffusion rate: %f",my_info.dr_cc); else sprintf(dr_cc_txt,"c/c diffusion rate: none"); @@ -1122,19 +1348,33 @@ int main(int argc,char **argv) arg_v[12]=diff_txt; arg_v[13]=zdiff_txt; arg_v[14]=r_txt; - arg_v[15]=ap_txt; - arg_v[16]=el_txt; - arg_v[17]=cd_txt; - arg_v[18]=cp_txt; - arg_v[19]=dr_ac_txt; - arg_v[20]=dr_cc_txt; - arg_v[21]=hpi_txt; - arg_v[22]=csat_txt; + arg_v[15]=ballistic_txt; + arg_v[16]=carbon_txt; + arg_v[17]=stress_txt; + arg_v[18]=dr_ac_txt; + arg_v[19]=dr_cc_txt; + arg_v[20]=hpi_txt; + arg_v[21]=csat_txt; + arg_v[22]=NULL; arg_v[23]=NULL; arg_v[24]=NULL; arg_v[25]=NULL; #endif + /* compute graphs for random number rejection method */ + if((c_profile=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) + { + puts("failed allocating memory for carbon profile graph"); + return -1; + } + if((n_e_loss=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) + { + puts("failed allocating memory for nuclear energy loss graph"); + return -1; + } + ip_max=get_reject_graph(&my_info,&d3_l,p_file,c_profile); + ne_max=get_reject_graph(&my_info,&d3_l,n_e_file,n_e_loss); + if((!strcmp(l_file,""))||(c_step)) { /* calculate ratio of c_simwindow / c_total */ @@ -1143,21 +1383,9 @@ int main(int argc,char **argv) puts("failed calculating ratio"); return -1; } - /* compute graphs for random number rejection method */ - if((c_profile=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) - { - puts("failed allocating memory for carbon profile graph"); - return -1; - } - if((n_e_loss=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) - { - puts("failed allocating memory for nuclear energy loss graph"); - return -1; - } - ip_max=get_reject_graph(&my_info,&d3_l,p_file,c_profile); - ne_max=get_reject_graph(&my_info,&d3_l,n_e_file,n_e_loss); #ifdef DEBUG_RAND + i=0; while(1) { #ifdef DEBUG_CP @@ -1169,6 +1397,8 @@ int main(int argc,char **argv) #ifdef DEBUG_NORM printf("%d\n",get_rand(d3_l.max_z)); #endif + if(i==10000000) return 1; + i++; } #endif @@ -1179,8 +1409,10 @@ int main(int argc,char **argv) { x_c=get_rand(d3_l.max_x); y_c=get_rand(d3_l.max_y); - z_c=get_rand_reject(d3_l.max_z,ne_max,n_e_loss); - process_cell(&d3_l,x_c,y_c,z_c,&my_info); + // z_c=get_rand_reject(d3_l.max_z,ne_max,n_e_loss); + z_c=get_rand(d3_l.max_z); + nel_z=URAND_MAX*(1.0*n_e_loss[z_c]/ne_max); + process_cell(&d3_l,x_c,y_c,z_c,&my_info,nel_z); } distrib_c(&d3_l,&my_info,i,c_ratio,ip_max,c_profile); #ifdef USE_DFB_API @@ -1191,7 +1423,7 @@ int main(int argc,char **argv) sprintf(conc_txt,"conc: %d",*(d3_l.extra+x+y*d3_l.max_x+z*d3_l.max_x*d3_l.max_y)); sprintf(steps_txt,"step: %d",i); sprintf(cc_txt,"total c: %d",my_info.cc); - d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,0); + d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,0,NULL,0,NULL,0); } #endif if(i%resave==0 && strcmp(s_file,"") && resave!=0 && i!=0) @@ -1229,6 +1461,7 @@ int main(int argc,char **argv) if(switchmode==0) mode=0; if(switchmode==1) mode=1; if(switchmode==2) mode=2; + if(switchmode==3) mode=3; /* end of bahh! */ sprintf(xyz_txt,"x: %d y: %d z: %d",x+1,y+1,z+1); 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'); @@ -1238,7 +1471,8 @@ int main(int argc,char **argv) if(switchmode==0) strcpy(mode_txt,"view: a/c mode"); if(switchmode==1) strcpy(mode_txt,"view: c conc mode"); if(switchmode==2) strcpy(mode_txt,"view: a pressure mode"); - d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,max_extra); + if(switchmode==3) strcpy(mode_txt,"view: a/c + profiles mode"); + d3_lattice_draw(&d3_l,x,y,z,25,arg_v,mode,max_extra,c_profile,ip_max,n_e_loss,ne_max); bmp=0; ac_distr=0; scan_event(&d3_l,&x,&y,&z,&quit,&escape,&switchmode,&bmp,&ac_distr);