nel_z bugfix ; new runs

[physik/nlsop.git] / nlsop.c

diff --git a/nlsop.c b/nlsop.c
index 7603876..deb65c9 100644 (file)
--- a/nlsop.c
+++ b/nlsop.c
@@ -87,11 +87,11 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info,u32 nel_z)
  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*nel_z;
+ 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++)
@@ -110,6 +110,17 @@ int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info,u32 nel_z)
  } 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);
  }
  
@@ -1158,7 +1169,7 @@ int main(int argc,char **argv)
     y_c=get_rand(d3_l.max_y);
     // 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*n_e_loss[z_c])/ne_max;
+    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);