Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 0a1513c6634a
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 4a1eeff37e08
Choose a head ref
  • 3 commits
  • 1 file changed
  • 1 contributor

Commits on Jun 17, 2019

  1. Copy the full SHA
    c2ca6f2 View commit details

Commits on Jun 18, 2019

  1. After egg's review.

    pleroy committed Jun 18, 2019
    Copy the full SHA
    9457015 View commit details
  2. Merge pull request #2211 from pleroy/FunctionRewrite2

     Rewrite BulirschCel and FukushimaEllipticBDJ
    pleroy authored Jun 18, 2019

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    4a1eeff View commit details
Showing with 56 additions and 47 deletions.
  1. +56 −47 numerics/elliptic_integrals.cpp
103 changes: 56 additions & 47 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#include "glog/logging.h"

#include <limits>

#include "numerics/elliptic_integrals.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
@@ -14,6 +16,8 @@
// 3. Figure something for the uninitialized variables.

// Bibliography:
// [Buli69] Bulirsch (1969), Numerical Calculation of Elliptic Integrals and
// Elliptic Fuctions. III.
// [Fuku11a] Fukushima (2011), Precise and fast computation of the general
// complete elliptic integral of the second kind.
// [Fuku11b] Fukushima (2011), Precise and fast computation of a general
@@ -27,6 +31,7 @@

namespace principia {

using quantities::Abs;
using quantities::Angle;
using quantities::Cos;
using quantities::Sin;
@@ -37,8 +42,8 @@ namespace numerics {

namespace {

// Bulirsch's cel function, [NIST10], 19.2(iii).
double BulirschCel(double kc0, double nc, double aa, double bb, int& err);
// Bulirsch's cel function, [Buli69], [NIST10], 19.2(iii).
double BulirschCel(double kc, double nc, double a, double b);

// Fukushima's complete elliptic integrals of the second kind [Fuku11a].
void FukushimaEllipticBD(double mc, double& elb, double& eld);
@@ -93,66 +98,70 @@ double FukushimaT(double t, double h);
// "Numerical computation of elliptic integrals and elliptic functions
// III"
//
// Inputs: kc0 = complimentary modulus 0 <= kc0 <= 1
// nc = complimentary characteristic 0 <= nc <= 1
// aa, bb = coefficients
// Inputs: kc = complementary modulus 0 <= kc <= 1
// nc = complementary characteristic 0 <= nc <= 1
// a, b = coefficients
//
// Outputs: cel = integral value
// err = integer error indicator
// Outputs: integral value
//
double BulirschCel(double const kc0,
double BulirschCel(double kc,
double const nc,
double const aa,
double const bb,
int& err) {
double out, a, b = 0, e, f, g; // TODO(phl): Initial value?
double p, q, kc, em, qc;
err = 0;
kc = kc0;
double a,
double b) {
// These values should give us 14 digits of accuracy, see [Buli69].
constexpr double ca = 1.0e-7;
constexpr double kc_nearly_0 = 1.0e-14;

// The identifiers below follow exactly [Buli69]. Note the (uncommon) use of
// non-const parameters to mimic [Buli69].
double p = nc;
if (kc == 0.0) {
if (b == 0.0) {
kc = 1.e-16;
kc = kc_nearly_0;
} else {
err = 1;
return 1.e99;
// "If in this case b ≠ 0 then cel is undefined."
DLOG(ERROR) << "kc = " << kc << " nc = " << nc << " a = " << a
<< " b = " << b;
return std::numeric_limits<double>::quiet_NaN();
}
}
qc = abs(kc);
a = aa;
b = bb;
p = nc;
e = qc;
em = 1.0;
kc = Abs(kc);
double e = kc;
double m = 1.0;

// Initial values for p, a, b.
if (p > 0.0) {
p = sqrt(p);
p = Sqrt(p);
b = b / p;
} else {
f = qc * qc;
q = 1.0 - f;
g = 1.0 - p;
double f = kc * kc;
double q = 1.0 - f;
double g = 1.0 - p;
f = f - p;
q = q * (b - a * p);
p = sqrt(f / g);
q = (b - a * p) * q;
p = Sqrt(f / g);
a = (a - b) / g;
b = a * p - q / (g * g * p);
}

// Bartky's algorithm.
for (;;) {
f = a;
a = a + b / p;
g = e / p;
b = b + f * g;
b = 2.0 * b;
p = p + g;
g = em;
em = em + qc;
if (abs(g - qc) < 1.e-7 * qc) {
double f = a;
a = b / p + a;
double g = e / p;
b = f * g + b;
b = b + b;
p = g + p;
g = m;
m = kc + m;
if (Abs(g - kc) <= g * ca) {
break;
}
qc = 2.0 * sqrt(e);
e = qc * em;
kc = sqrt(e);
kc = kc + kc;
e = kc * m;
}
out = 1.5707963267948966 * (b + a * em) / (em * (em + p));
return out;
return (π / 2) * (a * m + b) / (m * (m + p));
}

// Double precision general complete elliptic integrals of the second kind
@@ -617,11 +626,11 @@ void FukushimaEllipticBDJ(double const nc,
double& bc,
double& dc,
double& jc) {
double kc;
int err;
FukushimaEllipticBD(mc, bc, dc);
kc = sqrt(mc);
jc = BulirschCel(kc, nc, 0.0, 1.0, err);

// See [Buli69], special examples after equation (1.2.2).
double const kc = Sqrt(mc);
jc = BulirschCel(kc, nc, /*a=*/0.0, /*b=*/1.0);
}

void FukushimaEllipticBcDcJc(double const c₀,