Untitled

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