我想让程序实现通过产生 0 到 1 的随机数，从而有比阈值 0.026。（就是那个 s）小的随机数，使得 N_walkers 递减，可无法实现此功能。
#include <stdio.h>#include <stdlib.h>
#include <math.h>
#include <time.h>
#define re 0.000001
#define PI 3.1415926
void main()
{FILE *fp;
int l,w;
double s,t,N_decay;
double rho=20*(1E-6),D=2.5*(1E-9);
double p;
int a,b,c;
int i,j,k;
int index;
int count=0,response=100;
double rand_p;
double theta,phi;
int N_init=100,N_pore=0,N_walkers=0;
static int array[400][3];
static int cube[100][100][100];
static double X_old[100],Y_old[100],Z_old[100],X_new[100],Y_new[100],Z_new[100];
s=0.5*re;
p=(2*rho*s)/(3*D);
fp=fopen("D:\\test.txt","w");
for(i=0;i<N_init-1;i++)
{for(j=0;j<N_init-1;j++)
{for(k=0;k<N_init-1;k++)
{cube[i][j][k]=1;
if(pow((i-50),2.0)+pow((j-50),2.0)+pow((k-50),2.0)<=100)
{cube[i][j][k]=0;
array[N_pore][0]=i;
array[N_pore][1]=j;
array[N_pore][2]=k;
N_pore++;
}
}
}
}
N_walkers=N_pore;
if(N_init<N_pore)
{for(i=0;i<N_init-1;i++)
{index=(int) floor(rand()/(double)(RAND_MAX+1.0)*N_pore);
X_old[i]=((double)(array[index][0])+rand()/(RAND_MAX+1.0))*re;
Y_old[i]=((double)(array[index][1])+rand()/(RAND_MAX+1.0))*re;
Z_old[i]=((double)(array[index][2])+rand()/(RAND_MAX+1.0))*re;
}
N_walkers=N_init;
}
else
{N_init=N_pore;
for(w=0;w<N_pore-1;w++)
{X_old[w]=((double)(array[w][0])+rand()/(RAND_MAX+1.0))*re;
Y_old[w]=((double)(array[w][1])+rand()/(RAND_MAX+1.0))*re;
Z_old[w]=((double)(array[w][2])+rand()/(RAND_MAX+1.0))*re;
}
N_walkers=N_pore;
}
do
{for(i=0;i<N_walkers-1;i++)
{theta=(rand()/(RAND_MAX+1.0))*2*PI;
phi=(rand()/(RAND_MAX+1.0))*PI;
X_new[i]=X_old[i]+s*sin(phi)*cos(theta);
Y_new[i]=Y_old[i]+s*sin(phi)*sin(theta);
Z_new[i]=Z_old[i]+s*cos(phi);
a=(int)floor(X_new[i]/re);
b=(int)floor(Y_new[i]/re);
c=(int)floor(Z_new[i]/re);
if(cube[a][b][c]==0)
{X_old[i]=X_new[i];
Y_old[i]=Y_new[i];
Z_old[i]=Z_new[i];
}
else
{srand((unsigned)time(0));
rand_p=rand()/(RAND_MAX+1.0);
if(rand_p<p)
{for(l=i;l<N_walkers-1;l++)
{X_old[l]=X_old[l+1];
Y_old[l]=Y_old[l+1];
Z_old[l]=Z_old[l+1];
}
N_walkers--;
}
}
}
count++;
t=(pow(s,2)/(6*D))*count;
if(count%response==0)
{N_decay=N_walkers/N_init;
printf("%lf 弛豫时间是%lf 衰减率是%lf\n",rand_p,t,N_decay);
fprintf (fp,"弛豫时间是%lf 衰减率是%lf \n",t,N_decay);
}
}while((double)(N_walkers/N_init)>0.01);
fclose(fp);
}
