core: remove analyze_cell_gpb

Closes: #257

src/core/statistics.cpp

@@ -442,131 +442,6 @@ double distto(PartCfg &partCfg, double p[3], int pid) {
   return std::sqrt(mindist);
 }
 
-void calc_cell_gpb(double xi_m, double Rc, double ro, double gacc, int maxtry,
-                   double *result) {
-  double LOG, xi_min, RM, gamma, g1, g2, gmid = 0, dg, ig, f, fmid, rtb;
-  int i;
-  LOG = log(Rc / ro);
-  xi_min = LOG / (1 + LOG);
-  if (maxtry < 1)
-    maxtry = 1;
-
-  /* determine which of the regimes we are in: */
-  if (xi_m > 1) {
-    ig = 1.0;
-    g1 = PI / LOG;
-    g2 = PI / (LOG + xi_m / (xi_m - 1.0));
-  } else if (xi_m == 1) {
-    ig = 1.0;
-    g1 = (PI / 2.0) / LOG;
-    g2 = (PI / 2.0) / (LOG + 1.0);
-  } else if (xi_m == xi_min) {
-    ig = 1.0;
-    g1 = g2 = 0.0;
-  } else if (xi_m > xi_min) {
-    ig = 1.0;
-    g1 = (PI / 2.0) / LOG;
-    g2 =
-        sqrt(3.0 * (LOG - xi_m / (1.0 - xi_m)) / (1 - pow((1.0 - xi_m), -3.0)));
-  } else if (xi_m > 0.0) {
-    ig = -1.0;
-    g1 = 1 - xi_m;
-    g2 = xi_m * (6.0 - (3.0 - xi_m) * xi_m) / (3.0 * LOG);
-  } else if (xi_m == 0.0) {
-    ig = -1.0;
-    g1 = g2 = 1 - xi_m;
-  } else {
-    result[2] = -5.0;
-    return;
-  }
-
-  /* decide which method to use (if any): */
-  if (xi_m == xi_min) {
-    gamma = 0.0;
-    RM = 0.0;
-  } else if (xi_m == 0.0) {
-    gamma = 1 - xi_m;
-    RM = -1.0;
-  } else if (ig == 1.0) {
-    /* determine gamma via a bisection-search: */
-    f = atan(1.0 / g1) + atan((xi_m - 1.0) / g1) - g1 * LOG;
-    fmid = atan(1.0 / g2) + atan((xi_m - 1.0) / g2) - g2 * LOG;
-    if (f * fmid >= 0.0) {
-      /* failed to bracket function value with intial guess - abort: */
-      result[0] = f;
-      result[1] = fmid;
-      result[2] = -3.0;
-      return;
-    }
-
-    /* orient search such that the positive part of the function lies to the
-     * right of the zero */
-    rtb = f < 0.0 ? (dg = g2 - g1, g1) : (dg = g1 - g2, g2);
-    for (i = 1; i <= maxtry; i++) {
-      gmid = rtb + (dg *= 0.5);
-      fmid = atan(1.0 / gmid) + atan((xi_m - 1.0) / gmid) - gmid * LOG;
-      if (fmid <= 0.0)
-        rtb = gmid;
-      if (fabs(dg) < gacc || fmid == 0.0)
-        break;
-    }
-
-    if (fabs(dg) > gacc) {
-      /* too many iterations without success - abort: */
-      result[0] = gmid;
-      result[1] = dg;
-      result[2] = -2.0;
-      return;
-    }
-
-    /* So, these are the values for gamma and Manning-radius: */
-    gamma = gmid;
-    RM = Rc * exp(-(1.0 / gamma) * atan(1.0 / gamma));
-  } else if (ig == -1.0) {
-    /* determine -i*gamma: */
-    f = -1.0 * (atanh(g2) + atanh(g2 / (xi_m - 1))) - g2 * LOG;
-
-    /* modified orient search, this time starting from the upper bound
-     * backwards: */
-    if (f < 0.0) {
-      rtb = g1;
-      dg = g1 - g2;
-    } else {
-      fprintf(stderr, "WARNING: Lower boundary is actually larger than "
-                      "l.hpp.s, flipping!\n");
-      rtb = g1;
-      dg = g1;
-    }
-    for (i = 1; i <= maxtry; i++) {
-      gmid = rtb - (dg *= 0.5);
-      fmid = -1.0 * (atanh(gmid) + atanh(gmid / (xi_m - 1))) - gmid * LOG;
-      if (fmid >= 0.0)
-        rtb = gmid;
-      if (fabs(dg) < gacc || fmid == 0.0)
-        break;
-    }
-
-    if (fabs(dg) > gacc) {
-      /* too many iterations without success - abort: */
-      result[0] = gmid;
-      result[1] = dg;
-      result[2] = -2.0;
-      return;
-    }
-
-    /* So, these are the values for -i*gamma and Manning-radius: */
-    gamma = gmid;
-    RM = Rc * exp(atan(1.0 / gamma) / gamma);
-  } else {
-    result[2] = -5.0;
-    return;
-  }
-
-  result[0] = gamma;
-  result[1] = RM;
-  result[2] = ig;
-  return;
-}
 
 void calc_part_distribution(PartCfg &partCfg, int *p1_types, int n_p1,
                             int *p2_types, int n_p2, double r_min, double r_max,

src/core/statistics.hpp

@@ -162,32 +162,6 @@ void nbhood(PartCfg & partCfg, double pos[3], double r_catch, IntList *il, int p
    \<posy\>, \<posz\>). */
 double distto(PartCfg &, double pos[3], int pid);
 
-/** numerical solution for the integration constant \f$\gamma\f$ in the cell
-   model, determined by
-    \f[\gamma\,\ln\frac{R}{r_0}=\arctan\frac{1}{\gamma}+\arctan\frac{\xi_M-1}{\gamma}\f]
-    from which the second integration constant, the Manning radius \f$R_M\f$,
-   follows to
-    \f[R_M =
-   R\cdot\exp\left(-\frac{1}{\gamma}\cdot\arctan\frac{1}{\gamma}\right)\f]
-    Any value \f$\xi_M\geq 0\f$ is allowed, the function will automatically
-   ensure the
-    analytical continuation required for \f$\xi_M<\ln(R/r_0)/(1+\ln(R/r_0))\f$,
-   in which case
-    \f$\gamma\f$ becomes imaginary.
-    @param xi_m   Manning parameter \f$\xi_M=\ell_B/a\f$ (with Bjerrum-length
-   \f$\ell_B\f$ and charge distance \f$a\f$)
-    @param Rc     outer radius \f$R_C\f$ of the cylindrical cell around each
-   polyelectrolyte
-    @param ro     inner radius \f$r_0\f$ of the cylindrical cell around each
-   polyelectrolyte
-    @param gacc   the accuracy up to which \f$\gamma\f$ should be determined
-    @param maxtry maximum number of interations to find a solution
-    @param result pointer to double array containing \f$\gamma\f$ and \f$R_M\f$,
-                  and a third entry which is -1.0 if \f$\gamma\f$ is imaginary,
-   +1.0 else. */
-void calc_cell_gpb(double xi_m, double Rc, double ro, double gacc, int maxtry,
-                   double *result);
-
 /** appends particles' positions in 'partCfg' to onfigs */
 void analyze_append(PartCfg &);