/*********************************************************************/
/* */
/* Program: randnum.c */
/* Purpose: To generate random radii for particles sorted from */
/* largest to smallest */
/* Contains: genrand (random radii from 0.0 to 1.0) */
/* gensieve (radii based on sieve size classes) */
/* Programmer: Dale P. Bentz */
/* NIST */
/* Building 226 Room B-350 */
/* Gaithersburg, MD 20899-8621 */
/* Phone: (301) 975-5865 Fax: (301) 990-6891 */
/* E-mail: dale.bentz@nist.gov */
/* WWW page: http://ciks.cbt.nist.gov/~bentz */
/* */
/*********************************************************************/
/* Uniform random deviates can be transformed to other */
/* distributions (Weibull, etc.) as needed */
/* nsort is number of particle radii to be generated */
void genrand(nsort)
int nsort;
{
int ic,jc; /* loop index variables */
float temp; /* interchange variable */
/* Generate the random deviates */
for(ic=0;ic<nsort;ic++){
randrad[ic]=ran1(seed);
}
/* Sort the random deviates from largest to smallest */
/* for use in particle placement */
/* Simple bubble sort */
for(ic=0;ic<(nsort-1);ic++){
for(jc=(ic+1);jc<nsort;jc++){
if(randrad[ic]<randrad[jc]){
temp=randrad[ic];
randrad[ic]=randrad[jc];
randrad[jc]=temp;
}
}
}
}
/* Routine to generate particle radii based on a sieve analysis */
/* nsort is the number of particle radii to generate */
void gensieve(nsort)
long int nsort;
{
FILE *fsieve; /* FILE identifier for sieve file */
float dlow,dhigh; /* diameters for an individual sieve class */
float vlow,vhigh; /* volumes for an individual sieve class */
/* number fraction of particles for an individual sieve class */
float npart,
step, /* step size in volume space for a sieve */
vone, /* volume of an individual particle */
rcubed; /* radius cubed for an individual particle */
double vtot,stot; /* total volume and surface areas for system */
int nsieve,isv; /* number of sieve classes and loop index */
long int ngen, /* counter for number of particles generated */
itodo, /* number of particles needed */
ido; /* loop index for number of particles needed */
ngen=0; /* Counter for number of particles generated */
vtot=stot=0.0; /* Sums for volume and surface area counts */
fsieve=fopen("sieve.dat","r");
/* Read in the number of sieves being considered */
fscanf(fsieve,"%d\n",&nsieve);
/* loop for each sieve */
for(isv=1;isv<=nsieve;isv++){
/* Read in bounding diameters and number fraction */
/* of particles for this sieve */
fscanf(fsieve,"%f %f %f\n",&dlow,&dhigh,&npart);
/* Compute particle volumes and distribute requested number */
/* of particles uniformly in volume space */
vlow=PIVAL*CUB(dlow)/6.;
vhigh=PIVAL*CUB(dhigh)/6.;
printf("Seive %d - Vlow= %f Vhigh= %f frac= %f \n",isv,vlow,vhigh,npart);
itodo=(int)((float)nsort*npart+0.5); /* No. of particles needed */
/* For last sieve, generate all remaining needed particles */
/* to handle roundoff error, etc. */
if(isv==nsieve){itodo=nsort-ngen;}
step=0.0;
/* Determine volume step for sequential particles */
if(itodo>1){step=(vhigh-vlow)/(float)(itodo-1);}
printf("For particle size %f no. of particles is %ld \n",vlow,itodo);
fflush(stdout);
/* loop for all needed particles */
for(ido=1;ido<=itodo;ido++){
ngen+=1;
/* Compute volume of this particle */
vone=(vhigh-step*(float)(ido-1));
if(itodo==1){vone=((vhigh+vlow)/2.);}
/* Convert to radius */
rcubed=3.*vone/(4.*PIVAL);
/* Convert to system coordinates */
randrad[ngen]=(SYSIZE/REALSIZE)*pow(rcubed,(1./3.));
/* Accumulate volume and surface area counts */
vtot+=4.*PIVAL*CUB(randrad[ngen])/3.;
stot+=4.*PIVAL*SQR(randrad[ngen]);
}
}
fclose(fsieve);
/* Output volume and surface area in total */
printf("Volume and surface area are %lf and %lf \n",vtot,stot);
fflush(stdout);
}