/*
* GeodesicLine.js
* Transcription of GeodesicLine.[ch]pp into JavaScript.
*
* See the documentation for the C++ class. The conversion is a literal
* conversion from C++.
*
* The algorithms are derived in
*
* Charles F. F. Karney,
* Algorithms for geodesics, J. Geodesy 87, 43-55 (2013);
* https://doi.org/10.1007/s00190-012-0578-z
* Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
*
* Copyright (c) Charles Karney (2011-2022) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*/
// Load AFTER geodesic/Math.js, geodesic/Geodesic.js
(function(
g,
/**
* @exports geodesic/GeodesicLine
* @description Solve geodesic problems on a single geodesic line via the
* {@link module:geodesic/GeodesicLine.GeodesicLine GeodesicLine}
* class.
*/
l, m) {
"use strict";
/**
* @class
* @property {number} a the equatorial radius (meters).
* @property {number} f the flattening.
* @property {number} lat1 the initial latitude (degrees).
* @property {number} lon1 the initial longitude (degrees).
* @property {number} azi1 the initial azimuth (degrees).
* @property {number} salp1 the sine of the azimuth at the first point.
* @property {number} calp1 the cosine the azimuth at the first point.
* @property {number} s13 the distance to point 3 (meters).
* @property {number} a13 the arc length to point 3 (degrees).
* @property {bitmask} caps the capabilities of the object.
* @summary Initialize a GeodesicLine object. For details on the caps
* parameter, see {@tutorial 2-interface}, "The outmask and caps
* parameters".
* @classdesc Performs geodesic calculations along a given geodesic line.
* This object is usually instantiated by
* {@link module:geodesic/Geodesic.Geodesic#Line Geodesic.Line}.
* The methods
* {@link module:geodesic/Geodesic.Geodesic#DirectLine
* Geodesic.DirectLine} and
* {@link module:geodesic/Geodesic.Geodesic#InverseLine
* Geodesic.InverseLine} set in addition the position of a reference point
* 3.
* @param {object} geod a {@link module:geodesic/Geodesic.Geodesic
* Geodesic} object.
* @param {number} lat1 the latitude of the first point in degrees.
* @param {number} lon1 the longitude of the first point in degrees.
* @param {number} azi1 the azimuth at the first point in degrees.
* @param {bitmask} [caps = STANDARD | DISTANCE_IN] which capabilities to
* include; LATITUDE | AZIMUTH are always included.
*/
l.GeodesicLine = function(geod, lat1, lon1, azi1, caps, salp1, calp1) {
var t, cbet1, sbet1, eps, s, c;
if (!caps) caps = g.STANDARD | g.DISTANCE_IN;
this.a = geod.a;
this.f = geod.f;
this._b = geod._b;
this._c2 = geod._c2;
this._f1 = geod._f1;
this.caps = caps | g.LATITUDE | g.AZIMUTH | g.LONG_UNROLL;
this.lat1 = m.LatFix(lat1);
this.lon1 = lon1;
if (typeof salp1 === 'undefined' || typeof calp1 === 'undefined') {
this.azi1 = m.AngNormalize(azi1);
t = m.sincosd(m.AngRound(this.azi1)); this.salp1 = t.s; this.calp1 = t.c;
} else {
this.azi1 = azi1; this.salp1 = salp1; this.calp1 = calp1;
}
t = m.sincosd(m.AngRound(this.lat1)); sbet1 = this._f1 * t.s; cbet1 = t.c;
// norm(sbet1, cbet1);
t = m.hypot(sbet1, cbet1); sbet1 /= t; cbet1 /= t;
// Ensure cbet1 = +epsilon at poles
cbet1 = Math.max(g.tiny_, cbet1);
this._dn1 = Math.sqrt(1 + geod._ep2 * m.sq(sbet1));
// Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
this._salp0 = this.salp1 * cbet1; // alp0 in [0, pi/2 - |bet1|]
// Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
// is slightly better (consider the case salp1 = 0).
this._calp0 = m.hypot(this.calp1, this.salp1 * sbet1);
// Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
// sig = 0 is nearest northward crossing of equator.
// With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
// With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
// With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
// Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
// With alp0 in (0, pi/2], quadrants for sig and omg coincide.
// No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
// With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
this._ssig1 = sbet1; this._somg1 = this._salp0 * sbet1;
this._csig1 = this._comg1 =
sbet1 !== 0 || this.calp1 !== 0 ? cbet1 * this.calp1 : 1;
// norm(this._ssig1, this._csig1); // sig1 in (-pi, pi]
t = m.hypot(this._ssig1, this._csig1);
this._ssig1 /= t; this._csig1 /= t;
// norm(this._somg1, this._comg1); -- don't need to normalize!
this._k2 = m.sq(this._calp0) * geod._ep2;
eps = this._k2 / (2 * (1 + Math.sqrt(1 + this._k2)) + this._k2);
if (this.caps & g.CAP_C1) {
this._A1m1 = g.A1m1f(eps);
this._C1a = new Array(g.nC1_ + 1);
g.C1f(eps, this._C1a);
this._B11 = g.SinCosSeries(true, this._ssig1, this._csig1, this._C1a);
s = Math.sin(this._B11); c = Math.cos(this._B11);
// tau1 = sig1 + B11
this._stau1 = this._ssig1 * c + this._csig1 * s;
this._ctau1 = this._csig1 * c - this._ssig1 * s;
// Not necessary because C1pa reverts C1a
// _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa);
}
if (this.caps & g.CAP_C1p) {
this._C1pa = new Array(g.nC1p_ + 1);
g.C1pf(eps, this._C1pa);
}
if (this.caps & g.CAP_C2) {
this._A2m1 = g.A2m1f(eps);
this._C2a = new Array(g.nC2_ + 1);
g.C2f(eps, this._C2a);
this._B21 = g.SinCosSeries(true, this._ssig1, this._csig1, this._C2a);
}
if (this.caps & g.CAP_C3) {
this._C3a = new Array(g.nC3_);
geod.C3f(eps, this._C3a);
this._A3c = -this.f * this._salp0 * geod.A3f(eps);
this._B31 = g.SinCosSeries(true, this._ssig1, this._csig1, this._C3a);
}
if (this.caps & g.CAP_C4) {
this._C4a = new Array(g.nC4_); // all the elements of _C4a are used
geod.C4f(eps, this._C4a);
// Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
this._A4 = m.sq(this.a) * this._calp0 * this._salp0 * geod._e2;
this._B41 = g.SinCosSeries(false, this._ssig1, this._csig1, this._C4a);
}
this.a13 = this.s13 = NaN;
};
/**
* @summary Find the position on the line (general case).
* @param {bool} arcmode is the next parameter an arc length?
* @param {number} s12_a12 the (arcmode ? arc length : distance) from the
* first point to the second in (arcmode ? degrees : meters).
* @param {bitmask} [outmask = STANDARD] which results to include; this is
* subject to the capabilities of the object.
* @return {object} the requested results.
* @description The lat1, lon1, azi1, and a12 fields of the result are
* always set; s12 is included if arcmode is false. For details on the
* outmask parameter, see {@tutorial 2-interface}, "The outmask and caps
* parameters".
*/
l.GeodesicLine.prototype.GenPosition = function(arcmode, s12_a12,
outmask) {
var vals = {},
sig12, ssig12, csig12, B12, AB1, ssig2, csig2, tau12, s, c, serr,
omg12, lam12, lon12, E, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2,
B22, AB2, J12, t, B42, salp12, calp12;
if (!outmask) outmask = g.STANDARD;
else if (outmask === g.LONG_UNROLL) outmask |= g.STANDARD;
outmask &= this.caps & g.OUT_MASK;
vals.lat1 = this.lat1; vals.azi1 = this.azi1;
vals.lon1 = outmask & g.LONG_UNROLL ?
this.lon1 : m.AngNormalize(this.lon1);
if (arcmode)
vals.a12 = s12_a12;
else
vals.s12 = s12_a12;
if (!( arcmode || (this.caps & g.DISTANCE_IN & g.OUT_MASK) )) {
// Uninitialized or impossible distance calculation requested
vals.a12 = NaN;
return vals;
}
// Avoid warning about uninitialized B12.
B12 = 0; AB1 = 0;
if (arcmode) {
// Interpret s12_a12 as spherical arc length
sig12 = s12_a12 * m.degree;
t = m.sincosd(s12_a12); ssig12 = t.s; csig12 = t.c;
} else {
// Interpret s12_a12 as distance
tau12 = s12_a12 / (this._b * (1 + this._A1m1));
s = Math.sin(tau12);
c = Math.cos(tau12);
// tau2 = tau1 + tau12
B12 = -g.SinCosSeries(true,
this._stau1 * c + this._ctau1 * s,
this._ctau1 * c - this._stau1 * s,
this._C1pa);
sig12 = tau12 - (B12 - this._B11);
ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
if (Math.abs(this.f) > 0.01) {
// Reverted distance series is inaccurate for |f| > 1/100, so correct
// sig12 with 1 Newton iteration. The following table shows the
// approximate maximum error for a = WGS_a() and various f relative to
// GeodesicExact.
// erri = the error in the inverse solution (nm)
// errd = the error in the direct solution (series only) (nm)
// errda = the error in the direct solution
// (series + 1 Newton) (nm)
//
// f erri errd errda
// -1/5 12e6 1.2e9 69e6
// -1/10 123e3 12e6 765e3
// -1/20 1110 108e3 7155
// -1/50 18.63 200.9 27.12
// -1/100 18.63 23.78 23.37
// -1/150 18.63 21.05 20.26
// 1/150 22.35 24.73 25.83
// 1/100 22.35 25.03 25.31
// 1/50 29.80 231.9 30.44
// 1/20 5376 146e3 10e3
// 1/10 829e3 22e6 1.5e6
// 1/5 157e6 3.8e9 280e6
ssig2 = this._ssig1 * csig12 + this._csig1 * ssig12;
csig2 = this._csig1 * csig12 - this._ssig1 * ssig12;
B12 = g.SinCosSeries(true, ssig2, csig2, this._C1a);
serr = (1 + this._A1m1) * (sig12 + (B12 - this._B11)) -
s12_a12 / this._b;
sig12 = sig12 - serr / Math.sqrt(1 + this._k2 * m.sq(ssig2));
ssig12 = Math.sin(sig12); csig12 = Math.cos(sig12);
// Update B12 below
}
}
// sig2 = sig1 + sig12
ssig2 = this._ssig1 * csig12 + this._csig1 * ssig12;
csig2 = this._csig1 * csig12 - this._ssig1 * ssig12;
dn2 = Math.sqrt(1 + this._k2 * m.sq(ssig2));
if (outmask & (g.DISTANCE | g.REDUCEDLENGTH | g.GEODESICSCALE)) {
if (arcmode || Math.abs(this.f) > 0.01)
B12 = g.SinCosSeries(true, ssig2, csig2, this._C1a);
AB1 = (1 + this._A1m1) * (B12 - this._B11);
}
// sin(bet2) = cos(alp0) * sin(sig2)
sbet2 = this._calp0 * ssig2;
// Alt: cbet2 = hypot(csig2, salp0 * ssig2);
cbet2 = m.hypot(this._salp0, this._calp0 * csig2);
if (cbet2 === 0)
// I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
cbet2 = csig2 = g.tiny_;
// tan(alp0) = cos(sig2)*tan(alp2)
salp2 = this._salp0; calp2 = this._calp0 * csig2; // No need to normalize
if (arcmode && (outmask & g.DISTANCE))
vals.s12 = this._b * ((1 + this._A1m1) * sig12 + AB1);
if (outmask & g.LONGITUDE) {
// tan(omg2) = sin(alp0) * tan(sig2)
somg2 = this._salp0 * ssig2; comg2 = csig2; // No need to normalize
E = m.copysign(1, this._salp0);
// omg12 = omg2 - omg1
omg12 = outmask & g.LONG_UNROLL ?
E * (sig12 -
(Math.atan2(ssig2, csig2) -
Math.atan2(this._ssig1, this._csig1)) +
(Math.atan2(E * somg2, comg2) -
Math.atan2(E * this._somg1, this._comg1))) :
Math.atan2(somg2 * this._comg1 - comg2 * this._somg1,
comg2 * this._comg1 + somg2 * this._somg1);
lam12 = omg12 + this._A3c *
( sig12 + (g.SinCosSeries(true, ssig2, csig2, this._C3a) -
this._B31));
lon12 = lam12 / m.degree;
vals.lon2 = outmask & g.LONG_UNROLL ? this.lon1 + lon12 :
m.AngNormalize(m.AngNormalize(this.lon1) + m.AngNormalize(lon12));
}
if (outmask & g.LATITUDE)
vals.lat2 = m.atan2d(sbet2, this._f1 * cbet2);
if (outmask & g.AZIMUTH)
vals.azi2 = m.atan2d(salp2, calp2);
if (outmask & (g.REDUCEDLENGTH | g.GEODESICSCALE)) {
B22 = g.SinCosSeries(true, ssig2, csig2, this._C2a);
AB2 = (1 + this._A2m1) * (B22 - this._B21);
J12 = (this._A1m1 - this._A2m1) * sig12 + (AB1 - AB2);
if (outmask & g.REDUCEDLENGTH)
// Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
// accurate cancellation in the case of coincident points.
vals.m12 = this._b * (( dn2 * (this._csig1 * ssig2) -
this._dn1 * (this._ssig1 * csig2)) -
this._csig1 * csig2 * J12);
if (outmask & g.GEODESICSCALE) {
t = this._k2 * (ssig2 - this._ssig1) * (ssig2 + this._ssig1) /
(this._dn1 + dn2);
vals.M12 = csig12 +
(t * ssig2 - csig2 * J12) * this._ssig1 / this._dn1;
vals.M21 = csig12 -
(t * this._ssig1 - this._csig1 * J12) * ssig2 / dn2;
}
}
if (outmask & g.AREA) {
B42 = g.SinCosSeries(false, ssig2, csig2, this._C4a);
if (this._calp0 === 0 || this._salp0 === 0) {
// alp12 = alp2 - alp1, used in atan2 so no need to normalize
salp12 = salp2 * this.calp1 - calp2 * this.salp1;
calp12 = calp2 * this.calp1 + salp2 * this.salp1;
} else {
// tan(alp) = tan(alp0) * sec(sig)
// tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
// = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
// If csig12 > 0, write
// csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
// else
// csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
// No need to normalize
salp12 = this._calp0 * this._salp0 *
(csig12 <= 0 ? this._csig1 * (1 - csig12) + ssig12 * this._ssig1 :
ssig12 * (this._csig1 * ssig12 / (1 + csig12) + this._ssig1));
calp12 = m.sq(this._salp0) + m.sq(this._calp0) * this._csig1 * csig2;
}
vals.S12 = this._c2 * Math.atan2(salp12, calp12) +
this._A4 * (B42 - this._B41);
}
if (!arcmode)
vals.a12 = sig12 / m.degree;
return vals;
};
/**
* @summary Find the position on the line given s12.
* @param {number} s12 the distance from the first point to the second in
* meters.
* @param {bitmask} [outmask = STANDARD] which results to include; this is
* subject to the capabilities of the object.
* @return {object} the requested results.
* @description The lat1, lon1, azi1, s12, and a12 fields of the result are
* always set; s12 is included if arcmode is false. For details on the
* outmask parameter, see {@tutorial 2-interface}, "The outmask and caps
* parameters".
*/
l.GeodesicLine.prototype.Position = function(s12, outmask) {
return this.GenPosition(false, s12, outmask);
};
/**
* @summary Find the position on the line given a12.
* @param {number} a12 the arc length from the first point to the second in
* degrees.
* @param {bitmask} [outmask = STANDARD] which results to include; this is
* subject to the capabilities of the object.
* @return {object} the requested results.
* @description The lat1, lon1, azi1, and a12 fields of the result are
* always set. For details on the outmask parameter, see {@tutorial
* 2-interface}, "The outmask and caps parameters".
*/
l.GeodesicLine.prototype.ArcPosition = function(a12, outmask) {
return this.GenPosition(true, a12, outmask);
};
/**
* @summary Specify position of point 3 in terms of either distance or arc
* length.
* @param {bool} arcmode boolean flag determining the meaning of the second
* parameter; if arcmode is false, then the GeodesicLine object must have
* been constructed with caps |= DISTANCE_IN.
* @param {number} s13_a13 if arcmode is false, this is the distance from
* point 1 to point 3 (meters); otherwise it is the arc length from
* point 1 to point 3 (degrees); it can be negative.
*/
l.GeodesicLine.prototype.GenSetDistance = function(arcmode, s13_a13) {
if (arcmode)
this.SetArc(s13_a13);
else
this.SetDistance(s13_a13);
};
/**
* @summary Specify position of point 3 in terms distance.
* @param {number} s13 the distance from point 1 to point 3 (meters); it
* can be negative.
*/
l.GeodesicLine.prototype.SetDistance = function(s13) {
var r;
this.s13 = s13;
r = this.GenPosition(false, this.s13, g.ARC);
this.a13 = 0 + r.a12; // the 0+ converts undefined into NaN
};
/**
* @summary Specify position of point 3 in terms of arc length.
* @param {number} a13 the arc length from point 1 to point 3 (degrees);
* it can be negative.
*/
l.GeodesicLine.prototype.SetArc = function(a13) {
var r;
this.a13 = a13;
r = this.GenPosition(true, this.a13, g.DISTANCE);
this.s13 = 0 + r.s12; // the 0+ converts undefined into NaN
};
})(geodesic.Geodesic, geodesic.GeodesicLine, geodesic.Math);