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: 9db980914fa1
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: a0db3bdfbd1b
Choose a head ref
  • 5 commits
  • 1 file changed
  • 1 contributor

Commits on Jun 16, 2019

  1. Copy the full SHA
    0885145 View commit details
  2. Use our Sqrt.

    pleroy committed Jun 16, 2019
    Copy the full SHA
    aeecef2 View commit details
  3. Copy the full SHA
    9101a53 View commit details
  4. Typo and a bit of reordering.

    pleroy committed Jun 16, 2019
    Copy the full SHA
    602cf98 View commit details

Commits on Jun 17, 2019

  1. Merge pull request #2209 from pleroy/Identifiers2

    Better identifiers and control structure for FukushimaEllipticBsDsJs
    pleroy authored Jun 17, 2019
    Copy the full SHA
    a0db3bd View commit details
Showing with 40 additions and 40 deletions.
  1. +40 −40 numerics/elliptic_integrals.cpp
80 changes: 40 additions & 40 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -674,52 +674,52 @@ void FukushimaEllipticBcDcJc(double const c0,
}
}

void FukushimaEllipticBsDsJs(double const s0,
void FukushimaEllipticBsDsJs(double const s₀,
double const n,
double const mc,
double& b,
double& d,
double& j) {
double m, h, del, s, y, c, sy, t;
double yy[12];
double ss[12];
double cd[12];

m = 1.0 - mc;
h = n * (1.0 - n) * (n - m);
del = 0.01622;
s = s0;
y = s * s;
if (y < del) {
FukushimaEllipticBsDsMaclaurinSeries(y, m, b, d);
b = s * b;
d = s * y * d;
j = s * FukushimaEllipticJsMaclaurinSeries(y, n, m);
return;
}
yy[1] = y;
ss[1] = s;
int i;
for (i = 1; i <= 10; ++i) {
c = sqrt(1.0 - y);
d = sqrt(1.0 - m * y);
y = y / ((1.0 + c) * (1.0 + d));
yy[i + 1] = y;
ss[i + 1] = sqrt(y);
cd[i] = c * d;
if (y < del) {
break;
}
LOG_IF(FATAL, i == 10) << "(elsbdj) too many iterations: s0,m=" << s0 << " "
<< m;
// See [Fuku11c] section 3.5 for the determination of yB.
constexpr double yB = 0.01622;
// The maximum number of argument transformations, related to yB. This is the
// maximum number of iterations in the first loop below.
constexpr int max_transformations = 10;

double y[max_transformations + 1];
double s[max_transformations + 1];
double cd[max_transformations + 1];

// Half and double argument transformations, [Fuku11c] section 3.3.
double const m = 1.0 - mc;
double const h = n * (1.0 - n) * (n - m);
double const y₀ = s₀ * s₀;

// Half argument transformation of s.
y[0] = y₀;
s[0] = s₀;
double yᵢ = y₀;
int i = 0; // Note that this variable is used after the loop.
for (; yᵢ >= yB; ++i) {
DCHECK_LT(i, max_transformations);
double const cᵢ = Sqrt(1.0 - yᵢ);
double const dᵢ = Sqrt(1.0 - m * yᵢ);
yᵢ = yᵢ / ((1.0 + cᵢ) * (1.0 + dᵢ));
y[i + 1] = yᵢ;
s[i + 1] = Sqrt(yᵢ);
cd[i] = cᵢ * dᵢ;
}
FukushimaEllipticBsDsMaclaurinSeries(y, m, b, d);
b = ss[i + 1] * b;
d = ss[i + 1] * y * d;
j = ss[i + 1] * FukushimaEllipticJsMaclaurinSeries(y, n, m);
for (int k = i; k >= 1; --k) {
sy = ss[k] * yy[k + 1];
t = sy / (1.0 - n * (yy[k] - yy[k + 1] * cd[k]));

// Maclaurin series.
FukushimaEllipticBsDsMaclaurinSeries(yᵢ, m, b, d);
b = s[i] * b;
d = s[i] * yᵢ * d;
j = s[i] * FukushimaEllipticJsMaclaurinSeries(yᵢ, n, m);

// Double argument transformation of J.
for (int k = i; k > 0; --k) {
double const sy = s[k - 1] * y[k];
double const t = sy / (1.0 - n * (y[k - 1] - y[k] * cd[k - 1]));
b = 2.0 * b - sy;
d = d + (d + sy);
j = j + (j + FukushimaT(t, h));