EllipsoidGeodesic-924f7301.js 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544
  1. /**
  2. * Cesium - https://github.com/CesiumGS/cesium
  3. *
  4. * Copyright 2011-2020 Cesium Contributors
  5. *
  6. * Licensed under the Apache License, Version 2.0 (the "License");
  7. * you may not use this file except in compliance with the License.
  8. * You may obtain a copy of the License at
  9. *
  10. * http://www.apache.org/licenses/LICENSE-2.0
  11. *
  12. * Unless required by applicable law or agreed to in writing, software
  13. * distributed under the License is distributed on an "AS IS" BASIS,
  14. * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  15. * See the License for the specific language governing permissions and
  16. * limitations under the License.
  17. *
  18. * Columbus View (Pat. Pend.)
  19. *
  20. * Portions licensed separately.
  21. * See https://github.com/CesiumGS/cesium/blob/main/LICENSE.md for full licensing details.
  22. */
  23. define(['exports', './Matrix2-d35cf4b5', './RuntimeError-8952249c', './defaultValue-81eec7ed', './ComponentDatatype-9e86ac8f'], (function (exports, Matrix2, RuntimeError, defaultValue, ComponentDatatype) { 'use strict';
  24. function setConstants(ellipsoidGeodesic) {
  25. const uSquared = ellipsoidGeodesic._uSquared;
  26. const a = ellipsoidGeodesic._ellipsoid.maximumRadius;
  27. const b = ellipsoidGeodesic._ellipsoid.minimumRadius;
  28. const f = (a - b) / a;
  29. const cosineHeading = Math.cos(ellipsoidGeodesic._startHeading);
  30. const sineHeading = Math.sin(ellipsoidGeodesic._startHeading);
  31. const tanU = (1 - f) * Math.tan(ellipsoidGeodesic._start.latitude);
  32. const cosineU = 1.0 / Math.sqrt(1.0 + tanU * tanU);
  33. const sineU = cosineU * tanU;
  34. const sigma = Math.atan2(tanU, cosineHeading);
  35. const sineAlpha = cosineU * sineHeading;
  36. const sineSquaredAlpha = sineAlpha * sineAlpha;
  37. const cosineSquaredAlpha = 1.0 - sineSquaredAlpha;
  38. const cosineAlpha = Math.sqrt(cosineSquaredAlpha);
  39. const u2Over4 = uSquared / 4.0;
  40. const u4Over16 = u2Over4 * u2Over4;
  41. const u6Over64 = u4Over16 * u2Over4;
  42. const u8Over256 = u4Over16 * u4Over16;
  43. const a0 =
  44. 1.0 +
  45. u2Over4 -
  46. (3.0 * u4Over16) / 4.0 +
  47. (5.0 * u6Over64) / 4.0 -
  48. (175.0 * u8Over256) / 64.0;
  49. const a1 = 1.0 - u2Over4 + (15.0 * u4Over16) / 8.0 - (35.0 * u6Over64) / 8.0;
  50. const a2 = 1.0 - 3.0 * u2Over4 + (35.0 * u4Over16) / 4.0;
  51. const a3 = 1.0 - 5.0 * u2Over4;
  52. const distanceRatio =
  53. a0 * sigma -
  54. (a1 * Math.sin(2.0 * sigma) * u2Over4) / 2.0 -
  55. (a2 * Math.sin(4.0 * sigma) * u4Over16) / 16.0 -
  56. (a3 * Math.sin(6.0 * sigma) * u6Over64) / 48.0 -
  57. (Math.sin(8.0 * sigma) * 5.0 * u8Over256) / 512;
  58. const constants = ellipsoidGeodesic._constants;
  59. constants.a = a;
  60. constants.b = b;
  61. constants.f = f;
  62. constants.cosineHeading = cosineHeading;
  63. constants.sineHeading = sineHeading;
  64. constants.tanU = tanU;
  65. constants.cosineU = cosineU;
  66. constants.sineU = sineU;
  67. constants.sigma = sigma;
  68. constants.sineAlpha = sineAlpha;
  69. constants.sineSquaredAlpha = sineSquaredAlpha;
  70. constants.cosineSquaredAlpha = cosineSquaredAlpha;
  71. constants.cosineAlpha = cosineAlpha;
  72. constants.u2Over4 = u2Over4;
  73. constants.u4Over16 = u4Over16;
  74. constants.u6Over64 = u6Over64;
  75. constants.u8Over256 = u8Over256;
  76. constants.a0 = a0;
  77. constants.a1 = a1;
  78. constants.a2 = a2;
  79. constants.a3 = a3;
  80. constants.distanceRatio = distanceRatio;
  81. }
  82. function computeC(f, cosineSquaredAlpha) {
  83. return (
  84. (f * cosineSquaredAlpha * (4.0 + f * (4.0 - 3.0 * cosineSquaredAlpha))) /
  85. 16.0
  86. );
  87. }
  88. function computeDeltaLambda(
  89. f,
  90. sineAlpha,
  91. cosineSquaredAlpha,
  92. sigma,
  93. sineSigma,
  94. cosineSigma,
  95. cosineTwiceSigmaMidpoint
  96. ) {
  97. const C = computeC(f, cosineSquaredAlpha);
  98. return (
  99. (1.0 - C) *
  100. f *
  101. sineAlpha *
  102. (sigma +
  103. C *
  104. sineSigma *
  105. (cosineTwiceSigmaMidpoint +
  106. C *
  107. cosineSigma *
  108. (2.0 * cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint - 1.0)))
  109. );
  110. }
  111. function vincentyInverseFormula(
  112. ellipsoidGeodesic,
  113. major,
  114. minor,
  115. firstLongitude,
  116. firstLatitude,
  117. secondLongitude,
  118. secondLatitude
  119. ) {
  120. const eff = (major - minor) / major;
  121. const l = secondLongitude - firstLongitude;
  122. const u1 = Math.atan((1 - eff) * Math.tan(firstLatitude));
  123. const u2 = Math.atan((1 - eff) * Math.tan(secondLatitude));
  124. const cosineU1 = Math.cos(u1);
  125. const sineU1 = Math.sin(u1);
  126. const cosineU2 = Math.cos(u2);
  127. const sineU2 = Math.sin(u2);
  128. const cc = cosineU1 * cosineU2;
  129. const cs = cosineU1 * sineU2;
  130. const ss = sineU1 * sineU2;
  131. const sc = sineU1 * cosineU2;
  132. let lambda = l;
  133. let lambdaDot = ComponentDatatype.CesiumMath.TWO_PI;
  134. let cosineLambda = Math.cos(lambda);
  135. let sineLambda = Math.sin(lambda);
  136. let sigma;
  137. let cosineSigma;
  138. let sineSigma;
  139. let cosineSquaredAlpha;
  140. let cosineTwiceSigmaMidpoint;
  141. do {
  142. cosineLambda = Math.cos(lambda);
  143. sineLambda = Math.sin(lambda);
  144. const temp = cs - sc * cosineLambda;
  145. sineSigma = Math.sqrt(
  146. cosineU2 * cosineU2 * sineLambda * sineLambda + temp * temp
  147. );
  148. cosineSigma = ss + cc * cosineLambda;
  149. sigma = Math.atan2(sineSigma, cosineSigma);
  150. let sineAlpha;
  151. if (sineSigma === 0.0) {
  152. sineAlpha = 0.0;
  153. cosineSquaredAlpha = 1.0;
  154. } else {
  155. sineAlpha = (cc * sineLambda) / sineSigma;
  156. cosineSquaredAlpha = 1.0 - sineAlpha * sineAlpha;
  157. }
  158. lambdaDot = lambda;
  159. cosineTwiceSigmaMidpoint = cosineSigma - (2.0 * ss) / cosineSquaredAlpha;
  160. if (!isFinite(cosineTwiceSigmaMidpoint)) {
  161. cosineTwiceSigmaMidpoint = 0.0;
  162. }
  163. lambda =
  164. l +
  165. computeDeltaLambda(
  166. eff,
  167. sineAlpha,
  168. cosineSquaredAlpha,
  169. sigma,
  170. sineSigma,
  171. cosineSigma,
  172. cosineTwiceSigmaMidpoint
  173. );
  174. } while (Math.abs(lambda - lambdaDot) > ComponentDatatype.CesiumMath.EPSILON12);
  175. const uSquared =
  176. (cosineSquaredAlpha * (major * major - minor * minor)) / (minor * minor);
  177. const A =
  178. 1.0 +
  179. (uSquared *
  180. (4096.0 + uSquared * (uSquared * (320.0 - 175.0 * uSquared) - 768.0))) /
  181. 16384.0;
  182. const B =
  183. (uSquared *
  184. (256.0 + uSquared * (uSquared * (74.0 - 47.0 * uSquared) - 128.0))) /
  185. 1024.0;
  186. const cosineSquaredTwiceSigmaMidpoint =
  187. cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint;
  188. const deltaSigma =
  189. B *
  190. sineSigma *
  191. (cosineTwiceSigmaMidpoint +
  192. (B *
  193. (cosineSigma * (2.0 * cosineSquaredTwiceSigmaMidpoint - 1.0) -
  194. (B *
  195. cosineTwiceSigmaMidpoint *
  196. (4.0 * sineSigma * sineSigma - 3.0) *
  197. (4.0 * cosineSquaredTwiceSigmaMidpoint - 3.0)) /
  198. 6.0)) /
  199. 4.0);
  200. const distance = minor * A * (sigma - deltaSigma);
  201. const startHeading = Math.atan2(
  202. cosineU2 * sineLambda,
  203. cs - sc * cosineLambda
  204. );
  205. const endHeading = Math.atan2(cosineU1 * sineLambda, cs * cosineLambda - sc);
  206. ellipsoidGeodesic._distance = distance;
  207. ellipsoidGeodesic._startHeading = startHeading;
  208. ellipsoidGeodesic._endHeading = endHeading;
  209. ellipsoidGeodesic._uSquared = uSquared;
  210. }
  211. const scratchCart1 = new Matrix2.Cartesian3();
  212. const scratchCart2 = new Matrix2.Cartesian3();
  213. function computeProperties(ellipsoidGeodesic, start, end, ellipsoid) {
  214. const firstCartesian = Matrix2.Cartesian3.normalize(
  215. ellipsoid.cartographicToCartesian(start, scratchCart2),
  216. scratchCart1
  217. );
  218. const lastCartesian = Matrix2.Cartesian3.normalize(
  219. ellipsoid.cartographicToCartesian(end, scratchCart2),
  220. scratchCart2
  221. );
  222. //>>includeStart('debug', pragmas.debug);
  223. RuntimeError.Check.typeOf.number.greaterThanOrEquals(
  224. "value",
  225. Math.abs(
  226. Math.abs(Matrix2.Cartesian3.angleBetween(firstCartesian, lastCartesian)) - Math.PI
  227. ),
  228. 0.0125
  229. );
  230. //>>includeEnd('debug');
  231. vincentyInverseFormula(
  232. ellipsoidGeodesic,
  233. ellipsoid.maximumRadius,
  234. ellipsoid.minimumRadius,
  235. start.longitude,
  236. start.latitude,
  237. end.longitude,
  238. end.latitude
  239. );
  240. ellipsoidGeodesic._start = Matrix2.Cartographic.clone(
  241. start,
  242. ellipsoidGeodesic._start
  243. );
  244. ellipsoidGeodesic._end = Matrix2.Cartographic.clone(end, ellipsoidGeodesic._end);
  245. ellipsoidGeodesic._start.height = 0;
  246. ellipsoidGeodesic._end.height = 0;
  247. setConstants(ellipsoidGeodesic);
  248. }
  249. /**
  250. * Initializes a geodesic on the ellipsoid connecting the two provided planetodetic points.
  251. *
  252. * @alias EllipsoidGeodesic
  253. * @constructor
  254. *
  255. * @param {Cartographic} [start] The initial planetodetic point on the path.
  256. * @param {Cartographic} [end] The final planetodetic point on the path.
  257. * @param {Ellipsoid} [ellipsoid=Ellipsoid.WGS84] The ellipsoid on which the geodesic lies.
  258. */
  259. function EllipsoidGeodesic(start, end, ellipsoid) {
  260. const e = defaultValue.defaultValue(ellipsoid, Matrix2.Ellipsoid.WGS84);
  261. this._ellipsoid = e;
  262. this._start = new Matrix2.Cartographic();
  263. this._end = new Matrix2.Cartographic();
  264. this._constants = {};
  265. this._startHeading = undefined;
  266. this._endHeading = undefined;
  267. this._distance = undefined;
  268. this._uSquared = undefined;
  269. if (defaultValue.defined(start) && defaultValue.defined(end)) {
  270. computeProperties(this, start, end, e);
  271. }
  272. }
  273. Object.defineProperties(EllipsoidGeodesic.prototype, {
  274. /**
  275. * Gets the ellipsoid.
  276. * @memberof EllipsoidGeodesic.prototype
  277. * @type {Ellipsoid}
  278. * @readonly
  279. */
  280. ellipsoid: {
  281. get: function () {
  282. return this._ellipsoid;
  283. },
  284. },
  285. /**
  286. * Gets the surface distance between the start and end point
  287. * @memberof EllipsoidGeodesic.prototype
  288. * @type {Number}
  289. * @readonly
  290. */
  291. surfaceDistance: {
  292. get: function () {
  293. //>>includeStart('debug', pragmas.debug);
  294. RuntimeError.Check.defined("distance", this._distance);
  295. //>>includeEnd('debug');
  296. return this._distance;
  297. },
  298. },
  299. /**
  300. * Gets the initial planetodetic point on the path.
  301. * @memberof EllipsoidGeodesic.prototype
  302. * @type {Cartographic}
  303. * @readonly
  304. */
  305. start: {
  306. get: function () {
  307. return this._start;
  308. },
  309. },
  310. /**
  311. * Gets the final planetodetic point on the path.
  312. * @memberof EllipsoidGeodesic.prototype
  313. * @type {Cartographic}
  314. * @readonly
  315. */
  316. end: {
  317. get: function () {
  318. return this._end;
  319. },
  320. },
  321. /**
  322. * Gets the heading at the initial point.
  323. * @memberof EllipsoidGeodesic.prototype
  324. * @type {Number}
  325. * @readonly
  326. */
  327. startHeading: {
  328. get: function () {
  329. //>>includeStart('debug', pragmas.debug);
  330. RuntimeError.Check.defined("distance", this._distance);
  331. //>>includeEnd('debug');
  332. return this._startHeading;
  333. },
  334. },
  335. /**
  336. * Gets the heading at the final point.
  337. * @memberof EllipsoidGeodesic.prototype
  338. * @type {Number}
  339. * @readonly
  340. */
  341. endHeading: {
  342. get: function () {
  343. //>>includeStart('debug', pragmas.debug);
  344. RuntimeError.Check.defined("distance", this._distance);
  345. //>>includeEnd('debug');
  346. return this._endHeading;
  347. },
  348. },
  349. });
  350. /**
  351. * Sets the start and end points of the geodesic
  352. *
  353. * @param {Cartographic} start The initial planetodetic point on the path.
  354. * @param {Cartographic} end The final planetodetic point on the path.
  355. */
  356. EllipsoidGeodesic.prototype.setEndPoints = function (start, end) {
  357. //>>includeStart('debug', pragmas.debug);
  358. RuntimeError.Check.defined("start", start);
  359. RuntimeError.Check.defined("end", end);
  360. //>>includeEnd('debug');
  361. computeProperties(this, start, end, this._ellipsoid);
  362. };
  363. /**
  364. * Provides the location of a point at the indicated portion along the geodesic.
  365. *
  366. * @param {Number} fraction The portion of the distance between the initial and final points.
  367. * @param {Cartographic} [result] The object in which to store the result.
  368. * @returns {Cartographic} The location of the point along the geodesic.
  369. */
  370. EllipsoidGeodesic.prototype.interpolateUsingFraction = function (
  371. fraction,
  372. result
  373. ) {
  374. return this.interpolateUsingSurfaceDistance(
  375. this._distance * fraction,
  376. result
  377. );
  378. };
  379. /**
  380. * Provides the location of a point at the indicated distance along the geodesic.
  381. *
  382. * @param {Number} distance The distance from the inital point to the point of interest along the geodesic
  383. * @param {Cartographic} [result] The object in which to store the result.
  384. * @returns {Cartographic} The location of the point along the geodesic.
  385. *
  386. * @exception {DeveloperError} start and end must be set before calling function interpolateUsingSurfaceDistance
  387. */
  388. EllipsoidGeodesic.prototype.interpolateUsingSurfaceDistance = function (
  389. distance,
  390. result
  391. ) {
  392. //>>includeStart('debug', pragmas.debug);
  393. RuntimeError.Check.defined("distance", this._distance);
  394. //>>includeEnd('debug');
  395. const constants = this._constants;
  396. const s = constants.distanceRatio + distance / constants.b;
  397. const cosine2S = Math.cos(2.0 * s);
  398. const cosine4S = Math.cos(4.0 * s);
  399. const cosine6S = Math.cos(6.0 * s);
  400. const sine2S = Math.sin(2.0 * s);
  401. const sine4S = Math.sin(4.0 * s);
  402. const sine6S = Math.sin(6.0 * s);
  403. const sine8S = Math.sin(8.0 * s);
  404. const s2 = s * s;
  405. const s3 = s * s2;
  406. const u8Over256 = constants.u8Over256;
  407. const u2Over4 = constants.u2Over4;
  408. const u6Over64 = constants.u6Over64;
  409. const u4Over16 = constants.u4Over16;
  410. let sigma =
  411. (2.0 * s3 * u8Over256 * cosine2S) / 3.0 +
  412. s *
  413. (1.0 -
  414. u2Over4 +
  415. (7.0 * u4Over16) / 4.0 -
  416. (15.0 * u6Over64) / 4.0 +
  417. (579.0 * u8Over256) / 64.0 -
  418. (u4Over16 - (15.0 * u6Over64) / 4.0 + (187.0 * u8Over256) / 16.0) *
  419. cosine2S -
  420. ((5.0 * u6Over64) / 4.0 - (115.0 * u8Over256) / 16.0) * cosine4S -
  421. (29.0 * u8Over256 * cosine6S) / 16.0) +
  422. (u2Over4 / 2.0 -
  423. u4Over16 +
  424. (71.0 * u6Over64) / 32.0 -
  425. (85.0 * u8Over256) / 16.0) *
  426. sine2S +
  427. ((5.0 * u4Over16) / 16.0 -
  428. (5.0 * u6Over64) / 4.0 +
  429. (383.0 * u8Over256) / 96.0) *
  430. sine4S -
  431. s2 *
  432. ((u6Over64 - (11.0 * u8Over256) / 2.0) * sine2S +
  433. (5.0 * u8Over256 * sine4S) / 2.0) +
  434. ((29.0 * u6Over64) / 96.0 - (29.0 * u8Over256) / 16.0) * sine6S +
  435. (539.0 * u8Over256 * sine8S) / 1536.0;
  436. const theta = Math.asin(Math.sin(sigma) * constants.cosineAlpha);
  437. const latitude = Math.atan((constants.a / constants.b) * Math.tan(theta));
  438. // Redefine in terms of relative argument of latitude.
  439. sigma = sigma - constants.sigma;
  440. const cosineTwiceSigmaMidpoint = Math.cos(2.0 * constants.sigma + sigma);
  441. const sineSigma = Math.sin(sigma);
  442. const cosineSigma = Math.cos(sigma);
  443. const cc = constants.cosineU * cosineSigma;
  444. const ss = constants.sineU * sineSigma;
  445. const lambda = Math.atan2(
  446. sineSigma * constants.sineHeading,
  447. cc - ss * constants.cosineHeading
  448. );
  449. const l =
  450. lambda -
  451. computeDeltaLambda(
  452. constants.f,
  453. constants.sineAlpha,
  454. constants.cosineSquaredAlpha,
  455. sigma,
  456. sineSigma,
  457. cosineSigma,
  458. cosineTwiceSigmaMidpoint
  459. );
  460. if (defaultValue.defined(result)) {
  461. result.longitude = this._start.longitude + l;
  462. result.latitude = latitude;
  463. result.height = 0.0;
  464. return result;
  465. }
  466. return new Matrix2.Cartographic(this._start.longitude + l, latitude, 0.0);
  467. };
  468. exports.EllipsoidGeodesic = EllipsoidGeodesic;
  469. }));
  470. //# sourceMappingURL=EllipsoidGeodesic-924f7301.js.map