جافا اسكريبت بالتفصيل
/*
2 * Geodesic routines from
GeographicLib translated to JavaScript.
See
3 *
https://geographiclib.sourceforge.io/JavaScript/doc
4 *
5 * The algorithms are derived
in
6 *
7 * Charles F. F. Karney,
8 * Algorithms for geodesics, J. Geodesy 87,
43-55 (2013),
9 * https://doi.org/10.1007/s00190-012-0578-z
10 * Addenda:
https://geographiclib.sourceforge.io/geod-addenda.html
11 *
12 * This file is the
concatenation and compression of the JavaScript files in
13 * src/geodesic in the source
tree for geographiclib-js.
14 *
15 * Copyright (c) Charles
Karney (2011-2022) <karney@alum.mit.edu> and licensed
16 * under the MIT/X11
License. For more information, see
17 *
https://geographiclib.sourceforge.io/
18 *
19 * Version: 2.1.1
20 * Date: 2024-08-18
21 */
22
23(function(cb) {
24/*
25 * Math.js
26 * Transcription of Math.hpp,
Constants.hpp, and Accumulator.hpp into
27 * JavaScript.
28 *
29 * Copyright (c) Charles
Karney (2011-2021) <karney@alum.mit.edu> and licensed
30 * under the MIT/X11
License. For more information, see
31 *
https://geographiclib.sourceforge.io/
32 */
33
34/**
35 * @namespace geodesic
36 * @description The parent
namespace for the following modules:
37 * - {@link
module:geodesic/Geodesic geodesic/Geodesic} The main
38 * engine for solving geodesic problems via the
39 * {@link module:geodesic/Geodesic.Geodesic
Geodesic} class.
40 * - {@link
module:geodesic/GeodesicLine geodesic/GeodesicLine}
41 * computes points along a single geodesic line
via the
42 * {@link module:geodesic/GeodesicLine.GeodesicLine
GeodesicLine}
43 * class.
44 * - {@link
module:geodesic/PolygonArea geodesic/PolygonArea}
45 * computes the area of a geodesic polygon via
the
46 * {@link
module:geodesic/PolygonArea.PolygonArea PolygonArea}
47 * class.
48 * - {@link
module:geodesic/Constants geodesic/Constants} defines
49 * constants specifying the version numbers and
the parameters for the WGS84
50 * ellipsoid.
51 *
52 * The following modules are
used internally by the package:
53 * - {@link
module:geodesic/Math geodesic/Math} defines various
54 * mathematical functions.
55 * - {@link
module:geodesic/Accumulator geodesic/Accumulator}
56 * interally used by
57 * {@link
module:geodesic/PolygonArea.PolygonArea PolygonArea} (via the
58 * {@link
module:geodesic/Accumulator.Accumulator Accumulator} class)
59 * for summing the contributions to the area of
a polygon.
60 */
61
62// To allow swap via [y, x] =
[x, y]
63/* jshint esversion: 6 */
64
65var geodesic = {};
66geodesic.Constants = {};
67geodesic.Math = {};
68geodesic.Accumulator = {};
69
70(function(
71 /**
72 * @exports geodesic/Constants
73 * @description Define constants defining the
version and WGS84 parameters.
74 */
75 c) {
76 "use strict";
77
78 /**
79 * @constant
80 * @summary WGS84 parameters.
81 * @property {number} a the equatorial radius
(meters).
82 * @property {number} f the flattening.
83 */
84 c.WGS84 = { a: 6378137, f: 1/298.257223563 };
85 /**
86 * @constant
87 * @summary an array of version numbers.
88 * @property {number} major the major version
number.
89 * @property {number} minor the minor version
number.
90 * @property {number} patch the patch number.
91 */
92 c.version = { major: 2,
93 minor: 1,
94 patch: 1 };
95 /**
96 * @constant
97 * @summary version string
98 */
99 c.version_string = "2.1.1";
100})(geodesic.Constants);
101
102(function(
103 /**
104 * @exports geodesic/Math
105 * @description Some useful mathematical
constants and functions (mainly for
106 *
internal use).
107 */
108 m) {
109 "use strict";
110
111 /**
112 * @summary The number of digits of precision
in floating-point numbers.
113 * @constant {number}
114 */
115 m.digits = 53;
116 /**
117 * @summary The machine epsilon.
118 * @constant {number}
119 */
120 m.epsilon = Math.pow(0.5, m.digits - 1);
121 /**
122 * @summary The factor to convert degrees to
radians.
123 * @constant {number}
124 */
125 m.degree = Math.PI/180;
126
127 /**
128 * @summary Square a number.
129 * @param {number} x the number.
130 * @return {number} the square.
131 */
132 m.sq = function(x) { return x * x; };
133
134 /**
135 * @summary The hypotenuse function.
136 * @param {number} x the first side.
137 * @param {number} y the second side.
138 * @return {number} the hypotenuse.
139 */
140 m.hypot = function(x, y) {
141 // Built in Math.hypot give incorrect
results from GeodSolve92.
142 return Math.sqrt(x*x + y*y);
143 };
144
145 /**
146 * @summary Cube root function.
147 * @param {number} x the argument.
148 * @return {number} the real cube root.
149 */
150 m.cbrt = Math.cbrt || function(x) {
151 var y = Math.pow(Math.abs(x), 1/3);
152 return x > 0 ? y : (x < 0 ? -y : x);
153 };
154
155 /**
156 * @summary The log1p function.
157 * @param {number} x the argument.
158 * @return {number} log(1 + x).
159 */
160 m.log1p = Math.log1p || function(x) {
161 var y = 1 + x,
162 z = y - 1;
163 // Here's the explanation for this magic: y
= 1 + z, exactly, and z
164 // approx x, thus log(y)/z (which is nearly
constant near z = 0) returns
165 // a good approximation to the true log(1 +
x)/x. The multiplication x *
166 // (log(y)/z) introduces little additional
error.
167 return z === 0 ? x : x * Math.log(y) / z;
168 };
169
170 /**
171 * @summary Inverse hyperbolic tangent.
172 * @param {number} x the argument.
173 * @return {number}
tanh<sup>−1</sup> x.
174 */
175 m.atanh = Math.atanh || function(x) {
176 var y = Math.abs(x); // Enforce odd parity
177 y = m.log1p(2 * y/(1 - y))/2;
178 return x > 0 ? y : (x < 0 ? -y : x);
179 };
180
181 /**
182 * @summary Copy the sign.
183 * @param {number} x gives the magitude of
the result.
184 * @param {number} y gives the sign of the
result.
185 * @return {number} value with the magnitude
of x and with the sign of y.
186 */
187 m.copysign = function(x, y) {
188 return Math.abs(x) * (y < 0 || (y === 0
&& 1/y < 0) ? -1 : 1);
189 };
190
191 /**
192 * @summary An error-free sum.
193 * @param {number} u
194 * @param {number} v
195 * @return {object} sum with sum.s = round(u
+ v) and sum.t is u + v −
196 * round(u
+ v)
197 */
198 m.sum = function(u, v) {
199 var s = u + v,
200 up = s - v,
201 vpp = s - up,
202 t;
203 up -= u;
204 vpp -= v;
205 // if s = 0, then t = 0 and give t the same
sign as s
206 t = s ? 0 - (up + vpp) : s;
207 // u + v = s
+ t
208 //
= round(u + v) + t
209 return {s: s, t: t};
210 };
211
212 /**
213 * @summary Evaluate a polynomial.
214 * @param {integer} N the order of the
polynomial.
215 * @param {array} p the coefficient array (of
size N + 1) (leading
216 *
order coefficient first)
217 * @param {number} x the variable.
218 * @return {number} the value of the
polynomial.
219 */
220 m.polyval = function(N, p, s, x) {
221 var y = N < 0 ? 0 : p[s++];
222 while (--N >= 0) y = y * x + p[s++];
223 return y;
224 };
225
226 /**
227 * @summary Coarsen a value close to zero.
228 * @param {number} x
229 * @return {number} the coarsened value.
230 */
231 m.AngRound = function(x) {
232 // The makes the smallest gap in x = 1/16 -
nextafter(1/16, 0) = 1/2^57 for
233 // reals = 0.7 pm on the earth if x is an
angle in degrees. (This is about
234 // 1000 times more resolution than we get
with angles around 90 degrees.)
235 // We use this to avoid having to deal with
near singular cases when x is
236 // non-zero but tiny (e.g., 1.0e-200). This converts -0 to +0; however
237 // tiny negative numbers get converted to
-0.
238 var z = 1/16,
239 y = Math.abs(x);
240 // The compiler mustn't
"simplify" z - (z - y) to y
241 y = y < z ? z - (z - y) : y;
242 return m.copysign(y, x);
243 };
244
245 /**
246 * @summary The remainder function.
247 * @param {number} x the numerator of the
division
248 * @param {number} y the denominator of the
division
249 * @return {number} the remainder in the
range [−y/2, y/2].
250 * <p>
251 * The range of x is unrestricted; y must be
positive.
252 */
253 m.remainder = function(x, y) {
254 x %= y;
255 return x < -y/2 ? x + y : (x < y/2 ?
x : x - y);
256 };
257
258 /**
259 * @summary Normalize an angle.
260 * @param {number} x the angle in degrees.
261 * @return {number} the angle reduced to the
range [−180°,
262 *
180°].
263 *
264 * The range of x is unrestricted. If the result is ±0° or
265 * ±180° then the sign is
the sign of \e x.
266 */
267 m.AngNormalize = function(x) {
268 // Place angle in [-180, 180].
269 var y = m.remainder(x, 360);
270 return Math.abs(y) === 180 ?
m.copysign(180, x) : y;
271 };
272
273 /**
274 * @summary Normalize a latitude.
275 * @param {number} x the angle in degrees.
276 * @return {number} x if it is in the range
[−90°, 90°],
277 *
otherwise return NaN.
278 */
279 m.LatFix = function(x) {
280 // Replace angle with NaN if outside [-90,
90].
281 return Math.abs(x) > 90 ? NaN : x;
282 };
283
284 /**
285 * @summary The exact difference of two
angles reduced to [−180°,
286 *
180°]
287 * @param {number} x the first angle in
degrees.
288 * @param {number} y the second angle in
degrees.
289 * @return {object} diff the exact
difference, diff.d + diff.e = y − x
290 *
mod 360°.
291 *
292 * This computes z = y − x exactly,
reduced to
293 * [−180°, 180°];
and then sets z = d + e where d
294 * is the nearest representable number to z
and e is the truncation
295 * error.
If z = ±0° or ±180°, then the sign
of
296 * d is given by the sign of y −
x. The maximum absolute
297 * value of e is
2<sup>−26</sup> (for doubles).
298 */
299 m.AngDiff = function(x, y) {
300 // Compute y - x and reduce to [-180,180]
accurately.
301 var r = m.sum(m.remainder(-x, 360),
m.remainder(y, 360)), d, e;
302 r = m.sum(m.remainder(r.s, 360), r.t);
303 d = r.s;
304 e = r.t;
305 // Fix the sign if d = -180, 0, 180.
306 if (d === 0 || Math.abs(d) === 180)
307 // If e == 0, take sign from y - x
308 // else (e != 0, implies d = +/-180), d
and e must have opposite signs
309 d = m.copysign(d, e === 0 ? y - x : -e);
310 return {d: d, e: e};
311 };
312
313 /**
314 * @summary Evaluate the sine and cosine
function with the argument in
315 *
degrees
316 * @param {number} x in degrees.
317 * @return {object} r with r.s = sin(x) and
r.c = cos(x).
318 */
319 m.sincosd = function(x) {
320 // In order to minimize round-off errors,
this function exactly reduces
321 // the argument to the range [-45, 45]
before converting it to radians.
322 var d, r, q, s, c, sinx, cosx;
323 d = x % 360;
324 q = Math.round(d / 90); //
If d is NaN this returns NaN
325 d -= 90 * q;
326 // now abs(d) <= 45
327 r = d * this.degree;
328 // Possibly could call the gnu extension
sincos
329 s = Math.sin(r); c = Math.cos(r);
330 if (Math.abs(d) === 45) {
331 c = Math.sqrt(0.5);
332 s = m.copysign(c, r);
333 } else if (Math.abs(d) === 30) {
334 c = Math.sqrt(0.75);
335 s = m.copysign(0.5, r);
336 }
337 switch (q & 3) {
338 case 0:
sinx = s; cosx = c; break;
339 case 1:
sinx = c; cosx = -s; break;
340 case 2:
sinx = -s; cosx = -c; break;
341 default: sinx = -c; cosx = s; break; // case 3
342 }
343 cosx += 0;
344 if (sinx === 0) sinx = m.copysign(sinx, x);
345 return {s: sinx, c: cosx};
346 };
347
348 /**
349 * @summary Evaluate the sine and cosine with
reduced argument plus correction
350 *
351 * @param {number} x reduced angle in
degrees.
352 * @param {number} t correction in degrees.
353 * @return {object} r with r.s = sin(x + t)
and r.c = cos(x + t).
354 */
355 m.sincosde = function(x, t) {
356 // In order to minimize round-off errors,
this function exactly reduces
357 // the argument to the range [-45, 45]
before converting it to radians.
358 var d, r, q, s, c, sinx, cosx;
359 d = x % 360;
360 q = Math.round(d / 90); // If d is NaN this returns NaN
361 d = m.AngRound((d - 90 * q) + t);
362 // now abs(d) <= 45
363 r = d * this.degree;
364 // Possibly could call the gnu extension
sincos
365 s = Math.sin(r); c = Math.cos(r);
366 if (Math.abs(d) === 45) {
367 c = Math.sqrt(0.5);
368 s = m.copysign(c, r);
369 } else if (Math.abs(d) === 30) {
370 c = Math.sqrt(0.75);
371 s = m.copysign(0.5, r);
372 }
373 switch (q & 3) {
374 case 0:
sinx = s; cosx = c; break;
375 case 1:
sinx = c; cosx = -s; break;
376 case 2:
sinx = -s; cosx = -c; break;
377 default: sinx = -c; cosx = s; break; // case 3
378 }
379 cosx += 0;
380 if (sinx === 0) sinx = m.copysign(sinx,
x+t);
381 return {s: sinx, c: cosx};
382 };
383
384 /**
385 * @summary Evaluate the atan2 function with
the result in degrees
386 * @param {number} y
387 * @param {number} x
388 * @return atan2(y, x) in degrees, in the
range [−180°
389 *
180°].
390 */
391 m.atan2d = function(y, x) {
392 // In order to minimize round-off errors,
this function rearranges the
393 // arguments so that result of atan2 is in
the range [-pi/4, pi/4] before
394 // converting it to degrees and mapping the
result to the correct
395 // quadrant.
396 var q = 0, ang;
397 if (Math.abs(y) > Math.abs(x)) { [y, x]
= [x, y]; q = 2; } // swap(x, y)
398 if (m.copysign(1, x) < 0) { x = -x; ++q;
}
399 // here x >= 0 and x >= abs(y), so
angle is in [-pi/4, pi/4]
400 ang = Math.atan2(y, x) / this.degree;
401 switch (q) {
402 // Note that atan2d(-0.0, 1.0) will
return -0. However, we expect that
403 // atan2d will not be called with y =
-0. If need be, include
404 //
405 //
case 0: ang = 0 + ang; break;
406 //
407 // and handle mpfr as in AngRound.
408 case 1: ang = m.copysign(180, y) - ang;
break;
409 case 2: ang = 90 - ang; break;
410 case 3: ang = -90 + ang; break;
411 default: break;
412 }
413 return ang;
414 };
415})(geodesic.Math);
416
417(function(
418 /**
419 * @exports geodesic/Accumulator
420 * @description Accurate summation via the
421 *
{@link module:geodesic/Accumulator.Accumulator Accumulator} class
422 *
(mainly for internal use).
423 */
424 a, m) {
425 "use strict";
426
427 /**
428 * @class
429 * @summary Accurate summation of many
numbers.
430 * @classdesc This allows many numbers to be
added together with twice the
431 *
normal precision. In the
documentation of the member functions, sum
432 *
stands for the value currently held in the accumulator.
433 * @param {number | Accumulator} [y = 0] set sum = y.
434 */
435 a.Accumulator = function(y) {
436 this.Set(y);
437 };
438
439 /**
440 * @summary Set the accumulator to a number.
441 * @param {number | Accumulator} [y = 0] set
sum = y.
442 */
443 a.Accumulator.prototype.Set = function(y) {
444 if (!y) y = 0;
445 if (y.constructor === a.Accumulator) {
446 this._s = y._s;
447 this._t = y._t;
448 } else {
449 this._s = y;
450 this._t = 0;
451 }
452 };
453
454 /**
455 * @summary Add a number to the accumulator.
456 * @param {number} [y = 0] set sum += y.
457 */
458 a.Accumulator.prototype.Add = function(y) {
459 // Here's Shewchuk's solution...
460 // Accumulate starting at least significant
end
461 var u = m.sum(y, this._t),
462 v = m.sum(u.s, this._s);
463 u = u.t;
464 this._s = v.s;
465 this._t = v.t;
466 // Start is _s, _t decreasing and
non-adjacent. Sum is now (s + t + u)
467 // exactly with s, t, u non-adjacent and in
decreasing order (except
468 // for possible zeros). The following code tries to normalize the
469 // result.
Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
470 // _s).
The follow does an approximate job (and maintains the
471 // decreasing non-adjacent property). Here are two "failures" using
472 // 3-bit floats:
473 //
474 // Case 1: _s is not equal to round(s+t+u)
-- off by 1 ulp
475 // [12, -1] - 8 -> [4, 0, -1] -> [4,
-1] = 3 should be [3, 0] = 3
476 //
477 // Case 2: _s+_t is not as close to s+t+u
as it shold be
478 // [64, 5] + 4 -> [64, 8, 1] ->
[64, 8] = 72 (off by 1)
479 // should be [80, -7] = 73
(exact)
480 //
481 // "Fixing" these problems is
probably not worth the expense. The
482 // representation inevitably leads to small
errors in the accumulated
483 // values.
The additional errors illustrated here amount to 1 ulp of
484 // the less significant word during each
addition to the Accumulator
485 // and an additional possible error of 1
ulp in the reported sum.
486 //
487 // Incidentally, the "ideal"
representation described above is not
488 // canonical, because _s = round(_s + _t)
may not be true. For
489 // example, with 3-bit floats:
490 //
491 // [128, 16] + 1 -> [160, -16] -- 160 =
round(145).
492 // But [160, 0] - 16 -> [128, 16] -- 128
= round(144).
493 //
494 if (this._s === 0) // This implies t == 0,
495 this._s = u; // so result is u
496 else
497 this._t += u; // otherwise just accumulate u to
t.
498 };
499
500 /**
501 * @summary Return the result of adding a
number to sum (but
502 *
don't change sum).
503 * @param {number} [y = 0] the number to be
added to the sum.
504 * @return sum + y.
505 */
506 a.Accumulator.prototype.Sum = function(y) {
507 var b;
508 if (!y)
509 return this._s;
510 else {
511 b = new a.Accumulator(this);
512 b.Add(y);
513 return b._s;
514 }
515 };
516
517 /**
518 * @summary Set sum = −sum.
519 */
520 a.Accumulator.prototype.Negate = function() {
521 this._s *= -1;
522 this._t *= -1;
523 };
524
525 /**
526 * @summary Take the remainder
527 * @param {number} y the divisor of the
remainder operation.
528 * @return sum in range [−y/2,
y/2].
529 */
530 a.Accumulator.prototype.Remainder =
function(y) {
531 this._s = m.remainder(this._s, y);
532 this.Add(0);
533 };
534})(geodesic.Accumulator,
geodesic.Math);
535/*
536 * Geodesic.js
537 * Transcription of
Geodesic.[ch]pp into JavaScript.
538 *
539 * See the documentation for
the C++ class. The conversion is a
literal
540 * conversion from C++.
541 *
542 * The algorithms are
derived in
543 *
544 * Charles F. F. Karney,
545 * Algorithms for geodesics, J. Geodesy 87,
43-55 (2013);
546 * https://doi.org/10.1007/s00190-012-0578-z
547 * Addenda:
https://geographiclib.sourceforge.io/geod-addenda.html
548 *
549 * Copyright (c) Charles
Karney (2011-2022) <karney@alum.mit.edu> and licensed
550 * under the MIT/X11
License. For more information, see
551 *
https://geographiclib.sourceforge.io/
552 */
553
554// Load AFTER Math.js
555
556// To allow swap via [y, x]
= [x, y]
557/* jshint esversion: 6 */
558
559geodesic.Geodesic = {};
560geodesic.GeodesicLine = {};
561geodesic.PolygonArea = {};
562
563(function(
564 /**
565 * @exports geodesic/Geodesic
566 * @description Solve geodesic problems via
the
567 *
{@link module:geodesic/Geodesic.Geodesic Geodesic} class.
568 */
569 g, l, p, m, c) {
570 "use strict";
571
572 var GEOGRAPHICLIB_GEODESIC_ORDER = 6,
573 nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER,
574 nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER,
575 nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER,
576 nA3x_ = nA3_,
577 nC3x_, nC4x_,
578 maxit1_ = 20,
579 maxit2_ = maxit1_ + m.digits + 10,
580 tol0_ = m.epsilon,
581 tol1_ = 200 * tol0_,
582 tol2_ = Math.sqrt(tol0_),
583 tolb_ = tol0_,
584 xthresh_ = 1000 * tol2_,
585 CAP_NONE = 0,
586 CAP_ALL
= 0x1F,
587 OUT_ALL
= 0x7F80,
588 astroid,
589 A1m1f_coeff, C1f_coeff, C1pf_coeff,
590 A2m1f_coeff, C2f_coeff,
591 A3_coeff, C3_coeff, C4_coeff;
592
593 // N.B. Number.MIN_VALUE is denormalized;
divide by Number.EPSILON to get min
594 // normalized positive number.
595 g.tiny_ =
Math.sqrt(Number.MIN_VALUE/Number.EPSILON);
596 g.nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
597 g.nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
598 g.nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
599 g.nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
600 g.nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
601 nC3x_ = (g.nC3_ * (g.nC3_ - 1)) / 2;
602 nC4x_ = (g.nC4_ * (g.nC4_ + 1)) / 2;
603 g.CAP_C1
= 1<<0;
604 g.CAP_C1p
= 1<<1;
605 g.CAP_C2
= 1<<2;
606 g.CAP_C3
= 1<<3;
607 g.CAP_C4
= 1<<4;
608
609 g.NONE = 0;
610 g.ARC = 1<<6;
611 g.LATITUDE = 1<<7 | CAP_NONE;
612 g.LONGITUDE = 1<<8 | g.CAP_C3;
613 g.AZIMUTH = 1<<9 | CAP_NONE;
614 g.DISTANCE = 1<<10 | g.CAP_C1;
615 g.STANDARD = g.LATITUDE | g.LONGITUDE | g.AZIMUTH |
g.DISTANCE;
616 g.DISTANCE_IN = 1<<11 | g.CAP_C1 | g.CAP_C1p;
617 g.REDUCEDLENGTH = 1<<12 | g.CAP_C1 |
g.CAP_C2;
618 g.GEODESICSCALE = 1<<13 | g.CAP_C1 |
g.CAP_C2;
619 g.AREA = 1<<14 | g.CAP_C4;
620 g.ALL = OUT_ALL| CAP_ALL;
621 g.LONG_UNROLL = 1<<15;
622 g.OUT_MASK = OUT_ALL| g.LONG_UNROLL;
623
624 g.SinCosSeries = function(sinp, sinx, cosx,
c) {
625 // Evaluate
626 // y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
627 // sum(c[i] * cos((2*i+1) * x), i, 0,
n-1)
628 // using Clenshaw summation. N.B. c[0] is unused for sin series
629 // Approx operation count = (n + 5) mult
and (2 * n + 2) add
630 var k = c.length, // Point to one beyond last element
631 n = k - (sinp ? 1 : 0),
632 ar = 2 * (cosx - sinx) * (cosx + sinx),
// 2 * cos(2 * x)
633 y0 = n & 1 ? c[--k] : 0, y1 =
0; // accumulators for sum
634 // Now n is even
635 n = Math.floor(n/2);
636 while (n--) {
637 // Unroll loop x 2, so accumulators
return to their original role
638 y1 = ar * y0 - y1 + c[--k];
639 y0 = ar * y1 - y0 + c[--k];
640 }
641 return (sinp ? 2 * sinx * cosx * y0 : //
sin(2 * x) * y0
642 cosx * (y0 - y1)); // cos(x) * (y0 - y1)
643 };
644
645 astroid = function(x, y) {
646 // Solve
k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive
647 // root k.
This solution is adapted from Geocentric::Reverse.
648 var k,
649 p = m.sq(x),
650 q = m.sq(y),
651 r = (p + q - 1) / 6,
652 S, r2, r3, disc, u, T3, T, ang, v, uv,
w;
653 if ( !(q === 0 && r <= 0) ) {
654 // Avoid possible division by zero when r
= 0 by multiplying
655 // equations for s and t by r^3 and r,
resp.
656 S = p * q / 4; // S = r^3 * s
657 r2 = m.sq(r);
658 r3 = r * r2;
659 // The discriminant of the quadratic
equation for T3. This is
660 // zero on the evolute curve
p^(1/3)+q^(1/3) = 1
661 disc = S * (S + 2 * r3);
662 u = r;
663 if (disc >= 0) {
664 T3 = S + r3;
665 // Pick the sign on the sqrt to
maximize abs(T3). This
666 // minimizes loss of precision due to
cancellation. The
667 // result is unchanged because of the
way the T is used
668 // in definition of u.
669 T3 += T3 < 0 ? -Math.sqrt(disc) :
Math.sqrt(disc); // T3 = (r * t)^3
670 // N.B. cbrt always returns the real
root. cbrt(-8) = -2.
671 T = m.cbrt(T3); // T = r * t
672 // T can be zero; but then r2 / T ->
0.
673 u += T + (T !== 0 ? r2 / T : 0);
674 } else {
675 // T is complex, but the way u is
defined the result is real.
676 ang = Math.atan2(Math.sqrt(-disc), -(S
+ r3));
677 // There are three possible cube
roots. We choose the
678 // root which avoids cancellation. Note that disc < 0
679 // implies that r < 0.
680 u += 2 * r * Math.cos(ang / 3);
681 }
682 v = Math.sqrt(m.sq(u) + q); // guaranteed positive
683 // Avoid loss of accuracy when u < 0.
684 uv = u < 0 ? q / (v - u) : u + v; //
u+v, guaranteed positive
685 w = (uv - q) / (2 * v); // positive?
686 // Rearrange expression for k to avoid
loss of accuracy due to
687 // subtraction. Division by 0 not possible because uv > 0,
w >= 0.
688 k = uv / (Math.sqrt(uv + m.sq(w)) + w);
// guaranteed positive
689 } else { // q == 0
&& r <= 0
690 // y = 0 with |x| <= 1. Handle this case directly.
691 // for y small, positive root is k =
abs(y)/sqrt(1-x^2)
692 k = 0;
693 }
694 return k;
695 };
696
697 A1m1f_coeff = [
698 // (1-eps)*A1-1, polynomial in eps2 of
order 3
699 +1, 4, 64, 0, 256
700 ];
701
702 // The scale factor A1-1 = mean value of
(d/dsigma)I1 - 1
703 g.A1m1f = function(eps) {
704 var p = Math.floor(nA1_/2),
705 t = m.polyval(p, A1m1f_coeff, 0,
m.sq(eps)) / A1m1f_coeff[p + 1];
706 return (t + eps) / (1 - eps);
707 };
708
709 C1f_coeff = [
710 // C1[1]/eps^1, polynomial in eps2 of order
2
711 -1, 6, -16, 32,
712 // C1[2]/eps^2, polynomial in eps2 of order
2
713 -9, 64, -128, 2048,
714 // C1[3]/eps^3, polynomial in eps2 of order
1
715 +9, -16, 768,
716 // C1[4]/eps^4, polynomial in eps2 of order
1
717 +3, -5, 512,
718 // C1[5]/eps^5, polynomial in eps2 of order
0
719 -7, 1280,
720 // C1[6]/eps^6, polynomial in eps2 of order
0
721 -7, 2048
722 ];
723
724 // The coefficients C1[l] in the Fourier
expansion of B1
725 g.C1f = function(eps, c) {
726 var eps2 = m.sq(eps),
727 d = eps,
728 o = 0,
729 l, p;
730 for (l = 1; l <= g.nC1_; ++l) { // l is index of C1p[l]
731 p = Math.floor((g.nC1_ - l) / 2); //
order of polynomial in eps^2
732 c[l] = d * m.polyval(p, C1f_coeff, o,
eps2) / C1f_coeff[o + p + 1];
733 o += p + 2;
734 d *= eps;
735 }
736 };
737
738 C1pf_coeff = [
739 // C1p[1]/eps^1, polynomial in eps2 of
order 2
740 +205, -432, 768, 1536,
741 // C1p[2]/eps^2, polynomial in eps2 of
order 2
742 +4005, -4736, 3840, 12288,
743 // C1p[3]/eps^3, polynomial in eps2 of
order 1
744 -225, 116, 384,
745 // C1p[4]/eps^4, polynomial in eps2 of
order 1
746 -7173, 2695, 7680,
747 // C1p[5]/eps^5, polynomial in eps2 of
order 0
748 +3467, 7680,
749 // C1p[6]/eps^6, polynomial in eps2 of
order 0
750 +38081, 61440
751 ];
752
753 // The coefficients C1p[l] in the Fourier
expansion of B1p
754 g.C1pf = function(eps, c) {
755 var eps2 = m.sq(eps),
756 d = eps,
757 o = 0,
758 l, p;
759 for (l = 1; l <= g.nC1p_; ++l) { // l is index of C1p[l]
760 p = Math.floor((g.nC1p_ - l) / 2); //
order of polynomial in eps^2
761 c[l] = d * m.polyval(p, C1pf_coeff, o,
eps2) / C1pf_coeff[o + p + 1];
762 o += p + 2;
763 d *= eps;
764 }
765 };
766
767 A2m1f_coeff = [
768 // (eps+1)*A2-1, polynomial in eps2 of
order 3
769 -11, -28, -192, 0, 256
770 ];
771
772 // The scale factor A2-1 = mean value of
(d/dsigma)I2 - 1
773 g.A2m1f = function(eps) {
774 var p = Math.floor(nA2_/2),
775 t = m.polyval(p, A2m1f_coeff, 0,
m.sq(eps)) / A2m1f_coeff[p + 1];
776 return (t - eps) / (1 + eps);
777 };
778
779 C2f_coeff = [
780 // C2[1]/eps^1, polynomial in eps2 of order
2
781 +1, 2, 16, 32,
782 // C2[2]/eps^2, polynomial in eps2 of order
2
783 +35, 64, 384, 2048,
784 // C2[3]/eps^3, polynomial in eps2 of order
1
785 +15, 80, 768,
786 // C2[4]/eps^4, polynomial in eps2 of order
1
787 +7, 35, 512,
788 // C2[5]/eps^5, polynomial in eps2 of order
0
789 +63, 1280,
790 // C2[6]/eps^6, polynomial in eps2 of order
0
791 +77, 2048
792 ];
793
794 // The coefficients C2[l] in the Fourier
expansion of B2
795 g.C2f = function(eps, c) {
796 var eps2 = m.sq(eps),
797 d = eps,
798 o = 0,
799 l, p;
800 for (l = 1; l <= g.nC2_; ++l) { // l is index of C2[l]
801 p = Math.floor((g.nC2_ - l) / 2); //
order of polynomial in eps^2
802 c[l] = d * m.polyval(p, C2f_coeff, o,
eps2) / C2f_coeff[o + p + 1];
803 o += p + 2;
804 d *= eps;
805 }
806 };
807
808 /**
809 * @class
810 * @property {number} a the equatorial radius
(meters).
811 * @property {number} f the flattening.
812 * @summary Initialize a Geodesic object for
a specific ellipsoid.
813 * @classdesc Performs geodesic calculations
on an ellipsoid of revolution.
814 *
The routines for solving the direct and inverse problems return an
815 *
object with some of the following fields set: lat1, lon1, azi1, lat2,
816 *
lon2, azi2, s12, a12, m12, M12, M21, S12. See {@tutorial 2-interface},
817 *
section "The results".
818 * @example
819 * var geodesic =
require("geographiclib-geodesic"),
820 *
geod = geodesic.Geodesic.WGS84;
821 * var inv = geod.Inverse(1,2,3,4);
822 * console.log("lat1 = " + inv.lat1
+ ", lon1 = " + inv.lon1 +
823 *
", lat2 = " + inv.lat2 + ", lon2 = " + inv.lon2 +
824 *
",\nazi1 = " + inv.azi1 + ", azi2 = " + inv.azi2 +
825 *
", s12 = " + inv.s12);
826 * @param {number} a the equatorial radius of
the ellipsoid (meters).
827 * @param {number} f the flattening of the
ellipsoid. Setting f = 0 gives
828 * a
sphere (on which geodesics are great circles).
Negative f gives a
829 *
prolate ellipsoid.
830 * @throws an error if the parameters are
illegal.
831 */
832 g.Geodesic = function(a, f) {
833 this.a = a;
834 this.f = f;
835 this._f1 = 1 - this.f;
836 this._e2 = this.f * (2 - this.f);
837 this._ep2 = this._e2 / m.sq(this._f1); //
e2 / (1 - e2)
838 this._n = this.f / ( 2 - this.f);
839 this._b = this.a * this._f1;
840 // authalic radius squared
841 this._c2 = (m.sq(this.a) + m.sq(this._b) *
842 (this._e2 === 0 ? 1 :
843 (this._e2 > 0 ?
m.atanh(Math.sqrt(this._e2)) :
844
Math.atan(Math.sqrt(-this._e2))) /
845
Math.sqrt(Math.abs(this._e2))))/2;
846 // The sig12 threshold for "really
short". Using the auxiliary sphere
847 // solution with dnm computed at (bet1 +
bet2) / 2, the relative error in
848 // the azimuth consistency check is sig12^2
* abs(f) * min(1, 1-f/2) / 2.
849 // (Error measured for 1/100 < b/a <
100 and abs(f) >= 1/1000. For a given
850 // f and sig12, the max error occurs for
lines near the pole. If the old
851 // rule for computing dnm = (dn1 + dn2)/2
is used, then the error increases
852 // by a factor of 2.) Setting this equal to epsilon gives sig12 =
etol2.
853 // Here 0.1 is a safety factor (error
decreased by 100) and max(0.001,
854 // abs(f)) stops etol2 getting too large in
the nearly spherical case.
855 this._etol2 = 0.1 * tol2_ /
856 Math.sqrt( Math.max(0.001,
Math.abs(this.f)) *
857 Math.min(1, 1 - this.f/2) / 2
);
858 if (!(isFinite(this.a) && this.a
> 0))
859 throw new Error("Equatorial radius
is not positive");
860 if (!(isFinite(this._b) && this._b
> 0))
861 throw new Error("Polar semi-axis is
not positive");
862 this._A3x = new Array(nA3x_);
863 this._C3x = new Array(nC3x_);
864 this._C4x = new Array(nC4x_);
865 this.A3coeff();
866 this.C3coeff();
867 this.C4coeff();
868 };
869
870 A3_coeff = [
871 // A3, coeff of eps^5, polynomial in n of
order 0
872 -3, 128,
873 // A3, coeff of eps^4, polynomial in n of
order 1
874 -2, -3, 64,
875 // A3, coeff of eps^3, polynomial in n of
order 2
876 -1, -3, -1, 16,
877 // A3, coeff of eps^2, polynomial in n of
order 2
878 +3, -1, -2, 8,
879 // A3, coeff of eps^1, polynomial in n of
order 1
880 +1, -1, 2,
881 // A3, coeff of eps^0, polynomial in n of
order 0
882 +1, 1
883 ];
884
885 // The scale factor A3 = mean value of
(d/dsigma)I3
886 g.Geodesic.prototype.A3coeff = function() {
887 var o = 0, k = 0,
888 j, p;
889 for (j = nA3_ - 1; j >= 0; --j) { //
coeff of eps^j
890 p = Math.min(nA3_ - j - 1, j); // order of polynomial in n
891 this._A3x[k++] = m.polyval(p, A3_coeff,
o, this._n) /
892 A3_coeff[o + p + 1];
893 o += p + 2;
894 }
895 };
896
897 C3_coeff = [
898 // C3[1], coeff of eps^5, polynomial in n
of order 0
899 +3, 128,
900 // C3[1], coeff of eps^4, polynomial in n
of order 1
901 +2, 5, 128,
902 // C3[1], coeff of eps^3, polynomial in n
of order 2
903 -1, 3, 3, 64,
904 // C3[1], coeff of eps^2, polynomial in n
of order 2
905 -1, 0, 1, 8,
906 // C3[1], coeff of eps^1, polynomial in n
of order 1
907 -1, 1, 4,
908 // C3[2], coeff of eps^5, polynomial in n
of order 0
909 +5, 256,
910 // C3[2], coeff of eps^4, polynomial in n
of order 1
911 +1, 3, 128,
912 // C3[2], coeff of eps^3, polynomial in n
of order 2
913 -3, -2, 3, 64,
914 // C3[2], coeff of eps^2, polynomial in n
of order 2
915 +1, -3, 2, 32,
916 // C3[3], coeff of eps^5, polynomial in n
of order 0
917 +7, 512,
918 // C3[3], coeff of eps^4, polynomial in n
of order 1
919 -10, 9, 384,
920 // C3[3], coeff of eps^3, polynomial in n
of order 2
921 +5, -9, 5, 192,
922 // C3[4], coeff of eps^5, polynomial in n
of order 0
923 +7, 512,
924 // C3[4], coeff of eps^4, polynomial in n
of order 1
925 -14, 7, 512,
926 // C3[5], coeff of eps^5, polynomial in n
of order 0
927 +21, 2560
928 ];
929
930 // The coefficients C3[l] in the Fourier
expansion of B3
931 g.Geodesic.prototype.C3coeff = function() {
932 var o = 0, k = 0,
933 l, j, p;
934 for (l = 1; l < g.nC3_; ++l) { // l is index of C3[l]
935 for (j = g.nC3_ - 1; j >= l; --j) { //
coeff of eps^j
936 p = Math.min(g.nC3_ - j - 1, j); // order of polynomial in n
937 this._C3x[k++] = m.polyval(p, C3_coeff,
o, this._n) /
938 C3_coeff[o + p + 1];
939 o += p + 2;
940 }
941 }
942 };
943
944 C4_coeff = [
945 // C4[0], coeff of eps^5, polynomial in n
of order 0
946 +97, 15015,
947 // C4[0], coeff of eps^4, polynomial in n
of order 1
948 +1088, 156, 45045,
949 // C4[0], coeff of eps^3, polynomial in n
of order 2
950 -224, -4784, 1573, 45045,
951 // C4[0], coeff of eps^2, polynomial in n
of order 3
952 -10656, 14144, -4576, -858, 45045,
953 // C4[0], coeff of eps^1, polynomial in n
of order 4
954 +64, 624, -4576, 6864, -3003, 15015,
955 // C4[0], coeff of eps^0, polynomial in n
of order 5
956 +100, 208, 572, 3432, -12012, 30030,
45045,
957 // C4[1], coeff of eps^5, polynomial in n
of order 0
958 +1, 9009,
959 // C4[1], coeff of eps^4, polynomial in n
of order 1
960 -2944, 468, 135135,
961 // C4[1], coeff of eps^3, polynomial in n
of order 2
962 +5792, 1040, -1287, 135135,
963 // C4[1], coeff of eps^2, polynomial in n
of order 3
964 +5952, -11648, 9152, -2574, 135135,
965 // C4[1], coeff of eps^1, polynomial in n
of order 4
966 -64, -624, 4576, -6864, 3003, 135135,
967 // C4[2], coeff of eps^5, polynomial in n
of order 0
968 +8, 10725,
969 // C4[2], coeff of eps^4, polynomial in n
of order 1
970 +1856, -936, 225225,
971 // C4[2], coeff of eps^3, polynomial in n
of order 2
972 -8448, 4992, -1144, 225225,
973 // C4[2], coeff of eps^2, polynomial in n
of order 3
974 -1440, 4160, -4576, 1716, 225225,
975 // C4[3], coeff of eps^5, polynomial in n
of order 0
976 -136, 63063,
977 // C4[3], coeff of eps^4, polynomial in n
of order 1
978 +1024, -208, 105105,
979 // C4[3], coeff of eps^3, polynomial in n
of order 2
980 +3584, -3328, 1144, 315315,
981 // C4[4], coeff of eps^5, polynomial in n
of order 0
982 -128, 135135,
983 // C4[4], coeff of eps^4, polynomial in n
of order 1
984 -2560, 832, 405405,
985 // C4[5], coeff of eps^5, polynomial in n
of order 0
986 +128, 99099
987 ];
988
989 g.Geodesic.prototype.C4coeff = function() {
990 var o = 0, k = 0,
991 l, j, p;
992 for (l = 0; l < g.nC4_; ++l) { // l is index of C4[l]
993 for (j = g.nC4_ - 1; j >= l; --j) { //
coeff of eps^j
994 p = g.nC4_ - j - 1; // order of polynomial in n
995 this._C4x[k++] = m.polyval(p, C4_coeff,
o, this._n) /
996 C4_coeff[o + p + 1];
997 o += p + 2;
998 }
999 }
1000 };
1001
1002 g.Geodesic.prototype.A3f = function(eps) {
1003 // Evaluate A3
1004 return m.polyval(nA3x_ - 1, this._A3x, 0,
eps);
1005 };
1006
1007 g.Geodesic.prototype.C3f = function(eps, c) {
1008 // Evaluate C3 coeffs
1009 // Elements c[1] thru c[nC3_ - 1] are set
1010 var mult = 1,
1011 o = 0,
1012 l, p;
1013 for (l = 1; l < g.nC3_; ++l) { // l is
index of C3[l]
1014 p = g.nC3_ - l - 1; // order of polynomial in eps
1015 mult *= eps;
1016 c[l] = mult * m.polyval(p, this._C3x, o,
eps);
1017 o += p + 1;
1018 }
1019 };
1020
1021 g.Geodesic.prototype.C4f = function(eps, c) {
1022 // Evaluate C4 coeffs
1023 // Elements c[0] thru c[g.nC4_ - 1] are set
1024 var mult = 1,
1025 o = 0,
1026 l, p;
1027 for (l = 0; l < g.nC4_; ++l) { // l is
index of C4[l]
1028 p = g.nC4_ - l - 1; // order of polynomial in eps
1029 c[l] = mult * m.polyval(p, this._C4x, o,
eps);
1030 o += p + 1;
1031 mult *= eps;
1032 }
1033 };
1034
1035 // return s12b, m12b, m0, M12, M21
1036 g.Geodesic.prototype.Lengths = function(eps,
sig12,
1037
ssig1, csig1, dn1, ssig2, csig2, dn2,
1038
cbet1, cbet2, outmask,
1039 C1a,
C2a) {
1040 // Return m12b = (reduced length)/_b; also
calculate s12b =
1041 // distance/_b, and m0 = coefficient of
secular term in
1042 // expression for reduced length.
1043 outmask &= g.OUT_MASK;
1044 var vals = {},
1045 m0x = 0, J12 = 0, A1 = 0, A2 = 0,
1046 B1, B2, l, csig12, t;
1047 if (outmask & (g.DISTANCE |
g.REDUCEDLENGTH | g.GEODESICSCALE)) {
1048 A1 = g.A1m1f(eps);
1049 g.C1f(eps, C1a);
1050 if (outmask & (g.REDUCEDLENGTH |
g.GEODESICSCALE)) {
1051 A2 = g.A2m1f(eps);
1052 g.C2f(eps, C2a);
1053 m0x = A1 - A2;
1054 A2 = 1 + A2;
1055 }
1056 A1 = 1 + A1;
1057 }
1058 if (outmask & g.DISTANCE) {
1059 B1 = g.SinCosSeries(true, ssig2, csig2,
C1a) -
1060 g.SinCosSeries(true, ssig1, csig1,
C1a);
1061 // Missing a factor of _b
1062 vals.s12b = A1 * (sig12 + B1);
1063 if (outmask & (g.REDUCEDLENGTH |
g.GEODESICSCALE)) {
1064 B2 = g.SinCosSeries(true, ssig2, csig2,
C2a) -
1065 g.SinCosSeries(true, ssig1, csig1,
C2a);
1066 J12 = m0x * sig12 + (A1 * B1 - A2 *
B2);
1067 }
1068 } else if (outmask & (g.REDUCEDLENGTH |
g.GEODESICSCALE)) {
1069 // Assume here that nC1_ >= nC2_
1070 for (l = 1; l <= g.nC2_; ++l)
1071 C2a[l] = A1 * C1a[l] - A2 * C2a[l];
1072 J12 = m0x * sig12 + (g.SinCosSeries(true,
ssig2, csig2, C2a) -
1073 g.SinCosSeries(true,
ssig1, csig1, C2a));
1074 }
1075 if (outmask & g.REDUCEDLENGTH) {
1076 vals.m0 = m0x;
1077 // Missing a factor of _b.
1078 // Add parens around (csig1 * ssig2) and
(ssig1 * csig2) to ensure
1079 // accurate cancellation in the case of
coincident points.
1080 vals.m12b = dn2 * (csig1 * ssig2) - dn1 *
(ssig1 * csig2) -
1081 csig1 * csig2 * J12;
1082 }
1083 if (outmask & g.GEODESICSCALE) {
1084 csig12 = csig1 * csig2 + ssig1 * ssig2;
1085 t = this._ep2 * (cbet1 - cbet2) * (cbet1
+ cbet2) / (dn1 + dn2);
1086 vals.M12 = csig12 + (t * ssig2 - csig2 *
J12) * ssig1 / dn1;
1087 vals.M21 = csig12 - (t * ssig1 - csig1 *
J12) * ssig2 / dn2;
1088 }
1089 return vals;
1090 };
1091
1092 // return sig12, salp1, calp1, salp2, calp2,
dnm
1093 g.Geodesic.prototype.InverseStart =
function(sbet1, cbet1, dn1,
1094
sbet2, cbet2, dn2,
1095
lam12, slam12, clam12,
1096
C1a, C2a) {
1097 // Return a starting point for Newton's
method in salp1 and calp1
1098 // (function value is -1). If Newton's method doesn't need to be
1099 // used, return also salp2 and calp2 and
function value is sig12.
1100 // salp2, calp2 only updated if return val
>= 0.
1101 var vals = {},
1102 // bet12 = bet2 - bet1 in [0, pi);
bet12a = bet2 + bet1 in (-pi, 0]
1103 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
1104 cbet12 = cbet2 * cbet1 + sbet2 * sbet1,
1105 sbet12a, shortline, omg12, sbetm2,
somg12, comg12, t, ssig12, csig12,
1106 x, y, lamscale, betscale, k2, eps,
cbet12a, bet12a, m12b, m0, nvals,
1107 k, omg12a, lam12x;
1108 vals.sig12 = -1; // Return value
1109 // Volatile declaration needed to fix
inverse cases
1110 // 88.202499451857 0 -88.202499451857
179.981022032992859592
1111 // 89.262080389218 0 -89.262080389218
179.992207982775375662
1112 // 89.333123580033 0 -89.333123580032997687
179.99295812360148422
1113 // which otherwise fail with g++ 4.4.4 x86
-O3
1114 sbet12a = sbet2 * cbet1;
1115 sbet12a += cbet2 * sbet1;
1116
1117 shortline = cbet12 >= 0 &&
sbet12 < 0.5 && cbet2 * lam12 < 0.5;
1118 if (shortline) {
1119 sbetm2 = m.sq(sbet1 + sbet2);
1120 // sin((bet1+bet2)/2)^2
1121 // =
(sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
1122 sbetm2 /= sbetm2 + m.sq(cbet1 + cbet2);
1123 vals.dnm = Math.sqrt(1 + this._ep2 *
sbetm2);
1124 omg12 = lam12 / (this._f1 * vals.dnm);
1125 somg12 = Math.sin(omg12); comg12 =
Math.cos(omg12);
1126 } else {
1127 somg12 = slam12; comg12 = clam12;
1128 }
1129
1130 vals.salp1 = cbet2 * somg12;
1131 vals.calp1 = comg12 >= 0 ?
1132 sbet12 + cbet2 * sbet1 * m.sq(somg12) /
(1 + comg12) :
1133 sbet12a - cbet2 * sbet1 * m.sq(somg12) /
(1 - comg12);
1134
1135 ssig12 = m.hypot(vals.salp1, vals.calp1);
1136 csig12 = sbet1 * sbet2 + cbet1 * cbet2 *
comg12;
1137 if (shortline && ssig12 <
this._etol2) {
1138 // really short lines
1139 vals.salp2 = cbet1 * somg12;
1140 vals.calp2 = sbet12 - cbet1 * sbet2 *
1141 (comg12 >= 0 ? m.sq(somg12) / (1 +
comg12) : 1 - comg12);
1142 // norm(vals.salp2, vals.calp2);
1143 t = m.hypot(vals.salp2, vals.calp2);
vals.salp2 /= t; vals.calp2 /= t;
1144 // Set return value
1145 vals.sig12 = Math.atan2(ssig12, csig12);
1146 } else if (Math.abs(this._n) > 0.1 || //
Skip astroid calc if too eccentric
1147 csig12 >= 0 ||
1148 ssig12 >= 6 *
Math.abs(this._n) * Math.PI * m.sq(cbet1)) {
1149 // Nothing to do, zeroth order spherical
approximation is OK
1150 } else {
1151 // Scale lam12 and bet2 to x, y
coordinate system where antipodal
1152 // point is at origin and singular point
is at y = 0, x = -1.
1153 lam12x = Math.atan2(-slam12, -clam12); //
lam12 - pi
1154 if (this.f >= 0) { // In fact f == 0 does not get here
1155 // x = dlong, y = dlat
1156 k2 = m.sq(sbet1) * this._ep2;
1157 eps = k2 / (2 * (1 + Math.sqrt(1 + k2))
+ k2);
1158 lamscale = this.f * cbet1 *
this.A3f(eps) * Math.PI;
1159 betscale = lamscale * cbet1;
1160
1161 x = lam12x / lamscale;
1162 y = sbet12a / betscale;
1163 } else { // f < 0
1164 // x = dlat, y = dlong
1165 cbet12a = cbet2 * cbet1 - sbet2 *
sbet1;
1166 bet12a = Math.atan2(sbet12a, cbet12a);
1167 // In the case of lon12 = 180, this
repeats a calculation made
1168 // in Inverse.
1169 nvals = this.Lengths(this._n, Math.PI +
bet12a,
1170 sbet1, -cbet1,
dn1, sbet2, cbet2, dn2,
1171 cbet1, cbet2,
g.REDUCEDLENGTH, C1a, C2a);
1172 m12b = nvals.m12b; m0 = nvals.m0;
1173 x = -1 + m12b / (cbet1 * cbet2 * m0 *
Math.PI);
1174 betscale = x < -0.01 ? sbet12a / x :
1175 -this.f * m.sq(cbet1) * Math.PI;
1176 lamscale = betscale / cbet1;
1177 y = lam12 / lamscale;
1178 }
1179
1180 if (y > -tol1_ && x > -1 -
xthresh_) {
1181 // strip near cut
1182 if (this.f >= 0) {
1183 vals.salp1 = Math.min(1, -x);
1184 vals.calp1 = -Math.sqrt(1 -
m.sq(vals.salp1));
1185 } else {
1186 vals.calp1 = Math.max(x > -tol1_ ?
0 : -1, x);
1187 vals.salp1 = Math.sqrt(1 -
m.sq(vals.calp1));
1188 }
1189 } else {
1190 // Estimate alp1, by solving the
astroid problem.
1191 //
1192 // Could estimate alpha1 = theta +
pi/2, directly, i.e.,
1193 //
calp1 = y/k; salp1 = -x/(1+k);
for f >= 0
1194 //
calp1 = x/(1+k); salp1 = -y/k;
for f < 0 (need to check)
1195 //
1196 // However, it's better to estimate
omg12 from astroid and use
1197 // spherical formula to compute
alp1. This reduces the mean number of
1198 // Newton iterations for astroid cases
from 2.24 (min 0, max 6) to 2.12
1199 // (min 0 max 5). The changes in the number of iterations are
as
1200 // follows:
1201 //
1202 // change percent
1203 //
1 5
1204 //
0 78
1205 //
-1 16
1206 //
-2 0.6
1207 //
-3 0.04
1208 //
-4 0.002
1209 //
1210 // The histogram of iterations is (m =
number of iterations estimating
1211 // alp1 directly, n = number of
iterations estimating via omg12, total
1212 // number of trials = 148605):
1213 //
1214 //
iter m n
1215 //
0 148 186
1216 //
1 13046 13845
1217 //
2 93315 102225
1218 //
3 36189 32341
1219 //
4 5396 7
1220 //
5 455 1
1221 //
6 56 0
1222 //
1223 // Because omg12 is near pi, estimate
work with omg12a = pi - omg12
1224 k = astroid(x, y);
1225 omg12a = lamscale * ( this.f >= 0 ?
-x * k/(1 + k) : -y * (1 + k)/k );
1226 somg12 = Math.sin(omg12a); comg12 =
-Math.cos(omg12a);
1227 // Update spherical estimate of alp1
using omg12 instead of
1228 // lam12
1229 vals.salp1 = cbet2 * somg12;
1230 vals.calp1 = sbet12a -
1231 cbet2 * sbet1 * m.sq(somg12) / (1 -
comg12);
1232 }
1233 }
1234 // Sanity check on starting guess. Backwards check allows NaN through.
1235 // jshint -W018
1236 if (!(vals.salp1 <= 0)) {
1237 // norm(vals.salp1, vals.calp1);
1238 t = m.hypot(vals.salp1, vals.calp1);
vals.salp1 /= t; vals.calp1 /= t;
1239 } else {
1240 vals.salp1 = 1; vals.calp1 = 0;
1241 }
1242 return vals;
1243 };
1244
1245 // return lam12, salp2, calp2, sig12, ssig1,
csig1, ssig2, csig2, eps,
1246 // domg12, dlam12,
1247 g.Geodesic.prototype.Lambda12 =
function(sbet1, cbet1, dn1,
1248
sbet2, cbet2, dn2,
1249
salp1, calp1, slam120, clam120,
1250
diffp, C1a, C2a, C3a) {
1251 var vals = {},
1252 t, salp0, calp0,
1253 somg1, comg1, somg2, comg2, somg12,
comg12, B312, eta, k2, nvals;
1254 if (sbet1 === 0 && calp1 === 0)
1255 // Break degeneracy of equatorial
line. This case has already been
1256 // handled.
1257 calp1 = -g.tiny_;
1258
1259 // sin(alp1) * cos(bet1) = sin(alp0)
1260 salp0 = salp1 * cbet1;
1261 calp0 = m.hypot(calp1, salp1 * sbet1); //
calp0 > 0
1262
1263 // tan(bet1) = tan(sig1) * cos(alp1)
1264 // tan(omg1) = sin(alp0) * tan(sig1) =
tan(omg1)=tan(alp1)*sin(bet1)
1265 vals.ssig1 = sbet1; somg1 = salp0 * sbet1;
1266 vals.csig1 = comg1 = calp1 * cbet1;
1267 // norm(vals.ssig1, vals.csig1);
1268 t = m.hypot(vals.ssig1, vals.csig1);
vals.ssig1 /= t; vals.csig1 /= t;
1269 // norm(somg1, comg1); -- don't need to
normalize!
1270
1271 // Enforce symmetries in the case abs(bet2)
= -bet1. Need to be careful
1272 // about this case, since this can yield
singularities in the Newton
1273 // iteration.
1274 // sin(alp2) * cos(bet2) = sin(alp0)
1275 vals.salp2 = cbet2 !== cbet1 ? salp0 /
cbet2 : salp1;
1276 // calp2 = sqrt(1 - sq(salp2))
1277 //
= sqrt(sq(calp0) - sq(sbet2)) / cbet2
1278 // and subst for calp0 and rearrange to
give (choose positive sqrt
1279 // to give alp2 in [0, pi/2]).
1280 vals.calp2 = cbet2 !== cbet1 ||
Math.abs(sbet2) !== -sbet1 ?
1281 Math.sqrt(m.sq(calp1 * cbet1) + (cbet1
< -sbet1 ?
1282 (cbet2 -
cbet1) * (cbet1 + cbet2) :
1283 (sbet1 -
sbet2) * (sbet1 + sbet2))) /
1284 cbet2 : Math.abs(calp1);
1285 // tan(bet2) = tan(sig2) * cos(alp2)
1286 // tan(omg2) = sin(alp0) * tan(sig2).
1287 vals.ssig2 = sbet2; somg2 = salp0 * sbet2;
1288 vals.csig2 = comg2 = vals.calp2 * cbet2;
1289 // norm(vals.ssig2, vals.csig2);
1290 t = m.hypot(vals.ssig2, vals.csig2);
vals.ssig2 /= t; vals.csig2 /= t;
1291 // norm(somg2, comg2); -- don't need to
normalize!
1292
1293 // sig12 = sig2 - sig1, limit to [0, pi]
1294 vals.sig12 = Math.atan2(Math.max(0,
vals.csig1 * vals.ssig2 -
1295
vals.ssig1 * vals.csig2),
1296
vals.csig1 * vals.csig2 +
1297
vals.ssig1 * vals.ssig2);
1298
1299 // omg12 = omg2 - omg1, limit to [0, pi]
1300 somg12 = Math.max(0, comg1 * somg2 - somg1
* comg2);
1301 comg12 = comg1 * comg2 + somg1 * somg2;
1302 // eta = omg12 - lam120
1303 eta = Math.atan2(somg12 * clam120 - comg12
* slam120,
1304 comg12 * clam120 + somg12
* slam120);
1305 k2 = m.sq(calp0) * this._ep2;
1306 vals.eps = k2 / (2 * (1 + Math.sqrt(1 +
k2)) + k2);
1307 this.C3f(vals.eps, C3a);
1308 B312 = (g.SinCosSeries(true, vals.ssig2,
vals.csig2, C3a) -
1309 g.SinCosSeries(true, vals.ssig1,
vals.csig1, C3a));
1310 vals.domg12 = -this.f * this.A3f(vals.eps) * salp0 *
(vals.sig12 + B312);
1311 vals.lam12 = eta + vals.domg12;
1312 if (diffp) {
1313 if (vals.calp2 === 0)
1314 vals.dlam12 = -2 * this._f1 * dn1 /
sbet1;
1315 else {
1316 nvals = this.Lengths(vals.eps,
vals.sig12,
1317 vals.ssig1,
vals.csig1, dn1,
1318 vals.ssig2,
vals.csig2, dn2,
1319 cbet1, cbet2,
g.REDUCEDLENGTH, C1a, C2a);
1320 vals.dlam12 = nvals.m12b;
1321 vals.dlam12 *= this._f1 / (vals.calp2 *
cbet2);
1322 }
1323 }
1324 return vals;
1325 };
1326
1327 /**
1328 * @summary Solve the inverse geodesic
problem.
1329 * @param {number} lat1 the latitude of the
first point in degrees.
1330 * @param {number} lon1 the longitude of the
first point in degrees.
1331 * @param {number} lat2 the latitude of the
second point in degrees.
1332 * @param {number} lon2 the longitude of the
second point in degrees.
1333 * @param {bitmask} [outmask = STANDARD]
which results to include.
1334 * @return {object} the requested results
1335 * @description The lat1, lon1, lat2, lon2,
and a12 fields of the result are
1336 *
always set. For details on the
outmask parameter, see {@tutorial
1337 *
2-interface}, "The outmask and caps parameters".
1338 */
1339 g.Geodesic.prototype.Inverse = function(lat1,
lon1, lat2, lon2, outmask) {
1340 var r, vals;
1341 if (!outmask) outmask = g.STANDARD;
1342 if (outmask === g.LONG_UNROLL) outmask |=
g.STANDARD;
1343 outmask &= g.OUT_MASK;
1344 r = this.InverseInt(lat1, lon1, lat2, lon2,
outmask);
1345 vals = r.vals;
1346 if (outmask & g.AZIMUTH) {
1347 vals.azi1 = m.atan2d(r.salp1, r.calp1);
1348 vals.azi2 = m.atan2d(r.salp2, r.calp2);
1349 }
1350 return vals;
1351 };
1352
1353 g.Geodesic.prototype.InverseInt =
function(lat1, lon1, lat2, lon2, outmask) {
1354 var vals = {},
1355 lon12, lon12s, lonsign, t, swapp,
latsign,
1356 sbet1, cbet1, sbet2, cbet2, s12x, m12x,
1357 dn1, dn2, lam12, slam12, clam12,
1358 sig12, calp1, salp1, calp2, salp2, C1a,
C2a, C3a, meridian, nvals,
1359 ssig1, csig1, ssig2, csig2, eps, omg12,
dnm,
1360 numit, salp1a, calp1a, salp1b, calp1b,
1361 tripn, tripb, v, dv, dalp1, sdalp1,
cdalp1, nsalp1,
1362 lengthmask, salp0, calp0, alp12, k2,
A4, C4a, B41, B42,
1363 somg12, comg12, domg12, dbet1, dbet2,
salp12, calp12, sdomg12, cdomg12;
1364 // Compute longitude difference (AngDiff
does this carefully). Result is
1365 // in [-180, 180] but -180 is only for
west-going geodesics. 180 is for
1366 // east-going and meridional geodesics.
1367 vals.lat1 = lat1 = m.LatFix(lat1);
vals.lat2 = lat2 = m.LatFix(lat2);
1368 // If really close to the equator, treat as
on equator.
1369 lat1 = m.AngRound(lat1);
1370 lat2 = m.AngRound(lat2);
1371 lon12 = m.AngDiff(lon1, lon2); lon12s =
lon12.e; lon12 = lon12.d;
1372 if (outmask & g.LONG_UNROLL) {
1373 vals.lon1 = lon1; vals.lon2 = (lon1 +
lon12) + lon12s;
1374 } else {
1375 vals.lon1 = m.AngNormalize(lon1);
vals.lon2 = m.AngNormalize(lon2);
1376 }
1377 // Make longitude difference positive.
1378 lonsign = m.copysign(1, lon12);
1379 lon12 *= lonsign; lon12s *= lonsign;
1380 lam12 = lon12 * m.degree;
1381 // Calculate sincos of lon12 + error (this
applies AngRound internally).
1382 t = m.sincosde(lon12, lon12s); slam12 =
t.s; clam12 = t.c;
1383 lon12s = (180 - lon12) - lon12s; // the
supplementary longitude difference
1384
1385 // Swap points so that point with higher
(abs) latitude is point 1
1386 // If one latitude is a nan, then it
becomes lat1.
1387 swapp = Math.abs(lat1) < Math.abs(lat2)
|| isNaN(lat2) ? -1 : 1;
1388 if (swapp < 0) {
1389 lonsign *= -1;
1390 [lat2, lat1] = [lat1, lat2]; //
swap(lat1, lat2);
1391 }
1392 // Make lat1 <= 0
1393 latsign = m.copysign(1, -lat1);
1394 lat1 *= latsign;
1395 lat2 *= latsign;
1396 // Now we have
1397 //
1398 //
0 <= lon12 <= 180
1399 //
-90 <= lat1 <= 0
1400 //
lat1 <= lat2 <= -lat1
1401 //
1402 // longsign, swapp, latsign register the
transformation to bring the
1403 // coordinates to this canonical form. In all cases, 1 means no change was
1404 // made.
We make these transformations so that there are few cases to
1405 // check, e.g., on verifying quadrants in
atan2. In addition, this
1406 // enforces some symmetries in the results
returned.
1407
1408 t = m.sincosd(lat1); sbet1 = this._f1 *
t.s; cbet1 = t.c;
1409 // norm(sbet1, cbet1);
1410 t = m.hypot(sbet1, cbet1); sbet1 /= t;
cbet1 /= t;
1411 // Ensure cbet1 = +epsilon at poles
1412 cbet1 = Math.max(g.tiny_, cbet1);
1413
1414 t = m.sincosd(lat2); sbet2 = this._f1 *
t.s; cbet2 = t.c;
1415 // norm(sbet2, cbet2);
1416 t = m.hypot(sbet2, cbet2); sbet2 /= t;
cbet2 /= t;
1417 // Ensure cbet2 = +epsilon at poles
1418 cbet2 = Math.max(g.tiny_, cbet2);
1419
1420 // If cbet1 < -sbet1, then cbet2 - cbet1
is a sensitive measure of the
1421 // |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1),
abs(sbet2) + sbet1 is
1422 // a better measure. This logic is used in assigning calp2 in
Lambda12.
1423 // Sometimes these quantities vanish and in
that case we force bet2 = +/-
1424 // bet1 exactly. An example where is is necessary is the
inverse problem
1425 // 48.522876735459 0 -48.52287673545898293
179.599720456223079643
1426 // which failed with Visual Studio 10
(Release and Debug)
1427
1428 if (cbet1 < -sbet1) {
1429 if (cbet2 === cbet1)
1430 sbet2 = m.copysign(sbet1, sbet2);
1431 } else {
1432 if (Math.abs(sbet2) === -sbet1)
1433 cbet2 = cbet1;
1434 }
1435
1436 dn1 = Math.sqrt(1 + this._ep2 *
m.sq(sbet1));
1437 dn2 = Math.sqrt(1 + this._ep2 *
m.sq(sbet2));
1438
1439 // index zero elements of these arrays are
unused
1440 C1a = new Array(g.nC1_ + 1);
1441 C2a = new Array(g.nC2_ + 1);
1442 C3a = new Array(g.nC3_);
1443
1444 meridian = lat1 === -90 || slam12 === 0;
1445 if (meridian) {
1446
1447 // Endpoints are on a single full
meridian, so the geodesic might
1448 // lie on a meridian.
1449
1450 calp1 = clam12; salp1 = slam12; // Head
to the target longitude
1451 calp2 = 1; salp2 = 0; // At the target we're heading north
1452
1453 // tan(bet) = tan(sig) * cos(alp)
1454 ssig1 = sbet1; csig1 = calp1 * cbet1;
1455 ssig2 = sbet2; csig2 = calp2 * cbet2;
1456
1457 // sig12 = sig2 - sig1
1458 sig12 = Math.atan2(Math.max(0, csig1 *
ssig2 - ssig1 * csig2),
1459 csig1 *
csig2 + ssig1 * ssig2);
1460 nvals = this.Lengths(this._n, sig12,
1461 ssig1, csig1, dn1,
ssig2, csig2, dn2, cbet1, cbet2,
1462 outmask | g.DISTANCE
| g.REDUCEDLENGTH,
1463 C1a, C2a);
1464 s12x = nvals.s12b;
1465 m12x = nvals.m12b;
1466 // Ignore m0
1467 if (outmask & g.GEODESICSCALE) {
1468 vals.M12 = nvals.M12;
1469 vals.M21 = nvals.M21;
1470 }
1471 // Add the check for sig12 since zero
length geodesics might yield
1472 // m12 < 0. Test case was
1473 //
1474 //
echo 20.001 0 20.001 0 | GeodSolve -i
1475 //
1476 // In fact, we will have sig12 > pi/2
for meridional geodesic
1477 // which is not a shortest path.
1478 if (sig12 < 1 || m12x >= 0) {
1479 // Need at least 2, to handle 90 0 90
180
1480 if (sig12 < 3 * g.tiny_ ||
1481 // Prevent negative s12 or m12 for
short lines
1482 (sig12 < tol0_ && (s12x
< 0 || m12x < 0)))
1483 sig12 = m12x = s12x = 0;
1484 m12x *= this._b;
1485 s12x *= this._b;
1486 vals.a12 = sig12 / m.degree;
1487 } else
1488 // m12 < 0, i.e., prolate and too
close to anti-podal
1489 meridian = false;
1490 }
1491
1492 somg12 = 2;
1493 if (!meridian &&
1494 sbet1 === 0 && // and sbet2 == 0
1495 (this.f <= 0 || lon12s >= this.f
* 180)) {
1496
1497 // Geodesic runs along equator
1498 calp1 = calp2 = 0; salp1 = salp2 = 1;
1499 s12x = this.a * lam12;
1500 sig12 = omg12 = lam12 / this._f1;
1501 m12x = this._b * Math.sin(sig12);
1502 if (outmask & g.GEODESICSCALE)
1503 vals.M12 = vals.M21 = Math.cos(sig12);
1504 vals.a12 = lon12 / this._f1;
1505
1506 } else if (!meridian) {
1507
1508 // Now point1 and point2 belong within a
hemisphere bounded by a
1509 // meridian and geodesic is neither
meridional or equatorial.
1510
1511 // Figure a starting point for Newton's
method
1512 nvals = this.InverseStart(sbet1, cbet1,
dn1, sbet2, cbet2, dn2,
1513 lam12, slam12,
clam12, C1a, C2a);
1514 sig12 = nvals.sig12;
1515 salp1 = nvals.salp1;
1516 calp1 = nvals.calp1;
1517
1518 if (sig12 >= 0) {
1519 salp2 = nvals.salp2;
1520 calp2 = nvals.calp2;
1521 // Short lines (InverseStart sets
salp2, calp2, dnm)
1522
1523 dnm = nvals.dnm;
1524 s12x = sig12 * this._b * dnm;
1525 m12x = m.sq(dnm) * this._b *
Math.sin(sig12 / dnm);
1526 if (outmask & g.GEODESICSCALE)
1527 vals.M12 = vals.M21 = Math.cos(sig12
/ dnm);
1528 vals.a12 = sig12 / m.degree;
1529 omg12 = lam12 / (this._f1 * dnm);
1530 } else {
1531
1532 // Newton's method. This is a straightforward solution of f(alp1)
=
1533 // lambda12(alp1) - lam12 = 0 with one
wrinkle. f(alp) has exactly one
1534 // root in the interval (0, pi) and its
derivative is positive at the
1535 // root. Thus f(alp) is positive for alp > alp1 and
negative for alp <
1536 // alp1. During the course of the iteration, a range
(alp1a, alp1b) is
1537 // maintained which brackets the root
and with each evaluation of
1538 // f(alp) the range is shrunk if
possible. Newton's method is
1539 // restarted whenever the derivative of
f is negative (because the new
1540 // value of alp1 is then further from
the solution) or if the new
1541 // estimate of alp1 lies outside
(0,pi); in this case, the new starting
1542 // guess is taken to be (alp1a + alp1b)
/ 2.
1543 numit = 0;
1544 // Bracketing range
1545 salp1a = g.tiny_; calp1a = 1; salp1b =
g.tiny_; calp1b = -1;
1546 for (tripn = false, tripb = false;;
++numit) {
1547 // the WGS84 test set: mean = 1.47,
sd = 1.25, max = 16
1548 // WGS84 and random input: mean =
2.85, sd = 0.60
1549 nvals = this.Lambda12(sbet1, cbet1,
dn1, sbet2, cbet2, dn2,
1550 salp1, calp1,
slam12, clam12, numit < maxit1_,
1551 C1a, C2a, C3a);
1552 v = nvals.lam12;
1553 salp2 = nvals.salp2;
1554 calp2 = nvals.calp2;
1555 sig12 = nvals.sig12;
1556 ssig1 = nvals.ssig1;
1557 csig1 = nvals.csig1;
1558 ssig2 = nvals.ssig2;
1559 csig2 = nvals.csig2;
1560 eps = nvals.eps;
1561 domg12 = nvals.domg12;
1562 dv = nvals.dlam12;
1563
1564 // Reversed test to allow escape with
NaNs
1565 // jshint -W018
1566 if (tripb || !(Math.abs(v) >=
(tripn ? 8 : 1) * tol0_) ||
1567 numit == maxit2_)
1568 break;
1569 // Update bracketing values
1570 if (v > 0 && (numit <
maxit1_ || calp1/salp1 > calp1b/salp1b)) {
1571 salp1b = salp1; calp1b = calp1;
1572 } else if (v < 0 &&
1573 (numit < maxit1_ ||
calp1/salp1 < calp1a/salp1a)) {
1574 salp1a = salp1; calp1a = calp1;
1575 }
1576 if (numit < maxit1_ && dv
> 0) {
1577 dalp1 = -v/dv;
1578 if (Math.abs(dalp1) < Math.PI) {
1579 sdalp1 = Math.sin(dalp1); cdalp1
= Math.cos(dalp1);
1580 nsalp1 = salp1 * cdalp1 + calp1 *
sdalp1;
1581 if (nsalp1 > 0) {
1582 calp1 = calp1 * cdalp1 - salp1
* sdalp1;
1583 salp1 = nsalp1;
1584 // norm(salp1, calp1);
1585 t = m.hypot(salp1, calp1);
salp1 /= t; calp1 /= t;
1586 // In some regimes we don't get
quadratic convergence because
1587 // slope -> 0. So use convergence conditions based on
epsilon
1588 // instead of sqrt(epsilon).
1589 tripn = Math.abs(v) <= 16 *
tol0_;
1590 continue;
1591 }
1592 }
1593 }
1594 // Either dv was not positive or
updated value was outside legal
1595 // range. Use the midpoint of the bracket as the next
estimate.
1596 // This mechanism is not needed for
the WGS84 ellipsoid, but it does
1597 // catch problems with more eccentric
ellipsoids. Its efficacy is
1598 // such for the WGS84 test set with
the starting guess set to alp1 =
1599 // 90deg:
1600 // the WGS84 test set: mean = 5.21,
sd = 3.93, max = 24
1601 // WGS84 and random input: mean =
4.74, sd = 0.99
1602 salp1 = (salp1a + salp1b)/2;
1603 calp1 = (calp1a + calp1b)/2;
1604 // norm(salp1, calp1);
1605 t = m.hypot(salp1, calp1); salp1 /=
t; calp1 /= t;
1606 tripn = false;
1607 tripb = (Math.abs(salp1a - salp1) +
(calp1a - calp1) < tolb_ ||
1608 Math.abs(salp1 - salp1b) +
(calp1 - calp1b) < tolb_);
1609 }
1610 lengthmask = outmask |
1611 (outmask & (g.REDUCEDLENGTH |
g.GEODESICSCALE) ?
1612 g.DISTANCE : g.NONE);
1613 nvals = this.Lengths(eps, sig12,
1614 ssig1, csig1, dn1,
ssig2, csig2, dn2,
1615 cbet1, cbet2,
1616 lengthmask, C1a,
C2a);
1617 s12x = nvals.s12b;
1618 m12x = nvals.m12b;
1619 // Ignore m0
1620 if (outmask & g.GEODESICSCALE) {
1621 vals.M12 = nvals.M12;
1622 vals.M21 = nvals.M21;
1623 }
1624 m12x *= this._b;
1625 s12x *= this._b;
1626 vals.a12 = sig12 / m.degree;
1627 if (outmask & g.AREA) {
1628 // omg12 = lam12 - domg12
1629 sdomg12 = Math.sin(domg12); cdomg12 =
Math.cos(domg12);
1630 somg12 = slam12 * cdomg12 - clam12 *
sdomg12;
1631 comg12 = clam12 * cdomg12 + slam12 *
sdomg12;
1632 }
1633 }
1634 }
1635
1636 if (outmask & g.DISTANCE)
1637 vals.s12 = 0 + s12x; //
Convert -0 to 0
1638
1639 if (outmask & g.REDUCEDLENGTH)
1640 vals.m12 = 0 + m12x; // Convert -0 to 0
1641
1642 if (outmask & g.AREA) {
1643 // From Lambda12: sin(alp1) * cos(bet1) =
sin(alp0)
1644 salp0 = salp1 * cbet1;
1645 calp0 = m.hypot(calp1, salp1 * sbet1); //
calp0 > 0
1646 if (calp0 !== 0 && salp0 !== 0) {
1647 // From Lambda12: tan(bet) = tan(sig) *
cos(alp)
1648 ssig1 = sbet1; csig1 = calp1 * cbet1;
1649 ssig2 = sbet2; csig2 = calp2 * cbet2;
1650 k2 = m.sq(calp0) * this._ep2;
1651 eps = k2 / (2 * (1 + Math.sqrt(1 + k2))
+ k2);
1652 // Multiplier = a^2 * e^2 * cos(alpha0)
* sin(alpha0).
1653 A4 = m.sq(this.a) * calp0 * salp0 *
this._e2;
1654 // norm(ssig1, csig1);
1655 t = m.hypot(ssig1, csig1); ssig1 /= t;
csig1 /= t;
1656 // norm(ssig2, csig2);
1657 t = m.hypot(ssig2, csig2); ssig2 /= t;
csig2 /= t;
1658 C4a = new Array(g.nC4_);
1659 this.C4f(eps, C4a);
1660 B41 = g.SinCosSeries(false, ssig1,
csig1, C4a);
1661 B42 = g.SinCosSeries(false, ssig2,
csig2, C4a);
1662 vals.S12 = A4 * (B42 - B41);
1663 } else
1664 // Avoid problems with indeterminate
sig1, sig2 on equator
1665 vals.S12 = 0;
1666 if (!meridian && somg12 == 2) {
1667 somg12 = Math.sin(omg12); comg12 =
Math.cos(omg12);
1668 }
1669 if (!meridian &&
1670 comg12 > -0.7071 && // Long difference not too big
1671 sbet2 - sbet1 < 1.75) { // Lat
difference not too big
1672 // Use tan(Gamma/2) = tan(omg12/2)
1673 // *
(tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
1674 // with tan(x/2) = sin(x)/(1+cos(x))
1675 domg12 = 1 + comg12; dbet1 = 1 + cbet1;
dbet2 = 1 + cbet2;
1676 alp12 = 2 * Math.atan2( somg12 *
(sbet1*dbet2 + sbet2*dbet1),
1677 domg12 *
(sbet1*sbet2 + dbet1*dbet2) );
1678 } else {
1679 // alp12 = alp2 - alp1, used in atan2
so no need to normalize
1680 salp12 = salp2 * calp1 - calp2 * salp1;
1681 calp12 = calp2 * calp1 + salp2 * salp1;
1682 // The right thing appears to happen if
alp1 = +/-180 and alp2 = 0, viz
1683 // salp12 = -0 and alp12 = -180. However this depends on the sign
1684 // being attached to 0 correctly. The following ensures the correct
1685 // behavior.
1686 if (salp12 === 0 && calp12 <
0) {
1687 salp12 = g.tiny_ * calp1;
1688 calp12 = -1;
1689 }
1690 alp12 = Math.atan2(salp12, calp12);
1691 }
1692 vals.S12 += this._c2 * alp12;
1693 vals.S12 *= swapp * lonsign * latsign;
1694 // Convert -0 to 0
1695 vals.S12 += 0;
1696 }
1697
1698 // Convert calp, salp to azimuth accounting
for lonsign, swapp, latsign.
1699 if (swapp < 0) {
1700 [salp2, salp1] = [salp1, salp2]; //
swap(salp1, salp2);
1701 [calp2, calp1] = [calp1, calp2]; //
swap(calp1, calp2);
1702 if (outmask & g.GEODESICSCALE) {
1703 [vals.M21, vals.M12] = [vals.M12,
vals.M21]; //swap(vals.M12, vals.M21);
1704 }
1705 }
1706
1707 salp1 *= swapp * lonsign; calp1 *= swapp *
latsign;
1708 salp2 *= swapp * lonsign; calp2 *= swapp *
latsign;
1709
1710 return {vals: vals,
1711 salp1: salp1, calp1: calp1,
1712 salp2: salp2, calp2: calp2};
1713 };
1714
1715 /**
1716 * @summary Solve the general direct geodesic
problem.
1717 * @param {number} lat1 the latitude of the
first point in degrees.
1718 * @param {number} lon1 the longitude of the
first point in degrees.
1719 * @param {number} azi1 the azimuth at the
first point in degrees.
1720 * @param {bool} arcmode is the next
parameter an arc length?
1721 * @param {number} s12_a12 the (arcmode ? arc
length : distance) from the
1722 *
first point to the second in (arcmode ? degrees : meters).
1723 * @param {bitmask} [outmask = STANDARD]
which results to include.
1724 * @return {object} the requested results.
1725 * @description The lat1, lon1, azi1, and a12
fields of the result are always
1726 *
set; s12 is included if arcmode is false. For details on the outmask
1727 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1728 *
parameters".
1729 */
1730 g.Geodesic.prototype.GenDirect =
function(lat1, lon1, azi1,
1731
arcmode, s12_a12, outmask) {
1732 var line;
1733 if (!outmask) outmask = g.STANDARD;
1734 else if (outmask === g.LONG_UNROLL) outmask
|= g.STANDARD;
1735 // Automatically supply DISTANCE_IN if
necessary
1736 if (!arcmode) outmask |= g.DISTANCE_IN;
1737 line = new l.GeodesicLine(this, lat1, lon1,
azi1, outmask);
1738 return line.GenPosition(arcmode, s12_a12,
outmask);
1739 };
1740
1741 /**
1742 * @summary Solve the direct geodesic
problem.
1743 * @param {number} lat1 the latitude of the
first point in degrees.
1744 * @param {number} lon1 the longitude of the
first point in degrees.
1745 * @param {number} azi1 the azimuth at the
first point in degrees.
1746 * @param {number} s12 the distance from the
first point to the second in
1747 *
meters.
1748 * @param {bitmask} [outmask = STANDARD]
which results to include.
1749 * @return {object} the requested results.
1750 * @description The lat1, lon1, azi1, s12,
and a12 fields of the result are
1751 *
always set. For details on the
outmask parameter, see {@tutorial
1752 *
2-interface}, "The outmask and caps parameters".
1753 */
1754 g.Geodesic.prototype.Direct = function(lat1,
lon1, azi1, s12, outmask) {
1755 return this.GenDirect(lat1, lon1, azi1,
false, s12, outmask);
1756 };
1757
1758 /**
1759 * @summary Solve the direct geodesic problem
with arc length.
1760 * @param {number} lat1 the latitude of the
first point in degrees.
1761 * @param {number} lon1 the longitude of the
first point in degrees.
1762 * @param {number} azi1 the azimuth at the
first point in degrees.
1763 * @param {number} a12 the arc length from
the first point to the second in
1764 *
degrees.
1765 * @param {bitmask} [outmask = STANDARD]
which results to include.
1766 * @return {object} the requested results.
1767 * @description The lat1, lon1, azi1, and a12
fields of the result are
1768 *
always set. For details on the
outmask parameter, see {@tutorial
1769 *
2-interface}, "The outmask and caps parameters".
1770 */
1771 g.Geodesic.prototype.ArcDirect =
function(lat1, lon1, azi1, a12, outmask) {
1772 return this.GenDirect(lat1, lon1, azi1,
true, a12, outmask);
1773 };
1774
1775 /**
1776 * @summary Create a {@link
module:geodesic/GeodesicLine.GeodesicLine
1777 * GeodesicLine}
object.
1778 * @param {number} lat1 the latitude of the
first point in degrees.
1779 * @param {number} lon1 the longitude of the
first point in degrees.
1780 * @param {number} azi1 the azimuth at the
first point in degrees.
1781 *
degrees.
1782 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1783 *
include.
1784 * @return {object} the
1785 *
{@link module:geodesic/GeodesicLine.GeodesicLine
1786 *
GeodesicLine} object
1787 * @description For details on the caps
parameter, see {@tutorial
1788 *
2-interface}, "The outmask and caps parameters".
1789 */
1790 g.Geodesic.prototype.Line = function(lat1,
lon1, azi1, caps) {
1791 return new l.GeodesicLine(this, lat1, lon1,
azi1, caps);
1792 };
1793
1794 /**
1795 * @summary Define a {@link
module:geodesic/GeodesicLine.GeodesicLine
1796 *
GeodesicLine} in terms of the direct geodesic problem specified in terms
1797 * of
distance.
1798 * @param {number} lat1 the latitude of the
first point in degrees.
1799 * @param {number} lon1 the longitude of the
first point in degrees.
1800 * @param {number} azi1 the azimuth at the
first point in degrees.
1801 *
degrees.
1802 * @param {number} s12 the distance between
point 1 and point 2 (meters); it
1803 *
can be negative.
1804 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1805 *
include.
1806 * @return {object} the
1807 *
{@link module:geodesic/GeodesicLine.GeodesicLine
1808 *
GeodesicLine} object
1809 * @description This function sets point 3 of
the GeodesicLine to correspond
1810 * to
point 2 of the direct geodesic problem.
For details on the caps
1811 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1812 *
parameters".
1813 */
1814 g.Geodesic.prototype.DirectLine =
function(lat1, lon1, azi1, s12, caps) {
1815 return this.GenDirectLine(lat1, lon1, azi1,
false, s12, caps);
1816 };
1817
1818 /**
1819 * @summary Define a {@link
module:geodesic/GeodesicLine.GeodesicLine
1820 *
GeodesicLine} in terms of the direct geodesic problem specified in terms
1821 * of
arc length.
1822 * @param {number} lat1 the latitude of the
first point in degrees.
1823 * @param {number} lon1 the longitude of the
first point in degrees.
1824 * @param {number} azi1 the azimuth at the
first point in degrees.
1825 *
degrees.
1826 * @param {number} a12 the arc length between
point 1 and point 2 (degrees);
1827 * it
can be negative.
1828 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1829 *
include.
1830 * @return {object} the
1831 *
{@link module:geodesic/GeodesicLine.GeodesicLine
1832 *
GeodesicLine} object
1833 * @description This function sets point 3 of
the GeodesicLine to correspond
1834 * to
point 2 of the direct geodesic problem.
For details on the caps
1835 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1836 *
parameters".
1837 */
1838 g.Geodesic.prototype.ArcDirectLine =
function(lat1, lon1, azi1, a12, caps) {
1839 return this.GenDirectLine(lat1, lon1, azi1,
true, a12, caps);
1840 };
1841
1842 /**
1843 * @summary Define a {@link
module:geodesic/GeodesicLine.GeodesicLine
1844 *
GeodesicLine} in terms of the direct geodesic problem specified in terms
1845 * of
either distance or arc length.
1846 * @param {number} lat1 the latitude of the
first point in degrees.
1847 * @param {number} lon1 the longitude of the
first point in degrees.
1848 * @param {number} azi1 the azimuth at the
first point in degrees.
1849 *
degrees.
1850 * @param {bool} arcmode boolean flag
determining the meaning of the
1851 *
s12_a12.
1852 * @param {number} s12_a12 if arcmode is
false, this is the distance between
1853 *
point 1 and point 2 (meters); otherwise it is the arc length between
1854 *
point 1 and point 2 (degrees); it can be negative.
1855 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1856 *
include.
1857 * @return {object} the
1858 *
{@link module:geodesic/GeodesicLine.GeodesicLine
1859 *
GeodesicLine} object
1860 * @description This function sets point 3 of
the GeodesicLine to correspond
1861 * to
point 2 of the direct geodesic problem.
For details on the caps
1862 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1863 *
parameters".
1864 */
1865 g.Geodesic.prototype.GenDirectLine =
function(lat1, lon1, azi1,
1866
arcmode, s12_a12, caps) {
1867 var t;
1868 if (!caps) caps = g.STANDARD |
g.DISTANCE_IN;
1869 // Automatically supply DISTANCE_IN if
necessary
1870 if (!arcmode) caps |= g.DISTANCE_IN;
1871 t = new l.GeodesicLine(this, lat1, lon1,
azi1, caps);
1872 t.GenSetDistance(arcmode, s12_a12);
1873 return t;
1874 };
1875
1876 /**
1877 * @summary Define a {@link
module:geodesic/GeodesicLine.GeodesicLine
1878 *
GeodesicLine} in terms of the inverse geodesic problem.
1879 * @param {number} lat1 the latitude of the
first point in degrees.
1880 * @param {number} lon1 the longitude of the
first point in degrees.
1881 * @param {number} lat2 the latitude of the
second point in degrees.
1882 * @param {number} lon2 the longitude of the
second point in degrees.
1883 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1884 *
include.
1885 * @return {object} the
1886 *
{@link module:geodesic/GeodesicLine.GeodesicLine
1887 * GeodesicLine}
object
1888 * @description This function sets point 3 of
the GeodesicLine to correspond
1889 * to
point 2 of the inverse geodesic problem.
For details on the caps
1890 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1891 *
parameters".
1892 */
1893 g.Geodesic.prototype.InverseLine =
function(lat1, lon1, lat2, lon2, caps) {
1894 var r, t, azi1;
1895 if (!caps) caps = g.STANDARD |
g.DISTANCE_IN;
1896 r = this.InverseInt(lat1, lon1, lat2, lon2,
g.ARC);
1897 azi1 = m.atan2d(r.salp1, r.calp1);
1898 // Ensure that a12 can be converted to a
distance
1899 if (caps & (g.OUT_MASK &
g.DISTANCE_IN)) caps |= g.DISTANCE;
1900 t = new l.GeodesicLine(this, lat1, lon1,
azi1, caps, r.salp1, r.calp1);
1901 t.SetArc(r.vals.a12);
1902 return t;
1903 };
1904
1905 /**
1906 * @summary Create a {@link
module:geodesic/PolygonArea.PolygonArea
1907 *
PolygonArea} object.
1908 * @param {bool} [polyline = false] if true
the new PolygonArea object
1909 *
describes a polyline instead of a polygon.
1910 * @return {object} the
1911 *
{@link module:geodesic/PolygonArea.PolygonArea
1912 *
PolygonArea} object
1913 */
1914 g.Geodesic.prototype.Polygon =
function(polyline) {
1915 return new p.PolygonArea(this, polyline);
1916 };
1917
1918 /**
1919 * @summary a {@link
module:geodesic/Geodesic.Geodesic Geodesic} object
1920 *
initialized for the WGS84 ellipsoid.
1921 * @constant {object}
1922 */
1923 g.WGS84 = new g.Geodesic(c.WGS84.a,
c.WGS84.f);
1924})(geodesic.Geodesic,
geodesic.GeodesicLine,
1925 geodesic.PolygonArea, geodesic.Math,
geodesic.Constants);
1926/*
1927 * GeodesicLine.js
1928 * Transcription of
GeodesicLine.[ch]pp into JavaScript.
1929 *
1930 * See the documentation
for the C++ class. The conversion is a
literal
1931 * conversion from C++.
1932 *
1933 * The algorithms are
derived in
1934 *
1935 * Charles F. F. Karney,
1936 * Algorithms for geodesics, J. Geodesy 87,
43-55 (2013);
1937 * https://doi.org/10.1007/s00190-012-0578-z
1938 * Addenda:
https://geographiclib.sourceforge.io/geod-addenda.html
1939 *
1940 * Copyright (c) Charles
Karney (2011-2022) <karney@alum.mit.edu> and licensed
1941 * under the MIT/X11
License. For more information, see
1942 *
https://geographiclib.sourceforge.io/
1943 */
1944
1945// Load AFTER
geodesic/Math.js, geodesic/Geodesic.js
1946
1947(function(
1948 g,
1949 /**
1950 * @exports geodesic/GeodesicLine
1951 * @description Solve geodesic problems on a
single geodesic line via the
1952 *
{@link module:geodesic/GeodesicLine.GeodesicLine GeodesicLine}
1953 *
class.
1954 */
1955 l, m) {
1956 "use strict";
1957
1958 /**
1959 * @class
1960 * @property {number} a the equatorial radius
(meters).
1961 * @property {number} f the flattening.
1962 * @property {number} lat1 the initial
latitude (degrees).
1963 * @property {number} lon1 the initial
longitude (degrees).
1964 * @property {number} azi1 the initial
azimuth (degrees).
1965 * @property {number} salp1 the sine of the
azimuth at the first point.
1966 * @property {number} calp1 the cosine the
azimuth at the first point.
1967 * @property {number} s13 the distance to
point 3 (meters).
1968 * @property {number} a13 the arc length to
point 3 (degrees).
1969 * @property {bitmask} caps the capabilities
of the object.
1970 * @summary Initialize a GeodesicLine
object. For details on the caps
1971 *
parameter, see {@tutorial 2-interface}, "The outmask and caps
1972 *
parameters".
1973 * @classdesc Performs geodesic calculations
along a given geodesic line.
1974 *
This object is usually instantiated by
1975 *
{@link module:geodesic/Geodesic.Geodesic#Line Geodesic.Line}.
1976 *
The methods
1977 *
{@link module:geodesic/Geodesic.Geodesic#DirectLine
1978 *
Geodesic.DirectLine} and
1979 *
{@link module:geodesic/Geodesic.Geodesic#InverseLine
1980 * Geodesic.InverseLine}
set in addition the position of a reference point
1981 * 3.
1982 * @param {object} geod a {@link
module:geodesic/Geodesic.Geodesic
1983 *
Geodesic} object.
1984 * @param {number} lat1 the latitude of the
first point in degrees.
1985 * @param {number} lon1 the longitude of the
first point in degrees.
1986 * @param {number} azi1 the azimuth at the
first point in degrees.
1987 * @param {bitmask} [caps = STANDARD |
DISTANCE_IN] which capabilities to
1988 *
include; LATITUDE | AZIMUTH are always included.
1989 */
1990 l.GeodesicLine = function(geod, lat1, lon1,
azi1, caps, salp1, calp1) {
1991 var t, cbet1, sbet1, eps, s, c;
1992 if (!caps) caps = g.STANDARD |
g.DISTANCE_IN;
1993
1994 this.a = geod.a;
1995 this.f = geod.f;
1996 this._b = geod._b;
1997 this._c2 = geod._c2;
1998 this._f1 = geod._f1;
1999 this.caps = caps | g.LATITUDE | g.AZIMUTH |
g.LONG_UNROLL;
2000
2001 this.lat1 = m.LatFix(lat1);
2002 this.lon1 = lon1;
2003 if (typeof salp1 === 'undefined' || typeof
calp1 === 'undefined') {
2004 this.azi1 = m.AngNormalize(azi1);
2005 t = m.sincosd(m.AngRound(this.azi1));
this.salp1 = t.s; this.calp1 = t.c;
2006 } else {
2007 this.azi1 = azi1; this.salp1 = salp1;
this.calp1 = calp1;
2008 }
2009 t = m.sincosd(m.AngRound(this.lat1)); sbet1
= this._f1 * t.s; cbet1 = t.c;
2010 // norm(sbet1, cbet1);
2011 t = m.hypot(sbet1, cbet1); sbet1 /= t;
cbet1 /= t;
2012 // Ensure cbet1 = +epsilon at poles
2013 cbet1 = Math.max(g.tiny_, cbet1);
2014 this._dn1 = Math.sqrt(1 + geod._ep2 *
m.sq(sbet1));
2015
2016 // Evaluate alp0 from sin(alp1) * cos(bet1)
= sin(alp0),
2017 this._salp0 = this.salp1 * cbet1; // alp0
in [0, pi/2 - |bet1|]
2018 // Alt: calp0 = hypot(sbet1, calp1 *
cbet1). The following
2019 // is slightly better (consider the case
salp1 = 0).
2020 this._calp0 = m.hypot(this.calp1,
this.salp1 * sbet1);
2021 // Evaluate sig with tan(bet1) = tan(sig1)
* cos(alp1).
2022 // sig = 0 is nearest northward crossing of
equator.
2023 // With bet1 = 0, alp1 = pi/2, we have sig1
= 0 (equatorial line).
2024 // With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
2025 // With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
2026 // Evaluate omg1 with tan(omg1) = sin(alp0)
* tan(sig1).
2027 // With alp0 in (0, pi/2], quadrants for
sig and omg coincide.
2028 // No atan2(0,0) ambiguity at poles since
cbet1 = +epsilon.
2029 // With alp0 = 0, omg1 = 0 for alp1 = 0,
omg1 = pi for alp1 = pi.
2030 this._ssig1 = sbet1; this._somg1 =
this._salp0 * sbet1;
2031 this._csig1 = this._comg1 =
2032 sbet1 !== 0 || this.calp1 !== 0 ? cbet1 *
this.calp1 : 1;
2033 // norm(this._ssig1, this._csig1); // sig1
in (-pi, pi]
2034 t = m.hypot(this._ssig1, this._csig1);
2035 this._ssig1 /= t; this._csig1 /= t;
2036 // norm(this._somg1, this._comg1); -- don't
need to normalize!
2037
2038 this._k2 = m.sq(this._calp0) * geod._ep2;
2039 eps = this._k2 / (2 * (1 + Math.sqrt(1 +
this._k2)) + this._k2);
2040
2041 if (this.caps & g.CAP_C1) {
2042 this._A1m1 = g.A1m1f(eps);
2043 this._C1a = new Array(g.nC1_ + 1);
2044 g.C1f(eps, this._C1a);
2045 this._B11 = g.SinCosSeries(true,
this._ssig1, this._csig1, this._C1a);
2046 s = Math.sin(this._B11); c =
Math.cos(this._B11);
2047 // tau1 = sig1 + B11
2048 this._stau1 = this._ssig1 * c +
this._csig1 * s;
2049 this._ctau1 = this._csig1 * c -
this._ssig1 * s;
2050 // Not necessary because C1pa reverts C1a
2051 //
_B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa);
2052 }
2053
2054 if (this.caps & g.CAP_C1p) {
2055 this._C1pa = new Array(g.nC1p_ + 1);
2056 g.C1pf(eps, this._C1pa);
2057 }
2058
2059 if (this.caps & g.CAP_C2) {
2060 this._A2m1 = g.A2m1f(eps);
2061 this._C2a = new Array(g.nC2_ + 1);
2062 g.C2f(eps, this._C2a);
2063 this._B21 = g.SinCosSeries(true,
this._ssig1, this._csig1, this._C2a);
2064 }
2065
2066 if (this.caps & g.CAP_C3) {
2067 this._C3a = new Array(g.nC3_);
2068 geod.C3f(eps, this._C3a);
2069 this._A3c = -this.f * this._salp0 *
geod.A3f(eps);
2070 this._B31 = g.SinCosSeries(true, this._ssig1,
this._csig1, this._C3a);
2071 }
2072
2073 if (this.caps & g.CAP_C4) {
2074 this._C4a = new Array(g.nC4_); // all the
elements of _C4a are used
2075 geod.C4f(eps, this._C4a);
2076 // Multiplier = a^2 * e^2 * cos(alpha0) *
sin(alpha0)
2077 this._A4 = m.sq(this.a) * this._calp0 *
this._salp0 * geod._e2;
2078 this._B41 = g.SinCosSeries(false,
this._ssig1, this._csig1, this._C4a);
2079 }
2080
2081 this.a13 = this.s13 = NaN;
2082 };
2083
2084 /**
2085 * @summary Find the position on the line
(general case).
2086 * @param {bool} arcmode is the next
parameter an arc length?
2087 * @param {number} s12_a12 the (arcmode ? arc
length : distance) from the
2088 *
first point to the second in (arcmode ? degrees : meters).
2089 * @param {bitmask} [outmask = STANDARD]
which results to include; this is
2090 *
subject to the capabilities of the object.
2091 * @return {object} the requested results.
2092 * @description The lat1, lon1, azi1, and a12
fields of the result are
2093 *
always set; s12 is included if arcmode is false. For details on the
2094 *
outmask parameter, see {@tutorial 2-interface}, "The outmask and
caps
2095 *
parameters".
2096 */
2097 l.GeodesicLine.prototype.GenPosition =
function(arcmode, s12_a12,
2098
outmask) {
2099 var vals = {},
2100 sig12, ssig12, csig12, B12, AB1, ssig2,
csig2, tau12, s, c, serr,
2101 omg12, lam12, lon12, E, sbet2, cbet2,
somg2, comg2, salp2, calp2, dn2,
2102 B22, AB2, J12, t, B42, salp12, calp12;
2103 if (!outmask) outmask = g.STANDARD;
2104 else if (outmask === g.LONG_UNROLL) outmask
|= g.STANDARD;
2105 outmask &= this.caps & g.OUT_MASK;
2106 vals.lat1 = this.lat1; vals.azi1 =
this.azi1;
2107 vals.lon1 = outmask & g.LONG_UNROLL ?
2108 this.lon1 : m.AngNormalize(this.lon1);
2109 if (arcmode)
2110 vals.a12 = s12_a12;
2111 else
2112 vals.s12 = s12_a12;
2113 if (!( arcmode || (this.caps &
g.DISTANCE_IN & g.OUT_MASK) )) {
2114 // Uninitialized or impossible distance
calculation requested
2115 vals.a12 = NaN;
2116 return vals;
2117 }
2118
2119 // Avoid warning about uninitialized B12.
2120 B12 = 0; AB1 = 0;
2121 if (arcmode) {
2122 // Interpret s12_a12 as spherical arc
length
2123 sig12 = s12_a12 * m.degree;
2124 t = m.sincosd(s12_a12); ssig12 = t.s;
csig12 = t.c;
2125 } else {
2126 // Interpret s12_a12 as distance
2127 tau12 = s12_a12 / (this._b * (1 +
this._A1m1));
2128 s = Math.sin(tau12);
2129 c = Math.cos(tau12);
2130 // tau2 = tau1 + tau12
2131 B12 = -g.SinCosSeries(true,
2132 this._stau1 * c +
this._ctau1 * s,
2133 this._ctau1 * c -
this._stau1 * s,
2134 this._C1pa);
2135 sig12 = tau12 - (B12 - this._B11);
2136 ssig12 = Math.sin(sig12); csig12 =
Math.cos(sig12);
2137 if (Math.abs(this.f) > 0.01) {
2138 // Reverted distance series is
inaccurate for |f| > 1/100, so correct
2139 // sig12 with 1 Newton iteration. The following table shows the
2140 // approximate maximum error for a =
WGS_a() and various f relative to
2141 // GeodesicExact.
2142 //
erri = the error in the inverse solution (nm)
2143 //
errd = the error in the direct solution (series only) (nm)
2144 //
errda = the error in the direct solution
2145 // (series + 1 Newton) (nm)
2146 //
2147 //
f erri errd errda
2148 //
-1/5 12e6 1.2e9 69e6
2149 //
-1/10 123e3 12e6 765e3
2150 //
-1/20 1110 108e3 7155
2151 //
-1/50 18.63 200.9 27.12
2152 //
-1/100 18.63 23.78 23.37
2153 //
-1/150 18.63 21.05 20.26
2154 //
1/150 22.35 24.73 25.83
2155 //
1/100 22.35 25.03 25.31
2156 //
1/50 29.80 231.9 30.44
2157 //
1/20 5376 146e3 10e3
2158 //
1/10 829e3 22e6 1.5e6
2159 //
1/5 157e6 3.8e9 280e6
2160 ssig2 = this._ssig1 * csig12 +
this._csig1 * ssig12;
2161 csig2 = this._csig1 * csig12 -
this._ssig1 * ssig12;
2162 B12 = g.SinCosSeries(true, ssig2,
csig2, this._C1a);
2163 serr = (1 + this._A1m1) * (sig12 + (B12
- this._B11)) -
2164 s12_a12 / this._b;
2165 sig12 = sig12 - serr / Math.sqrt(1 +
this._k2 * m.sq(ssig2));
2166 ssig12 = Math.sin(sig12); csig12 =
Math.cos(sig12);
2167 // Update B12 below
2168 }
2169 }
2170
2171 // sig2 = sig1 + sig12
2172 ssig2 = this._ssig1 * csig12 + this._csig1
* ssig12;
2173 csig2 = this._csig1 * csig12 - this._ssig1
* ssig12;
2174 dn2 = Math.sqrt(1 + this._k2 *
m.sq(ssig2));
2175 if (outmask & (g.DISTANCE |
g.REDUCEDLENGTH | g.GEODESICSCALE)) {
2176 if (arcmode || Math.abs(this.f) >
0.01)
2177 B12 = g.SinCosSeries(true, ssig2,
csig2, this._C1a);
2178 AB1 = (1 + this._A1m1) * (B12 -
this._B11);
2179 }
2180 // sin(bet2) = cos(alp0) * sin(sig2)
2181 sbet2 = this._calp0 * ssig2;
2182 // Alt: cbet2 = hypot(csig2, salp0 *
ssig2);
2183 cbet2 = m.hypot(this._salp0, this._calp0 *
csig2);
2184 if (cbet2 === 0)
2185 // I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
2186 cbet2 = csig2 = g.tiny_;
2187 // tan(alp0) = cos(sig2)*tan(alp2)
2188 salp2 = this._salp0; calp2 = this._calp0 *
csig2; // No need to normalize
2189
2190 if (arcmode && (outmask &
g.DISTANCE))
2191 vals.s12 = this._b * ((1 + this._A1m1) *
sig12 + AB1);
2192
2193 if (outmask & g.LONGITUDE) {
2194 // tan(omg2) = sin(alp0) * tan(sig2)
2195 somg2 = this._salp0 * ssig2; comg2 =
csig2; // No need to normalize
2196 E = m.copysign(1, this._salp0);
2197 // omg12 = omg2 - omg1
2198 omg12 = outmask & g.LONG_UNROLL ?
2199 E * (sig12 -
2200 (Math.atan2(ssig2, csig2) -
2201 Math.atan2(this._ssig1,
this._csig1)) +
2202 (Math.atan2(E * somg2, comg2) -
2203 Math.atan2(E * this._somg1,
this._comg1))) :
2204 Math.atan2(somg2 * this._comg1 - comg2
* this._somg1,
2205 comg2 * this._comg1 +
somg2 * this._somg1);
2206 lam12 = omg12 + this._A3c *
2207 ( sig12 + (g.SinCosSeries(true, ssig2,
csig2, this._C3a) -
2208 this._B31));
2209 lon12 = lam12 / m.degree;
2210 vals.lon2 = outmask & g.LONG_UNROLL ?
this.lon1 + lon12 :
2211
m.AngNormalize(m.AngNormalize(this.lon1) + m.AngNormalize(lon12));
2212 }
2213
2214 if (outmask & g.LATITUDE)
2215 vals.lat2 = m.atan2d(sbet2, this._f1 *
cbet2);
2216
2217 if (outmask & g.AZIMUTH)
2218 vals.azi2 = m.atan2d(salp2, calp2);
2219
2220 if (outmask & (g.REDUCEDLENGTH |
g.GEODESICSCALE)) {
2221 B22 = g.SinCosSeries(true, ssig2, csig2,
this._C2a);
2222 AB2 = (1 + this._A2m1) * (B22 -
this._B21);
2223 J12 = (this._A1m1 - this._A2m1) * sig12 +
(AB1 - AB2);
2224 if (outmask & g.REDUCEDLENGTH)
2225 // Add parens around (_csig1 * ssig2)
and (_ssig1 * csig2) to ensure
2226 // accurate cancellation in the case of
coincident points.
2227 vals.m12 = this._b * (( dn2 * (this._csig1 * ssig2) -
2228 this._dn1 *
(this._ssig1 * csig2)) -
2229 this._csig1 *
csig2 * J12);
2230 if (outmask & g.GEODESICSCALE) {
2231 t = this._k2 * (ssig2 - this._ssig1) *
(ssig2 + this._ssig1) /
2232 (this._dn1 + dn2);
2233 vals.M12 = csig12 +
2234 (t * ssig2 - csig2 * J12) *
this._ssig1 / this._dn1;
2235 vals.M21 = csig12 -
2236 (t * this._ssig1 - this._csig1 * J12)
* ssig2 / dn2;
2237 }
2238 }
2239
2240 if (outmask & g.AREA) {
2241 B42 = g.SinCosSeries(false, ssig2, csig2,
this._C4a);
2242 if (this._calp0 === 0 || this._salp0 ===
0) {
2243 // alp12 = alp2 - alp1, used in atan2
so no need to normalize
2244 salp12 = salp2 * this.calp1 - calp2 *
this.salp1;
2245 calp12 = calp2 * this.calp1 + salp2 *
this.salp1;
2246 } else {
2247 // tan(alp) = tan(alp0) * sec(sig)
2248 // tan(alp2-alp1) = (tan(alp2)
-tan(alp1)) / (tan(alp2)*tan(alp1)+1)
2249 // = calp0 * salp0 * (csig1-csig2) /
(salp0^2 + calp0^2 * csig1*csig2)
2250 // If csig12 > 0, write
2251 //
csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
2252 // else
2253 //
csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
2254 // No need to normalize
2255 salp12 = this._calp0 * this._salp0 *
2256 (csig12 <= 0 ? this._csig1 * (1 -
csig12) + ssig12 * this._ssig1 :
2257 ssig12 * (this._csig1 * ssig12 / (1
+ csig12) + this._ssig1));
2258 calp12 = m.sq(this._salp0) +
m.sq(this._calp0) * this._csig1 * csig2;
2259 }
2260 vals.S12 = this._c2 * Math.atan2(salp12,
calp12) +
2261 this._A4 * (B42 - this._B41);
2262 }
2263
2264 if (!arcmode)
2265 vals.a12 = sig12 / m.degree;
2266 return vals;
2267 };
2268
2269 /**
2270 * @summary Find the position on the line
given s12.
2271 * @param {number} s12 the distance from the
first point to the second in
2272 *
meters.
2273 * @param {bitmask} [outmask = STANDARD]
which results to include; this is
2274 *
subject to the capabilities of the object.
2275 * @return {object} the requested results.
2276 * @description The lat1, lon1, azi1, s12,
and a12 fields of the result are
2277 *
always set; s12 is included if arcmode is false. For details on the
2278 *
outmask parameter, see {@tutorial 2-interface}, "The outmask and
caps
2279 *
parameters".
2280 */
2281 l.GeodesicLine.prototype.Position =
function(s12, outmask) {
2282 return this.GenPosition(false, s12,
outmask);
2283 };
2284
2285 /**
2286 * @summary Find the position on the line
given a12.
2287 * @param {number} a12 the arc length from
the first point to the second in
2288 *
degrees.
2289 * @param {bitmask} [outmask = STANDARD]
which results to include; this is
2290 *
subject to the capabilities of the object.
2291 * @return {object} the requested results.
2292 * @description The lat1, lon1, azi1, and a12
fields of the result are
2293 *
always set. For details on the
outmask parameter, see {@tutorial
2294 *
2-interface}, "The outmask and caps parameters".
2295 */
2296 l.GeodesicLine.prototype.ArcPosition =
function(a12, outmask) {
2297 return this.GenPosition(true, a12,
outmask);
2298 };
2299
2300 /**
2301 * @summary Specify position of point 3 in
terms of either distance or arc
2302 *
length.
2303 * @param {bool} arcmode boolean flag
determining the meaning of the second
2304 *
parameter; if arcmode is false, then the GeodesicLine object must have
2305 *
been constructed with caps |= DISTANCE_IN.
2306 * @param {number} s13_a13 if arcmode is
false, this is the distance from
2307 *
point 1 to point 3 (meters); otherwise it is the arc length from
2308 *
point 1 to point 3 (degrees); it can be negative.
2309 */
2310 l.GeodesicLine.prototype.GenSetDistance =
function(arcmode, s13_a13) {
2311 if (arcmode)
2312 this.SetArc(s13_a13);
2313 else
2314 this.SetDistance(s13_a13);
2315 };
2316
2317 /**
2318 * @summary Specify position of point 3 in
terms distance.
2319 * @param {number} s13 the distance from
point 1 to point 3 (meters); it
2320 *
can be negative.
2321 */
2322 l.GeodesicLine.prototype.SetDistance =
function(s13) {
2323 var r;
2324 this.s13 = s13;
2325 r = this.GenPosition(false, this.s13,
g.ARC);
2326 this.a13 = 0 + r.a12; // the 0+ converts undefined into NaN
2327 };
2328
2329 /**
2330 * @summary Specify position of point 3 in
terms of arc length.
2331 * @param {number} a13 the arc length from
point 1 to point 3 (degrees);
2332 * it
can be negative.
2333 */
2334 l.GeodesicLine.prototype.SetArc =
function(a13) {
2335 var r;
2336 this.a13 = a13;
2337 r = this.GenPosition(true, this.a13,
g.DISTANCE);
2338 this.s13 = 0 + r.s12; // the 0+ converts undefined into NaN
2339 };
2340
2341})(geodesic.Geodesic,
geodesic.GeodesicLine, geodesic.Math);
2342/*
2343 * PolygonArea.js
2344 * Transcription of
PolygonArea.[ch]pp into JavaScript.
2345 *
2346 * See the documentation
for the C++ class. The conversion is a
literal
2347 * conversion from C++.
2348 *
2349 * The algorithms are
derived in
2350 *
2351 * Charles F. F. Karney,
2352 * Algorithms for geodesics, J. Geodesy 87,
43-55 (2013);
2353 * https://doi.org/10.1007/s00190-012-0578-z
2354 * Addenda:
https://geographiclib.sourceforge.io/geod-addenda.html
2355 *
2356 * Copyright (c) Charles
Karney (2011-2022) <karney@alum.mit.edu> and licensed
2357 * under the MIT/X11
License. For more information, see
2358 *
https://geographiclib.sourceforge.io/
2359 */
2360
2361// Load AFTER
geodesic/Math.js and geodesic/Geodesic.js
2362
2363(function(
2364 /**
2365 * @exports geodesic/PolygonArea
2366 * @description Compute the area of geodesic
polygons via the
2367 *
{@link module:geodesic/PolygonArea.PolygonArea PolygonArea}
2368 *
class.
2369 */
2370 p, g, m, a) {
2371 "use strict";
2372
2373 var transit, transitdirect, AreaReduceA,
AreaReduceB;
2374 transit = function(lon1, lon2) {
2375 // Return 1 or -1 if crossing prime
meridian in east or west direction.
2376 // Otherwise return zero. longitude = +/-0 considered to be positive.
2377 // This is (should be?) compatible with
transitdirect which computes
2378 // exactly the parity of
2379 //
int(floor((lon1 + lon12) / 360)) - int(floor(lon1 / 360)))
2380 var lon12 = m.AngDiff(lon1, lon2).d;
2381 lon1 = m.AngNormalize(lon1);
2382 lon2 = m.AngNormalize(lon2);
2383 // edge case lon1 = 180, lon2 = 360->0,
lon12 = 180 to give 1
2384 return lon12 > 0 &&
2385 // lon12 > 0 && lon1 > 0
&& lon2 == 0 implies lon1 == 180
2386 ((lon1 < 0 && lon2 >= 0) ||
(lon1 > 0 && lon2 === 0)) ? 1 :
2387 // non edge case lon1 = -180, lon2 =
-360->-0, lon12 = -180
2388 (lon12 < 0 && lon1 >= 0
&& lon2 < 0 ? -1 : 0);
2389 };
2390
2391 // an alternate version of transit to deal
with longitudes in the direct
2392 // problem.
2393 transitdirect = function(lon1, lon2) {
2394 // We want to compute exactly
2395 //
int(floor(lon2 / 360)) - int(floor(lon1 / 360))
2396 lon1 = lon1 % 720; lon2 = lon2 % 720;
2397 return ((0 <= lon2 && lon2 <
360) || lon2 < -360 ? 0 : 1) -
2398 ((0 <= lon1 && lon1 < 360) || lon1 < -360 ? 0 : 1 );
2399 };
2400
2401 // Reduce Accumulator area
2402 AreaReduceA = function(area, area0,
crossings, reverse, sign) {
2403 area.Remainder(area0);
2404 if (crossings & 1)
2405 area.Add( (area.Sum() < 0 ? 1 : -1) *
area0/2 );
2406 // area is with the clockwise sense. If !reverse convert to
2407 // counter-clockwise convention.
2408 if (!reverse)
2409 area.Negate();
2410 // If sign put area in (-area0/2, area0/2],
else put area in [0, area0)
2411 if (sign) {
2412 if (area.Sum() > area0/2)
2413 area.Add( -area0 );
2414 else if (area.Sum() <= -area0/2)
2415 area.Add( +area0 );
2416 } else {
2417 if (area.Sum() >= area0)
2418 area.Add( -area0 );
2419 else if (area.Sum() < 0)
2420 area.Add( +area0 );
2421 }
2422 return 0 + area.Sum();
2423 };
2424
2425 // Reduce double area
2426 AreaReduceB = function(area, area0,
crossings, reverse, sign) {
2427 area = m.remainder(area, area0);
2428 if (crossings & 1)
2429 area += (area < 0 ? 1 : -1) * area0/2;
2430 // area is with the clockwise sense. If !reverse convert to
2431 // counter-clockwise convention.
2432 if (!reverse)
2433 area *= -1;
2434 // If sign put area in (-area0/2, area0/2],
else put area in [0, area0)
2435 if (sign) {
2436 if (area > area0/2)
2437 area -= area0;
2438 else if (area <= -area0/2)
2439 area += area0;
2440 } else {
2441 if (area >= area0)
2442 area -= area0;
2443 else if (area < 0)
2444 area += area0;
2445 }
2446 return 0 + area;
2447 };
2448
2449 /**
2450 * @class
2451 * @property {number} a the equatorial radius
(meters).
2452 * @property {number} f the flattening.
2453 * @property {bool} polyline whether the
PolygonArea object describes a
2454 *
polyline or a polygon.
2455 * @property {number} num the number of
vertices so far.
2456 * @property {number} lat the current
latitude (degrees).
2457 * @property {number} lon the current
longitude (degrees).
2458 * @summary Initialize a PolygonArea object.
2459 * @classdesc Computes the area and perimeter
of a geodesic polygon.
2460 *
This object is usually instantiated by
2461 *
{@link module:geodesic/Geodesic.Geodesic#Polygon Geodesic.Polygon}.
2462 * @param {object} geod a {@link module:geodesic/Geodesic.Geodesic
2463 *
Geodesic} object.
2464 * @param {bool} [polyline = false] if true
the new PolygonArea object
2465 *
describes a polyline instead of a polygon.
2466 */
2467 p.PolygonArea = function(geod, polyline) {
2468 this._geod = geod;
2469 this.a = this._geod.a;
2470 this.f = this._geod.f;
2471 this._area0 = 4 * Math.PI * geod._c2;
2472 this.polyline = !polyline ? false :
polyline;
2473 this._mask = g.LATITUDE | g.LONGITUDE |
g.DISTANCE |
2474 (this.polyline ? g.NONE : g.AREA |
g.LONG_UNROLL);
2475 if (!this.polyline)
2476 this._areasum = new a.Accumulator(0);
2477 this._perimetersum = new a.Accumulator(0);
2478 this.Clear();
2479 };
2480
2481 /**
2482 * @summary Clear the PolygonArea object,
setting the number of vertices to
2483 * 0.
2484 */
2485 p.PolygonArea.prototype.Clear = function() {
2486 this.num = 0;
2487 this._crossings = 0;
2488 if (!this.polyline)
2489 this._areasum.Set(0);
2490 this._perimetersum.Set(0);
2491 this._lat0 = this._lon0 = this.lat =
this.lon = NaN;
2492 };
2493
2494 /**
2495 * @summary Add the next vertex to the
polygon.
2496 * @param {number} lat the latitude of the
point (degrees).
2497 * @param {number} lon the longitude of the
point (degrees).
2498 * @description This adds an edge from the
current vertex to the new vertex.
2499 */
2500 p.PolygonArea.prototype.AddPoint =
function(lat, lon) {
2501 var t;
2502 if (this.num === 0) {
2503 this._lat0 = this.lat = lat;
2504 this._lon0 = this.lon = lon;
2505 } else {
2506 t = this._geod.Inverse(this.lat,
this.lon, lat, lon, this._mask);
2507 this._perimetersum.Add(t.s12);
2508 if (!this.polyline) {
2509 this._areasum.Add(t.S12);
2510 this._crossings += transit(this.lon,
lon);
2511 }
2512 this.lat = lat;
2513 this.lon = lon;
2514 }
2515 ++this.num;
2516 };
2517
2518 /**
2519 * @summary Add the next edge to the polygon.
2520 * @param {number} azi the azimuth at the
current the point (degrees).
2521 * @param {number} s the length of the edge
(meters).
2522 * @description This specifies the new vertex
in terms of the edge from the
2523 *
current vertex.
2524 */
2525 p.PolygonArea.prototype.AddEdge =
function(azi, s) {
2526 var t;
2527 if (this.num) {
2528 t = this._geod.Direct(this.lat, this.lon,
azi, s, this._mask);
2529 this._perimetersum.Add(s);
2530 if (!this.polyline) {
2531 this._areasum.Add(t.S12);
2532 this._crossings +=
transitdirect(this.lon, t.lon2);
2533 }
2534 this.lat = t.lat2;
2535 this.lon = t.lon2;
2536 }
2537 ++this.num;
2538 };
2539
2540 /**
2541 * @summary Compute the perimeter and area of
the polygon.
2542 * @param {bool} reverse if true then
clockwise (instead of
2543 *
counter-clockwise) traversal counts as a positive area.
2544 * @param {bool} sign if true then return a
signed result for the area if the
2545 *
polygon is traversed in the "wrong" direction instead of
returning the
2546 *
area for the rest of the earth.
2547 * @return {object} r where r.number is the
number of vertices, r.perimeter
2548 * is
the perimeter (meters), and r.area (only returned if polyline is
2549 *
false) is the area (meters<sup>2</sup>).
2550 * @description Arbitrarily complex polygons
are allowed. In the case of
2551 *
self-intersecting polygons the area is accumulated
"algebraically",
2552 *
e.g., the areas of the 2 loops in a figure-8 polygon will partially
2553 *
cancel. If the object is a
polygon (and not a polyline), the perimeter
2554 *
includes the length of a final edge connecting the current point to the
2555 *
initial point. If the object is a
polyline, then area is nan. More
2556 *
points can be added to the polygon after this call.
2557 */
2558 p.PolygonArea.prototype.Compute =
function(reverse, sign) {
2559 var vals = {number: this.num}, t, tempsum;
2560 if (this.num < 2) {
2561 vals.perimeter = 0;
2562 if (!this.polyline)
2563 vals.area = 0;
2564 return vals;
2565 }
2566 if (this.polyline) {
2567 vals.perimeter =
this._perimetersum.Sum();
2568 return vals;
2569 }
2570 t = this._geod.Inverse(this.lat, this.lon,
this._lat0, this._lon0,
2571 this._mask);
2572 vals.perimeter =
this._perimetersum.Sum(t.s12);
2573 tempsum = new a.Accumulator(this._areasum);
2574 tempsum.Add(t.S12);
2575 vals.area = AreaReduceA(tempsum,
this._area0,
2576 this._crossings +
transit(this.lon, this._lon0),
2577 reverse, sign);
2578 return vals;
2579 };
2580
2581 /**
2582 * @summary Compute the perimeter and area of
the polygon with a tentative
2583 *
new vertex.
2584 * @param {number} lat the latitude of the
point (degrees).
2585 * @param {number} lon the longitude of the
point (degrees).
2586 * @param {bool} reverse if true then
clockwise (instead of
2587 *
counter-clockwise) traversal counts as a positive area.
2588 * @param {bool} sign if true then return a
signed result for the area if the
2589 *
polygon is traversed in the "wrong" direction instead of
returning the
2590 *
area for the rest of the earth.
2591 * @return {object} r where r.number is the
number of vertices, r.perimeter
2592 * is
the perimeter (meters), and r.area (only returned if polyline is
2593 *
false) is the area (meters<sup>2</sup>).
2594 * @description A new vertex is *not* added
to the polygon.
2595 */
2596 p.PolygonArea.prototype.TestPoint =
function(lat, lon, reverse, sign) {
2597 var vals = {number: this.num + 1}, t,
tempsum, crossings, i;
2598 if (this.num === 0) {
2599 vals.perimeter = 0;
2600 if (!this.polyline)
2601 vals.area = 0;
2602 return vals;
2603 }
2604 vals.perimeter = this._perimetersum.Sum();
2605 tempsum = this.polyline ? 0 :
this._areasum.Sum();
2606 crossings = this._crossings;
2607 for (i = 0; i < (this.polyline ? 1 : 2);
++i) {
2608 t = this._geod.Inverse(
2609 i === 0 ? this.lat : lat, i === 0 ?
this.lon : lon,
2610 i !== 0 ? this._lat0 : lat, i !== 0 ?
this._lon0 : lon,
2611 this._mask);
2612 vals.perimeter += t.s12;
2613 if (!this.polyline) {
2614 tempsum += t.S12;
2615 crossings += transit(i === 0 ? this.lon
: lon,
2616 i !== 0 ?
this._lon0 : lon);
2617 }
2618 }
2619
2620 if (this.polyline)
2621 return vals;
2622
2623 vals.area = AreaReduceB(tempsum, this._area0,
crossings, reverse, sign);
2624 return vals;
2625 };
2626
2627 /**
2628 * @summary Compute the perimeter and area of
the polygon with a tentative
2629 *
new edge.
2630 * @param {number} azi the azimuth of the
edge (degrees).
2631 * @param {number} s the length of the edge
(meters).
2632 * @param {bool} reverse if true then
clockwise (instead of
2633 *
counter-clockwise) traversal counts as a positive area.
2634 * @param {bool} sign if true then return a
signed result for the area if the
2635 *
polygon is traversed in the "wrong" direction instead of
returning the
2636 *
area for the rest of the earth.
2637 * @return {object} r where r.number is the
number of vertices, r.perimeter
2638 * is
the perimeter (meters), and r.area (only returned if polyline is
2639 *
false) is the area (meters<sup>2</sup>).
2640 * @description A new vertex is *not* added
to the polygon.
2641 */
2642 p.PolygonArea.prototype.TestEdge =
function(azi, s, reverse, sign) {
2643 var vals = {number: this.num ? this.num + 1
: 0}, t, tempsum, crossings;
2644 if (this.num === 0)
2645 return vals;
2646 vals.perimeter = this._perimetersum.Sum() +
s;
2647 if (this.polyline)
2648 return vals;
2649
2650 tempsum = this._areasum.Sum();
2651 crossings = this._crossings;
2652 t = this._geod.Direct(this.lat, this.lon,
azi, s, this._mask);
2653 tempsum += t.S12;
2654 crossings += transitdirect(this.lon,
t.lon2);
2655 crossings += transit(t.lon2, this._lon0);
2656 t = this._geod.Inverse(t.lat2, t.lon2,
this._lat0, this._lon0, this._mask);
2657 vals.perimeter += t.s12;
2658 tempsum += t.S12;
2659
2660 vals.area = AreaReduceB(tempsum,
this._area0, crossings, reverse, sign);
2661 return vals;
2662 };
2663
2664})(geodesic.PolygonArea,
geodesic.Geodesic,
2665 geodesic.Math, geodesic.Accumulator);
2666cb(geodesic);
2667
2668})(function(geo) {
2669 if (typeof module === 'object' &&
module.exports) {
2670 /******** support loading with node's
require ********/
2671 module.exports = geo;
2672 } else if (typeof define === 'function'
&& define.amd) {
2673 /******** support loading with AMD
********/
2674 define('geographiclib-geodesic', [],
function() { return geo; });
2675 } else {
2676 /******** otherwise just pollute our global
namespace ********/
2677 window.geodesic = geo;
2678 }
2679});
2680
تعليقات
إرسال تعليق