From: hackbard Date: Wed, 2 Jun 2004 15:12:26 +0000 (+0000) Subject: new default values, first attempts of rejection method for random n generation X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=de3ae9b29a64eaf297b44665022665962eba4a14;p=physik%2Fnlsop.git new default values, first attempts of rejection method for random n generation --- diff --git a/nlsop.c b/nlsop.c index 5d44239..0104695 100644 --- a/nlsop.c +++ b/nlsop.c @@ -130,7 +130,9 @@ int distrib_c(d3_lattice *d3_l,info *my_info,int step,double c_ratio) { x=get_rand(d3_l->max_x); y=get_rand(d3_l->max_y); + // printf("nd: %d %d\n",x,y); z=get_rand_lgp(d3_l->max_z,my_info->a_cd,my_info->b_cd); + // printf("ld: %d\n",z); *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1; (my_info->cc)++; } @@ -692,6 +694,42 @@ int get_c_ratio(double *c_ratio,char *pfile,info *my_info,d3_lattice *d3_l) return 1; } +int get_reject_graph(info *my_info,d3_lattice *d3_l,char *file,u32 *graph) { + double a,b; + int i,j,k; + int fd; + char buf[32],*p; + + if((fd=open(file,O_RDONLY))<0) + { + puts("cannot open file to calculate rejection graph"); + return -1; + } + memset(graph,0,d3_l->max_z*sizeof(u32)); + /* get fixpoints */ + k=1; + while(k) + { + for(i=0;i<32;i++) + { + k=read(fd,&buf[i],1); + if((buf[i]=='\n')||(k==0)) break; + } + if(k) + { + p=strtok(buf," "); + a=atof(p)/10; /* nm */ + p=strtok(NULL," "); + b=atof(p); + if(a>d3_l->max_z*CELL_LENGTH) k=0; + else graph[(int)(a/CELL_LENGTH)]=(int)(URAND_MAX*b); + } + } + /* do interpolation here! */ + + return 1; +} + int main(int argc,char **argv) { u32 x,y,z,x_c,y_c,z_c; @@ -927,10 +965,10 @@ 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: %.3f b_ap: %.3f",my_info.a_ap,my_info.b_ap); + sprintf(ap_txt,"a_ap: %.4f b_ap: %.3f",my_info.a_ap,my_info.b_ap); sprintf(el_txt,"a_el: %.3f b_el: %.3f",my_info.a_el,my_info.b_el); sprintf(cd_txt,"a_cd: %.3f b_cd: %.3f",my_info.a_cd,my_info.b_cd); - sprintf(cp_txt,"a_cp: %.4f",my_info.a_cp); + sprintf(cp_txt,"a_cp: %.5f",my_info.a_cp); sprintf(dr_ac_txt,"a/c diffusion rate: %.4f",my_info.dr_ac); if(my_info.c_diff!=0) sprintf(dr_cc_txt,"c/c diffusion rate: %.4f",my_info.dr_cc); else sprintf(dr_cc_txt,"c/c diffusion rate: none"); diff --git a/nlsop.h b/nlsop.h index bb3307f..10d2a43 100644 --- a/nlsop.h +++ b/nlsop.h @@ -26,10 +26,10 @@ typedef struct __info #define Y 64 #define Z 100 -#define STEPS 30000000 +#define STEPS 100000000 #define RANGE 5 #define REFRESH 100000 -#define RESAVE 1000000 +#define RESAVE 10000000 #define A_CD 1. #define B_CD .0 @@ -38,12 +38,12 @@ typedef struct __info #define DR_AC .5 #define DR_CC .2 -#define DIFF_RATE 100 +#define DIFF_RATE 1000 -#define A_AP .003 +#define A_AP .0005 #define B_AP .0 -#define A_CP .0001 +#define A_CP .00001 #define MAX_CHARS 128 #define MAX_TXT 32 diff --git a/random.c b/random.c index 6690531..ec99a33 100644 --- a/random.c +++ b/random.c @@ -73,3 +73,18 @@ u32 get_rand_lgp(u32 max,double a,double b) return((u32)(1.0*max*(-1.0*b+sqrt(b*b+2*a*((b+a/2)*get_rand(URAND_MAX)/((long long unsigned int)URAND_MAX+1))))/a)); } +u32 get_rand_reject(u32 max_x,u32 max_y,u32 *graph) +{ + u32 x,y; + unsigned char ok; + + ok=0; + while(!ok) + { + x=get_rand(max_x); + y=get_rand(max_y); + if(y<=graph[x]) ok=1; + } + + return y; +} diff --git a/random.h b/random.h index e443f81..a578afb 100644 --- a/random.h +++ b/random.h @@ -17,6 +17,7 @@ int rand_init(char *rf); int rand_close(void); u32 get_rand(u32 max); u32 get_rand_lgp(u32 max,double a,double b); +u32 get_rand_reject(u32 max_x,u32 max_y,u32 *graph); #define BUFSIZE (16*1048576) /* 64 MByte */ #define URAND_MAX 0xffffffff