Use smallest encircling sphere to estimate particle radius in...

  if (_use_external_sphere) {

    result = _radii_of_spheres[0] * _rcf[0][_nshl_vec[0] - 1];

  } else {

    double x, y, z;

    int far_index = -1;

    double dist2, max_dist;

    double max_dist2 = 0.0;

    double avgX = 0.0, avgY = 0.0, avgZ = 0.0;

#pragma omp parallel for reduction(+: avgX, avgY, avgZ)

    for (int si = 0; si < _number_of_spheres; si++) {

      x = gc->get_sph_x(si);

      y = gc->get_sph_y(si);

      z = gc->get_sph_z(si);

      dist2 = x * x + y * y + z * z;

      avgX += gc->get_sph_x(si);

      avgY += gc->get_sph_y(si);

      avgZ += gc->get_sph_z(si);

    }

    avgX /= _number_of_spheres;

    avgY /= _number_of_spheres;

    avgZ /= _number_of_spheres;

    int configuration = 0;

    for (int si = 0; si < _number_of_spheres; si++) {

      int iogi = _iog_vec[si];

      if (iogi >= si + 1) configuration++;

      double dX = gc->get_sph_x(si) - avgX;

      double dY = gc->get_sph_y(si) - avgY;

      double dZ = gc->get_sph_z(si) - avgZ;

      double radius = _radii_of_spheres[configuration - 1];

      dist2 = dX * dX + dY * dY + dZ * dZ + radius * radius;

      if (dist2 > max_dist2) {

	max_dist2 = dist2;

	far_index = si;

      }

    }

    max_dist = sqrt(max_dist2);

    result = max_dist + _radii_of_spheres[far_index];

    result = max_dist;

  }

  return result;

}