/*
   ###########################################################################
   This function computes the circumference or perimeter of an ellipse with
   the given the parameters (r, R).  It is a numerical solution to the
   complete elliptic integral of the second kind.

   Arbitrary-precision arithmetic is used internally for high precision.

   ARGUMENTS:
   Rp = r = Polar radius = A numeric string, like '12.345'
   Re = R = Equatorial radius = A numeric string, like '12.345'

   Internally, the larger of the two arguments will automatically
   be interpreted as the equatorial radius (R), where (R >= r).

   If (r == R), then form = circle, otherwise, form = ellipse.

   ERRORS:
   FALSE is returned if either argument is non-numeric.

   -------
   CAVEAT:
   This function may take a long time to converge for large elliptic ratios
   (e.g. R:r > 4.5). It is meant for practical problems involving Earth,
   lunar and planetary spheroid meridian and elliptical orbit lengths,
   for which it should quickly converge on the numerical solution to
   16 decimals in a reasonable time.  If the R:r ratio is too large,
   it could result in a PHP script time-out, depending on the PHP
   INI settings.

   NO DEPENDENCIES
   ###########################################################################
*/

   function BC_Ellipse_Circum ($Rp='1.0', $Re='1.0')
{
/* -------------------------------------------------
   Read input (r, R) ellipsoid radii arguments. They
   must be purely numeric strings where: (r <= R).
   If (r > R), then their values will be swapped.
*/
   $r = trim($Rp);
   $R = trim($Re);

// ---------------------------------------------
// Return FALSE if either (r, R) is non-numeric.

   if (!Is_Numeric($r) or !Is_Numeric($R))  {return FALSE;}

// --------------------------------------
// Define constant (2*pi) to 64 decimals.

   $Twopi = '6.2831853071795864769252867665590057683943387987502116419498891846';

// ------------------------------------
// Set to 16 decimals output precision.

   $q = 16;

/* -----------------------------------------
   Use 32 decimals for internal computations
   to reduce accumulative rounding errors in
   the final 16-decimal result.
*/
   $Q = 32;

/* ----------------------------------------------
   If (r > R), then swap values so that (R >= r).
   R will always equate to the equatorial radius
   which is always the greater of the two values.
*/
   if (bcComp($R, $r, $Q) < 0)  {$w=$R;  $R=$r;  $r=$w;}

// ------------------------------------------------------
// Compute the square of the elliptic eccentricity value.

   $e2 = bcSub('1', bcDiv(bcMul($r,$r, $Q), bcMul($R,$R, $Q),$Q),$Q);

/* ---------------------------------------------
   Perform power series summation loop while the
   nth term is > the 16-decimal (PrecisionLimit).
   This limit could be changed, but 16 decimals
   should generally be sufficient.
*/
   $sum = 0;   $n = $epow2n = $PrecisionLevel = 1;

   $PrecisionLimit = '0.' . Str_Repeat('0', $Q-1) . '1';

   while (abs($PrecisionLevel) > $PrecisionLimit)
  {
   $N = 2*$n-1;  $X=1;  for($n=0;  $n <= ceil($N/2) - 1;  $n++)
  {$X = bcMul($X, $N - 2*$n);}  $X2 = bcMul($X,$X, $Q); // (2n-1)!!^2

   $N = 2*$n;    $Y=1;  for($n=0;  $n <= ceil($N/2) - 1;  $n++)
  {$Y = bcMul($Y, $N - 2*$n);} $Y2 = bcMul($Y,$Y,$Q); // (2n)!!^2

   $epow2n      = bcMul($epow2n, $e2, $Q); // e^2n
   $E           = bcMul($X2, $epow2n, $Q); // (2n-1)!!^2 * e^2n
   $F           = bcMul(2*$n-1, $Y2,  $Q); // (2n)!!^2 * 2n
   $CurrSumTerm = bcDiv($E, $F, $Q);
   $sum         = bcAdd($sum, $CurrSumTerm, $Q);
   $n++;
   $PrecisionLevel = bcSub($CurrSumTerm, $PrecisionLimit, $Q);
  }
   $CurrSumTerm = bcDiv(bcMul($X2, $epow2n, $Q), bcMul(2*$n-1, $Y2, $Q),$Q);

/* ----------------------------------------------------
   Finish computing the full elliptic circumference and
   round-off result to the given number of decimals (16).
*/
   $w = bcMul($Twopi, bcMul($R, bcSub('1.0', $sum, $Q),$Q),$Q);

   return bcAdd($w, '0.'.Str_Repeat('0', $q).'5', $q);

} // End of  BC_Ellipse_Circum (...)