#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);
}