added c ratio calculation

diff --git a/nlsop.c b/nlsop.c
index 3be95f930f6658e31d70dad91343804084e21247..d57cc04cb467fa2f4909ce935bb01ee89ca6097d 100644 (file)
--- a/nlsop.c
+++ b/nlsop.c
@@ -94,12 +94,10 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info)
   if(get_rand(URAND_MAX)>p) MAKE_CRYST(thiz);
  }
  
- my_info->cc+=1;
-
  return 1;
 }
 
-int distrib_c(d3_lattice *d3_l,info *my_info,int step)
+int distrib_c(d3_lattice *d3_l,info *my_info,int step,double c_ratio)
 {
  u32 x,y,z;
  int i,j,k,c;
@@ -107,10 +105,14 @@ int distrib_c(d3_lattice *d3_l,info *my_info,int step)
  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,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;
+ if(my_info->cc<c_ratio*step)
+ {
+  x=get_rand(d3_l->max_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)
  {
@@ -355,6 +357,46 @@ 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;
+
+ return 1;
+}
+
 int main(int argc,char **argv)
 {
  u32 x,y,z,x_c,y_c,z_c;
@@ -390,6 +432,7 @@ int main(int argc,char **argv)
  d3_lattice d3_l;
  info my_info;
  unsigned char mode;
+ double c_ratio;
 
  d3_l.max_x=X;
  d3_l.max_y=Y;
@@ -592,27 +635,33 @@ int main(int argc,char **argv)
  arg_v[12]=zdiff_txt;
  arg_v[13]=NULL;
  arg_v[14]=a_txt;
- arg_v[15]=s_txt;
- arg_v[16]=r_txt;
- arg_v[17]=ap_txt;
- arg_v[18]=el_txt;
- arg_v[19]=cd_txt;
- arg_v[20]=cp_txt;
- arg_v[21]=dr_ac_txt;
- arg_v[22]=dr_cc_txt;
+ arg_v[15]=NULL;
+ arg_v[16]=s_txt;
+ arg_v[17]=r_txt;
+ arg_v[18]=NULL;
+ arg_v[19]=ap_txt;
+ arg_v[20]=el_txt;
+ arg_v[21]=cd_txt;
+ arg_v[22]=cp_txt;
  arg_v[23]=NULL;
- arg_v[24]=NULL;
+ arg_v[24]=dr_ac_txt;
+ arg_v[25]=dr_cc_txt;
 #endif
 
  if(!strcmp(l_file,""))
  {
+  if(get_c_ratio(&c_ratio,p_file,&my_info,&d3_l)!=1)
+  {
+   puts("failed calculating ratio");
+   return -1;
+  }
   i=0;
   while((i<my_info.steps) && (quit==0) && (escape==0))
   {
    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,i);
+   distrib_c(&d3_l,&my_info,i,c_ratio);
    process_cell(&d3_l,x_c,y_c,z_c,&my_info);
 #ifdef USE_DFB_API
    if(i%refresh==0)

diff --git a/nlsop.h b/nlsop.h
index 2f25a6a3993eadb78904443a2a375f5a140f6524..9fec8cb7de87227617aa5d8efab53b2bae417340 100644 (file)
--- a/nlsop.h
+++ b/nlsop.h
@@ -52,5 +52,6 @@ typedef struct __info
 #define MAX_VPTR 50
 
 #define IMP_PROFILE "imp-profile.txt"
+#define CELL_LENGTH 3 /* nm */
 
 #endif /* NLSOP_H */