/*
* Math.js
* Transcription of Math.hpp, Constants.hpp, and Accumulator.hpp into
* JavaScript.
*
* Copyright (c) Charles Karney (2011-2021) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*/
/**
* @namespace geodesic
* @description The parent namespace for the following modules:
* - {@link module:geodesic/Geodesic geodesic/Geodesic} The main
* engine for solving geodesic problems via the
* {@link module:geodesic/Geodesic.Geodesic Geodesic} class.
* - {@link module:geodesic/GeodesicLine geodesic/GeodesicLine}
* computes points along a single geodesic line via the
* {@link module:geodesic/GeodesicLine.GeodesicLine GeodesicLine}
* class.
* - {@link module:geodesic/PolygonArea geodesic/PolygonArea}
* computes the area of a geodesic polygon via the
* {@link module:geodesic/PolygonArea.PolygonArea PolygonArea}
* class.
* - {@link module:geodesic/Constants geodesic/Constants} defines
* constants specifying the version numbers and the parameters for the WGS84
* ellipsoid.
*
* The following modules are used internally by the package:
* - {@link module:geodesic/Math geodesic/Math} defines various
* mathematical functions.
* - {@link module:geodesic/Accumulator geodesic/Accumulator}
* interally used by
* {@link module:geodesic/PolygonArea.PolygonArea PolygonArea} (via the
* {@link module:geodesic/Accumulator.Accumulator Accumulator} class)
* for summing the contributions to the area of a polygon.
*/
// To allow swap via [y, x] = [x, y]
/* jshint esversion: 6 */
var geodesic = {};
geodesic.Constants = {};
geodesic.Math = {};
geodesic.Accumulator = {};
(function(
/**
* @exports geodesic/Constants
* @description Define constants defining the version and WGS84 parameters.
*/
c) {
"use strict";
/**
* @constant
* @summary WGS84 parameters.
* @property {number} a the equatorial radius (meters).
* @property {number} f the flattening.
*/
c.WGS84 = { a: 6378137, f: 1/298.257223563 };
/**
* @constant
* @summary an array of version numbers.
* @property {number} major the major version number.
* @property {number} minor the minor version number.
* @property {number} patch the patch number.
*/
c.version = { major: 2,
minor: 0,
patch: 0 };
/**
* @constant
* @summary version string
*/
c.version_string = "2.0.0";
})(geodesic.Constants);
(function(
/**
* @exports geodesic/Math
* @description Some useful mathematical constants and functions (mainly for
* internal use).
*/
m) {
"use strict";
/**
* @summary The number of digits of precision in floating-point numbers.
* @constant {number}
*/
m.digits = 53;
/**
* @summary The machine epsilon.
* @constant {number}
*/
m.epsilon = Math.pow(0.5, m.digits - 1);
/**
* @summary The factor to convert degrees to radians.
* @constant {number}
*/
m.degree = Math.PI/180;
/**
* @summary Square a number.
* @param {number} x the number.
* @return {number} the square.
*/
m.sq = function(x) { return x * x; };
/**
* @summary The hypotenuse function.
* @param {number} x the first side.
* @param {number} y the second side.
* @return {number} the hypotenuse.
*/
m.hypot = function(x, y) {
// Built in Math.hypot give incorrect results from GeodSolve92.
return Math.sqrt(x*x + y*y);
};
/**
* @summary Cube root function.
* @param {number} x the argument.
* @return {number} the real cube root.
*/
m.cbrt = Math.cbrt || function(x) {
var y = Math.pow(Math.abs(x), 1/3);
return x > 0 ? y : (x < 0 ? -y : x);
};
/**
* @summary The log1p function.
* @param {number} x the argument.
* @return {number} log(1 + x).
*/
m.log1p = Math.log1p || function(x) {
var y = 1 + x,
z = y - 1;
// Here's the explanation for this magic: y = 1 + z, exactly, and z
// approx x, thus log(y)/z (which is nearly constant near z = 0) returns
// a good approximation to the true log(1 + x)/x. The multiplication x *
// (log(y)/z) introduces little additional error.
return z === 0 ? x : x * Math.log(y) / z;
};
/**
* @summary Inverse hyperbolic tangent.
* @param {number} x the argument.
* @return {number} tanh<sup>−1</sup> x.
*/
m.atanh = Math.atanh || function(x) {
var y = Math.abs(x); // Enforce odd parity
y = m.log1p(2 * y/(1 - y))/2;
return x > 0 ? y : (x < 0 ? -y : x);
};
/**
* @summary Copy the sign.
* @param {number} x gives the magitude of the result.
* @param {number} y gives the sign of the result.
* @return {number} value with the magnitude of x and with the sign of y.
*/
m.copysign = function(x, y) {
return Math.abs(x) * (y < 0 || (y === 0 && 1/y < 0) ? -1 : 1);
};
/**
* @summary An error-free sum.
* @param {number} u
* @param {number} v
* @return {object} sum with sum.s = round(u + v) and sum.t is u + v −
* round(u + v)
*/
m.sum = function(u, v) {
var s = u + v,
up = s - v,
vpp = s - up,
t;
up -= u;
vpp -= v;
// if s = 0, then t = 0 and give t the same sign as s
t = s ? 0 - (up + vpp) : s;
// u + v = s + t
// = round(u + v) + t
return {s: s, t: t};
};
/**
* @summary Evaluate a polynomial.
* @param {integer} N the order of the polynomial.
* @param {array} p the coefficient array (of size N + 1) (leading
* order coefficient first)
* @param {number} x the variable.
* @return {number} the value of the polynomial.
*/
m.polyval = function(N, p, s, x) {
var y = N < 0 ? 0 : p[s++];
while (--N >= 0) y = y * x + p[s++];
return y;
};
/**
* @summary Coarsen a value close to zero.
* @param {number} x
* @return {number} the coarsened value.
*/
m.AngRound = function(x) {
// The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 for
// reals = 0.7 pm on the earth if x is an angle in degrees. (This is about
// 1000 times more resolution than we get with angles around 90 degrees.)
// We use this to avoid having to deal with near singular cases when x is
// non-zero but tiny (e.g., 1.0e-200). This converts -0 to +0; however
// tiny negative numbers get converted to -0.
var z = 1/16,
y = Math.abs(x);
// The compiler mustn't "simplify" z - (z - y) to y
y = y < z ? z - (z - y) : y;
return m.copysign(y, x);
};
/**
* @summary The remainder function.
* @param {number} x the numerator of the division
* @param {number} y the denominator of the division
* @return {number} the remainder in the range [−y/2, y/2].
* <p>
* The range of x is unrestricted; y must be positive.
*/
m.remainder = function(x, y) {
x %= y;
return x < -y/2 ? x + y : (x < y/2 ? x : x - y);
};
/**
* @summary Normalize an angle.
* @param {number} x the angle in degrees.
* @return {number} the angle reduced to the range [−180°,
* 180°].
*
* The range of x is unrestricted. If the result is ±0° or
* ±180° then the sign is the sign of \e x.
*/
m.AngNormalize = function(x) {
// Place angle in [-180, 180].
var y = m.remainder(x, 360);
return Math.abs(y) === 180 ? m.copysign(180, x) : y;
};
/**
* @summary Normalize a latitude.
* @param {number} x the angle in degrees.
* @return {number} x if it is in the range [−90°, 90°],
* otherwise return NaN.
*/
m.LatFix = function(x) {
// Replace angle with NaN if outside [-90, 90].
return Math.abs(x) > 90 ? NaN : x;
};
/**
* @summary The exact difference of two angles reduced to [−180°,
* 180°]
* @param {number} x the first angle in degrees.
* @param {number} y the second angle in degrees.
* @return {object} diff the exact difference, diff.d + diff.e = y − x
* mod 360°.
*
* This computes z = y − x exactly, reduced to
* [−180°, 180°]; and then sets z = d + e where d
* is the nearest representable number to z and e is the truncation
* error. If z = ±0° or ±180°, then the sign of
* d is given by the sign of y − x. The maximum absolute
* value of e is 2<sup>−26</sup> (for doubles).
*/
m.AngDiff = function(x, y) {
// Compute y - x and reduce to [-180,180] accurately.
var r = m.sum(m.remainder(-x, 360), m.remainder(y, 360)), d, e;
r = m.sum(m.remainder(r.s, 360), r.t);
d = r.s;
e = r.t;
// Fix the sign if d = -180, 0, 180.
if (d === 0 || Math.abs(d) === 180)
// If e == 0, take sign from y - x
// else (e != 0, implies d = +/-180), d and e must have opposite signs
d = m.copysign(d, e === 0 ? y - x : -e);
return {d: d, e: e};
};
/**
* @summary Evaluate the sine and cosine function with the argument in
* degrees
* @param {number} x in degrees.
* @return {object} r with r.s = sin(x) and r.c = cos(x).
*/
m.sincosd = function(x) {
// In order to minimize round-off errors, this function exactly reduces
// the argument to the range [-45, 45] before converting it to radians.
var r, q, s, c, sinx, cosx;
r = x % 360;
q = Math.round(r / 90); // If r is NaN this returns NaN
r -= 90 * q;
// now abs(r) <= 45
r *= this.degree;
// Possibly could call the gnu extension sincos
s = Math.sin(r); c = Math.cos(r);
switch (q & 3) {
case 0: sinx = s; cosx = c; break;
case 1: sinx = c; cosx = -s; break;
case 2: sinx = -s; cosx = -c; break;
default: sinx = -c; cosx = s; break; // case 3
}
cosx += 0;
if (sinx === 0) sinx = m.copysign(sinx, x);
return {s: sinx, c: cosx};
};
/**
* @summary Evaluate the sine and cosine with reduced argument plus correction
*
* @param {number} x reduced angle in degrees.
* @param {number} t correction in degrees.
* @return {object} r with r.s = sin(x + t) and r.c = cos(x + t).
*/
m.sincosde = function(x, t) {
// In order to minimize round-off errors, this function exactly reduces
// the argument to the range [-45, 45] before converting it to radians.
var r, q, s, c, sinx, cosx;
r = x % 360;
q = Math.round(r / 90); // If r is NaN this returns NaN
r = m.AngRound((r - 90 * q) + t);
// now abs(r) <= 45
r *= this.degree;
// Possibly could call the gnu extension sincos
s = Math.sin(r); c = Math.cos(r);
switch (q & 3) {
case 0: sinx = s; cosx = c; break;
case 1: sinx = c; cosx = -s; break;
case 2: sinx = -s; cosx = -c; break;
default: sinx = -c; cosx = s; break; // case 3
}
cosx += 0;
if (sinx === 0) sinx = m.copysign(sinx, x);
return {s: sinx, c: cosx};
};
/**
* @summary Evaluate the atan2 function with the result in degrees
* @param {number} y
* @param {number} x
* @return atan2(y, x) in degrees, in the range [−180°
* 180°].
*/
m.atan2d = function(y, x) {
// In order to minimize round-off errors, this function rearranges the
// arguments so that result of atan2 is in the range [-pi/4, pi/4] before
// converting it to degrees and mapping the result to the correct
// quadrant.
var q = 0, ang;
if (Math.abs(y) > Math.abs(x)) { [y, x] = [x, y]; q = 2; } // swap(x, y)
if (m.copysign(1, x) < 0) { x = -x; ++q; }
// here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
ang = Math.atan2(y, x) / this.degree;
switch (q) {
// Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
// atan2d will not be called with y = -0. If need be, include
//
// case 0: ang = 0 + ang; break;
//
// and handle mpfr as in AngRound.
case 1: ang = m.copysign(180, y) - ang; break;
case 2: ang = 90 - ang; break;
case 3: ang = -90 + ang; break;
default: break;
}
return ang;
};
})(geodesic.Math);
(function(
/**
* @exports geodesic/Accumulator
* @description Accurate summation via the
* {@link module:geodesic/Accumulator.Accumulator Accumulator} class
* (mainly for internal use).
*/
a, m) {
"use strict";
/**
* @class
* @summary Accurate summation of many numbers.
* @classdesc This allows many numbers to be added together with twice the
* normal precision. In the documentation of the member functions, sum
* stands for the value currently held in the accumulator.
* @param {number | Accumulator} [y = 0] set sum = y.
*/
a.Accumulator = function(y) {
this.Set(y);
};
/**
* @summary Set the accumulator to a number.
* @param {number | Accumulator} [y = 0] set sum = y.
*/
a.Accumulator.prototype.Set = function(y) {
if (!y) y = 0;
if (y.constructor === a.Accumulator) {
this._s = y._s;
this._t = y._t;
} else {
this._s = y;
this._t = 0;
}
};
/**
* @summary Add a number to the accumulator.
* @param {number} [y = 0] set sum += y.
*/
a.Accumulator.prototype.Add = function(y) {
// Here's Shewchuk's solution...
// Accumulate starting at least significant end
var u = m.sum(y, this._t),
v = m.sum(u.s, this._s);
u = u.t;
this._s = v.s;
this._t = v.t;
// Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u)
// exactly with s, t, u non-adjacent and in decreasing order (except
// for possible zeros). The following code tries to normalize the
// result. Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
// _s). The follow does an approximate job (and maintains the
// decreasing non-adjacent property). Here are two "failures" using
// 3-bit floats:
//
// Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
// [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
//
// Case 2: _s+_t is not as close to s+t+u as it shold be
// [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
// should be [80, -7] = 73 (exact)
//
// "Fixing" these problems is probably not worth the expense. The
// representation inevitably leads to small errors in the accumulated
// values. The additional errors illustrated here amount to 1 ulp of
// the less significant word during each addition to the Accumulator
// and an additional possible error of 1 ulp in the reported sum.
//
// Incidentally, the "ideal" representation described above is not
// canonical, because _s = round(_s + _t) may not be true. For
// example, with 3-bit floats:
//
// [128, 16] + 1 -> [160, -16] -- 160 = round(145).
// But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
//
if (this._s === 0) // This implies t == 0,
this._s = u; // so result is u
else
this._t += u; // otherwise just accumulate u to t.
};
/**
* @summary Return the result of adding a number to sum (but
* don't change sum).
* @param {number} [y = 0] the number to be added to the sum.
* @return sum + y.
*/
a.Accumulator.prototype.Sum = function(y) {
var b;
if (!y)
return this._s;
else {
b = new a.Accumulator(this);
b.Add(y);
return b._s;
}
};
/**
* @summary Set sum = −sum.
*/
a.Accumulator.prototype.Negate = function() {
this._s *= -1;
this._t *= -1;
};
/**
* @summary Take the remainder
* @param {number} y the divisor of the remainder operation.
* @return sum in range [−y/2, y/2].
*/
a.Accumulator.prototype.Remainder = function(y) {
this._s = m.remainder(this._s, y);
this.Add(0);
};
})(geodesic.Accumulator, geodesic.Math);