void RadialDistrib2D (int t_init, int t_end, int interval, int max_x, int max_y, int max_z, char *work_dir) {
char filename[200];
FILE *stream;
double binWidth = 1.0;
int numBins = max_x/2 - 1; // g(r) calculation starts from r=1
int numFrames = 1 + (t_end - t_init) / interval;
double area = PI*(max_x/2)*(max_x/2);
double *bin = calloc (numBins, sizeof(double));
double *gr_avg = calloc (numBins, sizeof(double));
int numPHI=6;
int numCa=13;
double phi [6] = {2.3994, 2.9993, 3.4492, 3.8991, 4.9488, 5.9986};
int part[6] = {16, 20, 23, 26, 33, 40};
double ca [13] = {0.01, 0.02, 0.03, 0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2};
int trials[13];
for (int i=0; i < numPHI; i++) {
int numPart = part[i];
double *com_x = (double *) malloc (numPart*sizeof(double));
double *com_y = (double *) malloc (numPart*sizeof(double));
double *com_z = (double *) malloc (numPart*sizeof(double));
switch (i) {
case 0:
trials[0] =10;
trials[1] =10;
trials[2] =20;
trials[3] =20;
trials[4] =20;
trials[5] =20;
trials[6] =20;
trials[7] =20;
trials[8] =20;
trials[9] =10;
trials[10] =10;
trials[11] =20;
trials[12] =10;
break;
case 1:
trials[0] =10;
trials[1] =10;
trials[2] =20;
trials[3] =20;
trials[4] =20;
trials[5] =20;
trials[6] =10;
trials[7] =20;
trials[8] =20;
trials[9] =20;
trials[10] =10;
trials[11] =10;
trials[12] =10;
break;
case 2:
trials[0] =10;
trials[1] =10;
trials[2] =10;
trials[3] =20;
trials[4] =20;
trials[5] =10;
trials[6] =20;
trials[7] =10;
trials[8] =10;
trials[9] =10;
trials[10] =20;
trials[11] =20;
trials[12] =10;
break;
case 3:
trials[0] =10;
trials[1] =10;
trials[2] =10;
trials[3] =10;
trials[4] =20;
trials[5] =20;
trials[6] =20;
trials[7] =10;
trials[8] =10;
trials[9] =10;
trials[10] =10;
trials[11] =10;
trials[12] =10;
break;
case 4:
trials[0] =10;
trials[1] =10;
trials[2] =10;
trials[3] =10;
trials[4] =10;
trials[5] =10;
trials[6] =10;
trials[7] =10;
trials[8] =10;
trials[9] =10;
trials[10] =10;
trials[11] =10;
trials[12] =10;
break;
case 5:
trials[0] =10;
trials[1] =20;
trials[2] =10;
trials[3] =10;
trials[4] =10;
trials[5] =10;
trials[6] =10;
trials[7] =20;
trials[8] =10;
trials[9] =10;
trials[10] =10;
trials[11] =10;
trials[12] =10;
break;
}
for (int j=0; j < numCa; j++)
{
double nearNeighbors_r12 = 0; // zero variables for each condition
double nearNeighbors_r32 = 0;
double avgCountedParticles = 0;
for (int m=0; m < numBins; m++) {
gr_avg[m] = 0;
}
int numTrials = trials[j];
for (int k=1; k <= numTrials; k++) {
for (int step = t_init; step <= t_end; step += interval) {
sprintf (filename, "%s/h24_phi%.4lf_Re0.1_Ca%.2lf_WCA1_zero0.8-%d/data/COM_%d.dat", work_dir, phi[i], ca[j], k, step);
if (j==7 || j==12) {
sprintf (filename, "%s/h24_phi%.4lf_Re0.1_Ca%.1lf_WCA1_zero0.8-%d/data/COM_%d.dat", work_dir, phi[i], ca[j], k, step);
}
stream = fopen (filename, "r");
if (stream == NULL) {
printf("(phi, ca, index)=(%.4lf, %.2lf, %d)\n", phi[i], ca[j], k);
}
for (int n=0; n < numPart; n++) {
fscanf (stream, "%lf %lf %lf", &com_x[n], &com_y[n], &com_z[n]);
}
fclose (stream);
for (int n=0; n < numPart; n++)
{
int countedPart = 0; // zero containers (variables) for each trial
for(int m=0; m < numBins; m++) {
bin[m] = 0;
}
for (int k=0; k < numPart; k++) {
double dx = n_image((com_x[k]-com_x[n]), max_x); // take periodic conditions into account
double dz = n_image((com_z[k]-com_z[n]), max_z);
double dr = sqrt(dx*dx + dz*dz);
for (int m=0; m < numBins; m++) {
double lowerBound = (m+1)*binWidth;
double upperBound = lowerBound+binWidth;
if (dr >= lowerBound && dr < upperBound) {
bin[m] = bin[m]+1;
countedPart++;
avgCountedParticles++;
break;
}
}
}
for (int m=0; m < numBins; m++) {
double lowerBound = (m+1)*binWidth;
double upperBound = lowerBound + binWidth;
if (lowerBound <= 12.0) {
nearNeighbors_r12 += bin[m]; // # of neighbors within r <= 12
}
if (lowerBound <= 32.0) {
nearNeighbors_r32 += bin[m]; // # of neighbors within r <= 32
}
bin[m] = bin[m] / (PI*(upperBound*upperBound - lowerBound*lowerBound)) / (countedPart/area); // g(r) of a trial
gr_avg[m] += bin[m]; // average g(r)
}
}
}
}
for (int m=0; m < numBins; m++) {
gr_avg[m] /= (numPart*numFrames*numTrials);
}
nearNeighbors_r12 /= (numPart*numFrames*numTrials);
nearNeighbors_r32 /= (numPart*numFrames*numTrials);
avgCountedParticles /= (numPart*numFrames*numTrials);
sprintf (filename, "%s/RDF_2D_avg_phi%.4lf_Ca%.2lf.dat", work_dir, phi[i], ca[j]);
stream = fopen (filename, "w");
for (int n=0; n < numBins; n++) {
double r = (n+1)*binWidth;
double tick = r + 0.5*binWidth;
fprintf (stream, "%lf %lf\n", tick, gr_avg[n]);
}
fclose (stream);
sprintf (filename, "%s/nearNeighbors_2D_tAvg_phi%.4lf_Ca%.2lf.dat", work_dir, phi[i], ca[j]);
stream = fopen (filename, "w");
fprintf (stream, "%lf %lf %lf\n", nearNeighbors_r12, nearNeighbors_r32, avgCountedParticles);
fclose (stream);
}
free (com_x);
free (com_y);
free (com_z);
}
free (bin);
free (gr_avg);
}