#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 q=0.01,D=0.02;
 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;
 int array[999][3];
 int cube[100][100][100];
 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*q*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)+pow((j-50),2)+pow((k-50),2)<=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()*N_pore);
  X_old[i]=(array[index][0]+rand())*re;
  Y_old[i]=(array[index][1]+rand())*re;
  Z_old[i]=(array[index][2]+rand())*re;  
 }
 N_walkers=N_init;
 }
 else
 {N_init=N_pore;
  for(w=0;w<N_pore-1;w++)
  {X_old[w]=(array[w][0]+rand())*re;
   Y_old[w]=(array[w][1]+rand())*re;
   Z_old[w]=(array[w][2]+rand())*re;
  }
  N_walkers=N_pore;
 }
 do
 {for(i=0;i<N_walkers;i++)
 {theta=rand()*2*PI;
  phi=rand()*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
  {rand_p=rand();
   if(rand_p<p)
   {for(l=i;l<N_walkers;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=s*s/6/D*count;
 if(count%response==0)
 {N_decay=N_walkers/N_init;
  fprintf (fp,"弛豫时间是%lf 衰减率是%lf \n",t,N_decay);
 }
 }while(N_walkers/N_init>0.01);
 fclose(fp);
}

用code blocks编译后结果为Process terminated with status -1073741571

    int cube[100][100][100];
数组太大，栈上分配失败，使程序无法运行。可改为 malloc，用指针存取。或偷懒改 int 为 char，可减少3/4所需内存。

此外，
index=(int) floor(rand()*N_pore);
rand() 返回 0 ~ RAND_MAX （一般是32767）以内随机数，这有可能得到很大的 index，造成接下来数组越界。猜测语句的意思，可能应改为：
index=(int) floor(rand()/(double)(RAND_MAX+1)*N_pore);
其它调用 rand() 同样要注意这一点。


[ 本帖最后由 top398 于 2014-5-25 17:11 编辑 ]

谢谢您了，随机数部分我修改了。上面如果改做malloc会不会对程序有很大变动啊？