Untitled

void RadialDistrib2D (int numPart, 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 rho = (double)numPart / (double)(max_x*max_z);

  double *com_x = (double *) malloc (numPart*sizeof(double));
  double *com_y = (double *) malloc (numPart*sizeof(double));
  double *com_z = (double *) malloc (numPart*sizeof(double));
  double *bin      = calloc (numBins, sizeof(double));
  double *bin_tAvg = calloc (numBins, sizeof(double));

  for (int step = t_init; step <= t_end; step += interval) {
    sprintf (filename, "%s/data/COM_%d.dat", work_dir, step);
    stream = fopen (filename, "r");
    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 < numBins; n++) { // zero variables at every time frame
      bin[n]=0;
    }
    for (int n=0; n < numPart; n++) {
      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;
            bin_tAvg[m] = bin_tAvg[m]+1; // bins for computing time average 
            break;
          }
        }
      }
    }
    for (int m=0; m < numBins; m++) { // average over all particles
      bin[m] /= numPart;
      bin_tAvg[m] /= numPart;

    sprintf (filename, "%s/data/RDF_2D_%d.dat", work_dir, step);
    stream = fopen (filename, "w");
    for (int n=0; n < numBins; n++) {
      double r = (n+1)*binWidth;
      double tick = r + 0.5*binWidth;
      double dr = binWidth;
      double g_r = bin[n] / (2*PI*r*dr) / rho; // local density normalized by the global density
      fprintf (stream, "%lf %lf\n", tick, g_r);
    }
    fclose (stream);
  }

  sprintf (filename, "%s/data/RDF_2D_tAvg.dat", work_dir);
  stream = fopen (filename, "w");
  for (int n=0; n < numBins; n++) {
    double r = (n+1)*binWidth;
    double tick = r + 0.5*binWidth;
    double dr = binWidth;
    bin_tAvg[n] /= numFrames; // average over all time frames
    double g_r = bin_tAvg[n] / (2*PI*r*dr) / rho; // local density normalized by the global density
    fprintf (stream, "%lf %lf\n", tick, g_r);
  }
  fclose (stream);

  double nearNeighbors_r12 = 0;
  double nearNeighbors_r32 = 0;
  for (int n=0; n < numBins; n++) {
    double r = (n+1)*binWidth;
    if (r <= 12.0) {
      nearNeighbors_r12 += bin_tAvg[n];
    }
    if (r <= 32.0) {
      nearNeighbors_r32 += bin_tAvg[n];
    }
  }
  sprintf (filename, "%s/data/nearNeighbors_2D_tAvg.dat", work_dir);
  stream = fopen (filename, "w");
  fprintf (stream, "%lf %lf\n", nearNeighbors_r12, nearNeighbors_r32);
  fclose (stream);

  free (com_x);
  free (com_y);
  free (com_z);
  free (bin);
  free (bin_tAvg);
}