From 0ee31a02cd254d757c4750270d9423a049f25335 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 30 May 2003 15:38:59 +0000 Subject: [PATCH] added diffusion of c ions --- nlsop.c | 116 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 108 insertions(+), 8 deletions(-) diff --git a/nlsop.c b/nlsop.c index 87b551b..c4fc825 100644 --- a/nlsop.c +++ b/nlsop.c @@ -50,7 +50,10 @@ int usage(void) printf("-p \t pressure offset (default %f)\n",B_AP); 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",D_R); puts("-L \t load from file"); puts("-S \t save to file"); puts("-R \t read from random file"); @@ -83,22 +86,110 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,int r,double a,double b,int if(get_rand(URAND_MAX)<=p) { MAKE_AMORPH(thiz); - *t_c=*t_c+1-*conc; - } else *t_c+=1; + // *t_c=*t_c+1-*conc; + } // else *t_c+=1; } else { /* assume 1-p probability */ if(get_rand(URAND_MAX)>p) { MAKE_CRYST(thiz); - *t_c=*t_c+1+*conc; - } else *t_c+=1; + // *t_c=*t_c+1+*conc; + } // else *t_c+=1; } + + *t_c+=1; return 1; } -int distrib_c(d3_lattice *d3_l,int t_c,double a,double b) +int distrib_c(d3_lattice *d3_l,double d_r,double a,double b) +{ + u32 x,y,z; + int i,j,k,c; + int offset,off; + 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; + + /* diffusion in layer */ + for(i=0;imax_x;i++) + { + for(j=0;jmax_y;j++) + { + for(k=0;kmax_z;k++) + { + offset=i+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y; + /* case amorph */ + 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; + /* case neighbour not amorph */ + if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off))); + /* case neighbour amorph */ + else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; + if(carry!=0) + { + *(d3_l->extra+offset)+=carry; + *(d3_l->extra+off)-=carry; + } + } + } + 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; + /* case neighbour not amorph */ + if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(d_r*(*(d3_l->extra+off))); + /* case neighbour amorph */ + else carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; + if(carry!=0) + { + *(d3_l->extra+offset)+=carry; + *(d3_l->extra+off)-=carry; + } + } + } + } else + /* case not amorph */ + { + /* 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; + /* case neighbour not amorph */ + if(!(*(d3_l->status+off)&AMORPH)) + { + carry=(*(d3_l->extra+off)-*(d3_l->extra+offset))/2; + if(carry!=0) + { + *(d3_l->extra+offset)+=carry; + *(d3_l->extra+off)-=carry; + } + } + } + } + } + } /* for z */ + } /* for y */ + } /* for x */ + + return 1; +} + +int distrib_c_old(d3_lattice *d3_l,int t_c,double a,double b) { int i,j,k,total,area; double sum; @@ -264,6 +355,7 @@ int main(int argc,char **argv) char ap2_txt[MAX_TXT]; char cd2_txt[MAX_TXT]; char el2_txt[MAX_TXT]; + char dr_txt[MAX_TXT]; char *arg_v[MAX_ARGV]; d3_lattice d3_l; info my_info; @@ -281,6 +373,7 @@ int main(int argc,char **argv) my_info.a_ap=A_AP; my_info.b_ap=B_AP; my_info.cc=CC; + my_info.d_r=D_R; nowait=0; quit=0; escape=0; @@ -347,9 +440,14 @@ int main(int argc,char **argv) case 'B': my_info.b_cd=atof(argv[++i]); break; + /* case 'C': my_info.cc=atoi(argv[++i]); break; + */ + case 'D': + my_info.d_r=atof(argv[++i]); + break; case 'L': strcpy(l_file,argv[++i]); break; @@ -408,6 +506,7 @@ int main(int argc,char **argv) sprintf(el2_txt,"energy loss offset: %.2f",my_info.b_el); sprintf(cd_txt,"c distrib slope: %.2f",my_info.a_cd); sprintf(cd2_txt,"c distrib offset: %.2f",my_info.b_cd); + sprintf(dr_txt,"diffusion rate: %.2f",my_info.d_r); arg_v[1]=x_txt; arg_v[2]=y_txt; arg_v[3]=z_txt; @@ -431,6 +530,7 @@ int main(int argc,char **argv) arg_v[21]=el2_txt; arg_v[22]=cd_txt; arg_v[23]=cd2_txt; + arg_v[24]=dr_txt; if(!strcmp(l_file,"")) { @@ -440,7 +540,7 @@ 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_lgp(d3_l.max_z,my_info.a_el,my_info.b_el); - distrib_c(&d3_l,my_info.cc,my_info.a_cd,my_info.b_cd); + distrib_c(&d3_l,my_info.d_r,my_info.a_cd,my_info.b_cd); process_cell(&d3_l,x_c,y_c,z_c,my_info.range,my_info.a_ap,my_info.b_ap,&(my_info.cc)); if(i%refresh==0) { @@ -451,7 +551,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,23,arg_v); + d3_lattice_draw(&d3_l,x,y,z,24,arg_v); // scan_event(&d3_l,&x,&y,&z,&quit,&escape); } i++; @@ -469,7 +569,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)); strcpy(steps_txt,"step: end!"); sprintf(cc_txt,"total c: %d",my_info.cc); - d3_lattice_draw(&d3_l,x,y,z,23,arg_v); + d3_lattice_draw(&d3_l,x,y,z,24,arg_v); scan_event(&d3_l,&x,&y,&z,&quit,&escape); } -- 2.39.2