diff --git a/scripts/format-scope.json b/scripts/format-scope.json index 6a1b853dc..3278f4678 100644 --- a/scripts/format-scope.json +++ b/scripts/format-scope.json @@ -1148,5 +1148,6 @@ "src/celestialRing/view.js", "src/layers/satellites/render/navara.js", "src/layers/satellites/render/navara.test.mjs", - "src/layers/satellites/render/cesium.test.mjs" + "src/layers/satellites/render/cesium.test.mjs", + "src/geo/ellipsoidGeodesic.test.mjs" ] diff --git a/scripts/package-boundaries.json b/scripts/package-boundaries.json index 7aa10025d..366b5cbd1 100644 --- a/scripts/package-boundaries.json +++ b/scripts/package-boundaries.json @@ -528,6 +528,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -693,6 +694,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -816,6 +818,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -884,6 +887,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -943,6 +947,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -988,6 +993,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1032,6 +1038,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1094,6 +1101,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1136,6 +1144,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1179,6 +1188,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1231,6 +1241,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1274,6 +1285,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1331,6 +1343,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1370,6 +1383,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1427,6 +1441,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1487,6 +1502,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1523,6 +1539,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1572,6 +1589,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -1605,6 +1623,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -2560,6 +2579,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -3171,6 +3191,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -3270,6 +3291,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -3753,6 +3775,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", @@ -3807,6 +3830,7 @@ "src/geo/math.js", "src/geo/ellipsoid.js", "src/geo/cartographic.js", + "src/geo/ellipsoidGeodesic.js", "src/geo/cartesian2.js", "src/geo/cartesian3.js", "src/geo/matrix3.js", diff --git a/src/geo/cartesian3.js b/src/geo/cartesian3.js index 18f8ef3ad..aea1a5ba9 100644 --- a/src/geo/cartesian3.js +++ b/src/geo/cartesian3.js @@ -387,3 +387,19 @@ export function fromDegreesArrayHeights(coordinates, ellipsoid, result) { } return result; } + +const angleBetweenScratch = create(); +const angleBetweenScratch2 = create(); +export function angleBetween(left, right) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + } + normalize(left, angleBetweenScratch); + normalize(right, angleBetweenScratch2); + const cosine = dot(angleBetweenScratch, angleBetweenScratch2); + const sine = magnitude( + cross(angleBetweenScratch, angleBetweenScratch2, angleBetweenScratch), + ); + return Math.atan2(sine, cosine); +} diff --git a/src/geo/ellipsoidGeodesic.js b/src/geo/ellipsoidGeodesic.js new file mode 100644 index 000000000..ac87a6102 --- /dev/null +++ b/src/geo/ellipsoidGeodesic.js @@ -0,0 +1,527 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — EllipsoidGeodesic.js. +// Vincenty inverse and surface interpolation, preserving operation order. +import * as Cartesian3 from './cartesian3.js'; +import * as Cartographic from './cartographic.js'; +import * as Ellipsoid from './ellipsoid.js'; +import * as CesiumMath from './math.js'; +import { Check } from './check.js'; +import { DEBUG } from './debug.js'; +const defined = (value) => value !== null && value !== undefined; +function setConstants(ellipsoidGeodesic) { + const uSquared = ellipsoidGeodesic._uSquared; + const a = ellipsoidGeodesic._ellipsoid.maximumRadius; + const b = ellipsoidGeodesic._ellipsoid.minimumRadius; + const f = (a - b) / a; + + const cosineHeading = Math.cos(ellipsoidGeodesic._startHeading); + const sineHeading = Math.sin(ellipsoidGeodesic._startHeading); + + const tanU = (1 - f) * Math.tan(ellipsoidGeodesic._start.latitude); + + const cosineU = 1.0 / Math.sqrt(1.0 + tanU * tanU); + const sineU = cosineU * tanU; + + const sigma = Math.atan2(tanU, cosineHeading); + + const sineAlpha = cosineU * sineHeading; + const sineSquaredAlpha = sineAlpha * sineAlpha; + + const cosineSquaredAlpha = 1.0 - sineSquaredAlpha; + const cosineAlpha = Math.sqrt(cosineSquaredAlpha); + + const u2Over4 = uSquared / 4.0; + const u4Over16 = u2Over4 * u2Over4; + const u6Over64 = u4Over16 * u2Over4; + const u8Over256 = u4Over16 * u4Over16; + + const a0 = + 1.0 + + u2Over4 - + (3.0 * u4Over16) / 4.0 + + (5.0 * u6Over64) / 4.0 - + (175.0 * u8Over256) / 64.0; + const a1 = 1.0 - u2Over4 + (15.0 * u4Over16) / 8.0 - (35.0 * u6Over64) / 8.0; + const a2 = 1.0 - 3.0 * u2Over4 + (35.0 * u4Over16) / 4.0; + const a3 = 1.0 - 5.0 * u2Over4; + + const distanceRatio = + a0 * sigma - + (a1 * Math.sin(2.0 * sigma) * u2Over4) / 2.0 - + (a2 * Math.sin(4.0 * sigma) * u4Over16) / 16.0 - + (a3 * Math.sin(6.0 * sigma) * u6Over64) / 48.0 - + (Math.sin(8.0 * sigma) * 5.0 * u8Over256) / 512; + + const constants = ellipsoidGeodesic._constants; + + constants.a = a; + constants.b = b; + constants.f = f; + constants.cosineHeading = cosineHeading; + constants.sineHeading = sineHeading; + constants.tanU = tanU; + constants.cosineU = cosineU; + constants.sineU = sineU; + constants.sigma = sigma; + constants.sineAlpha = sineAlpha; + constants.sineSquaredAlpha = sineSquaredAlpha; + constants.cosineSquaredAlpha = cosineSquaredAlpha; + constants.cosineAlpha = cosineAlpha; + constants.u2Over4 = u2Over4; + constants.u4Over16 = u4Over16; + constants.u6Over64 = u6Over64; + constants.u8Over256 = u8Over256; + constants.a0 = a0; + constants.a1 = a1; + constants.a2 = a2; + constants.a3 = a3; + constants.distanceRatio = distanceRatio; +} + +function computeC(f, cosineSquaredAlpha) { + return ( + (f * cosineSquaredAlpha * (4.0 + f * (4.0 - 3.0 * cosineSquaredAlpha))) / + 16.0 + ); +} + +function computeDeltaLambda( + f, + sineAlpha, + cosineSquaredAlpha, + sigma, + sineSigma, + cosineSigma, + cosineTwiceSigmaMidpoint, +) { + const C = computeC(f, cosineSquaredAlpha); + + return ( + (1.0 - C) * + f * + sineAlpha * + (sigma + + C * + sineSigma * + (cosineTwiceSigmaMidpoint + + C * + cosineSigma * + (2.0 * cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint - 1.0))) + ); +} + +function vincentyInverseFormula( + ellipsoidGeodesic, + major, + minor, + firstLongitude, + firstLatitude, + secondLongitude, + secondLatitude, +) { + const eff = (major - minor) / major; + const l = secondLongitude - firstLongitude; + + const u1 = Math.atan((1 - eff) * Math.tan(firstLatitude)); + const u2 = Math.atan((1 - eff) * Math.tan(secondLatitude)); + + const cosineU1 = Math.cos(u1); + const sineU1 = Math.sin(u1); + const cosineU2 = Math.cos(u2); + const sineU2 = Math.sin(u2); + + const cc = cosineU1 * cosineU2; + const cs = cosineU1 * sineU2; + const ss = sineU1 * sineU2; + const sc = sineU1 * cosineU2; + + let lambda = l; + let lambdaDot = CesiumMath.TWO_PI; + + let cosineLambda = Math.cos(lambda); + let sineLambda = Math.sin(lambda); + + let sigma; + let cosineSigma; + let sineSigma; + let cosineSquaredAlpha; + let cosineTwiceSigmaMidpoint; + + do { + cosineLambda = Math.cos(lambda); + sineLambda = Math.sin(lambda); + + const temp = cs - sc * cosineLambda; + sineSigma = Math.sqrt( + cosineU2 * cosineU2 * sineLambda * sineLambda + temp * temp, + ); + cosineSigma = ss + cc * cosineLambda; + + sigma = Math.atan2(sineSigma, cosineSigma); + + let sineAlpha; + + if (sineSigma === 0.0) { + sineAlpha = 0.0; + cosineSquaredAlpha = 1.0; + } else { + sineAlpha = (cc * sineLambda) / sineSigma; + cosineSquaredAlpha = 1.0 - sineAlpha * sineAlpha; + } + + lambdaDot = lambda; + + cosineTwiceSigmaMidpoint = cosineSigma - (2.0 * ss) / cosineSquaredAlpha; + + if (!isFinite(cosineTwiceSigmaMidpoint)) { + cosineTwiceSigmaMidpoint = 0.0; + } + + lambda = + l + + computeDeltaLambda( + eff, + sineAlpha, + cosineSquaredAlpha, + sigma, + sineSigma, + cosineSigma, + cosineTwiceSigmaMidpoint, + ); + } while (Math.abs(lambda - lambdaDot) > CesiumMath.EPSILON12); + + const uSquared = + (cosineSquaredAlpha * (major * major - minor * minor)) / (minor * minor); + const A = + 1.0 + + (uSquared * + (4096.0 + uSquared * (uSquared * (320.0 - 175.0 * uSquared) - 768.0))) / + 16384.0; + const B = + (uSquared * + (256.0 + uSquared * (uSquared * (74.0 - 47.0 * uSquared) - 128.0))) / + 1024.0; + + const cosineSquaredTwiceSigmaMidpoint = + cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint; + const deltaSigma = + B * + sineSigma * + (cosineTwiceSigmaMidpoint + + (B * + (cosineSigma * (2.0 * cosineSquaredTwiceSigmaMidpoint - 1.0) - + (B * + cosineTwiceSigmaMidpoint * + (4.0 * sineSigma * sineSigma - 3.0) * + (4.0 * cosineSquaredTwiceSigmaMidpoint - 3.0)) / + 6.0)) / + 4.0); + + const distance = minor * A * (sigma - deltaSigma); + + const startHeading = Math.atan2( + cosineU2 * sineLambda, + cs - sc * cosineLambda, + ); + const endHeading = Math.atan2(cosineU1 * sineLambda, cs * cosineLambda - sc); + + ellipsoidGeodesic._distance = distance; + ellipsoidGeodesic._startHeading = startHeading; + ellipsoidGeodesic._endHeading = endHeading; + ellipsoidGeodesic._uSquared = uSquared; +} + +const scratchCart1 = Cartesian3.create(); +const scratchCart2 = Cartesian3.create(); +function computeProperties(ellipsoidGeodesic, start, end, ellipsoid) { + const firstCartesian = Cartesian3.normalize( + ellipsoid.cartographicToCartesian(start, scratchCart2), + scratchCart1, + ); + const lastCartesian = Cartesian3.normalize( + ellipsoid.cartographicToCartesian(end, scratchCart2), + scratchCart2, + ); + + if (DEBUG) { + Check.typeOf.number.greaterThanOrEquals( + 'value', + Math.abs( + Math.abs(Cartesian3.angleBetween(firstCartesian, lastCartesian)) - + Math.PI, + ), + 0.0125, + ); + } + + vincentyInverseFormula( + ellipsoidGeodesic, + ellipsoid.maximumRadius, + ellipsoid.minimumRadius, + start.longitude, + start.latitude, + end.longitude, + end.latitude, + ); + + ellipsoidGeodesic._start = Cartographic.clone( + start, + ellipsoidGeodesic._start, + ); + ellipsoidGeodesic._end = Cartographic.clone(end, ellipsoidGeodesic._end); + ellipsoidGeodesic._start.height = 0; + ellipsoidGeodesic._end.height = 0; + + setConstants(ellipsoidGeodesic); +} + +/** + * Initializes a geodesic on the ellipsoid connecting the two provided planetodetic points. + * + * @alias EllipsoidGeodesic + * @constructor + * + * @param {Cartographic} [start] The initial planetodetic point on the path. + * @param {Cartographic} [end] The final planetodetic point on the path. + * @param {Ellipsoid} [ellipsoid=Ellipsoid.WGS84] The ellipsoid on which the geodesic lies. + */ +function EllipsoidGeodesic(start, end, ellipsoid) { + const e = ellipsoid ?? Ellipsoid.WGS84; + this._ellipsoid = e; + this._start = Cartographic.create(); + this._end = Cartographic.create(); + + this._constants = {}; + this._startHeading = undefined; + this._endHeading = undefined; + this._distance = undefined; + this._uSquared = undefined; + + if (defined(start) && defined(end)) { + computeProperties(this, start, end, e); + } +} + +Object.defineProperties(EllipsoidGeodesic.prototype, { + /** + * Gets the ellipsoid. + * @memberof EllipsoidGeodesic.prototype + * @type {Ellipsoid} + * @readonly + */ + ellipsoid: { + get: function () { + return this._ellipsoid; + }, + }, + + /** + * Gets the surface distance between the start and end point + * @memberof EllipsoidGeodesic.prototype + * @type {number} + * @readonly + */ + surfaceDistance: { + get: function () { + if (DEBUG) { + Check.defined('distance', this._distance); + } + + return this._distance; + }, + }, + + /** + * Gets the initial planetodetic point on the path. + * @memberof EllipsoidGeodesic.prototype + * @type {Cartographic} + * @readonly + */ + start: { + get: function () { + return this._start; + }, + }, + + /** + * Gets the final planetodetic point on the path. + * @memberof EllipsoidGeodesic.prototype + * @type {Cartographic} + * @readonly + */ + end: { + get: function () { + return this._end; + }, + }, + + /** + * Gets the heading at the initial point. + * @memberof EllipsoidGeodesic.prototype + * @type {number} + * @readonly + */ + startHeading: { + get: function () { + if (DEBUG) { + Check.defined('distance', this._distance); + } + + return this._startHeading; + }, + }, + + /** + * Gets the heading at the final point. + * @memberof EllipsoidGeodesic.prototype + * @type {number} + * @readonly + */ + endHeading: { + get: function () { + if (DEBUG) { + Check.defined('distance', this._distance); + } + + return this._endHeading; + }, + }, +}); + +/** + * Sets the start and end points of the geodesic + * + * @param {Cartographic} start The initial planetodetic point on the path. + * @param {Cartographic} end The final planetodetic point on the path. + */ +EllipsoidGeodesic.prototype.setEndPoints = function (start, end) { + if (DEBUG) { + Check.defined('start', start); + Check.defined('end', end); + } + + computeProperties(this, start, end, this._ellipsoid); +}; + +/** + * Provides the location of a point at the indicated portion along the geodesic. + * + * @param {number} fraction The portion of the distance between the initial and final points. + * @param {Cartographic} [result] The object in which to store the result. + * @returns {Cartographic} The location of the point along the geodesic. + */ +EllipsoidGeodesic.prototype.interpolateUsingFraction = function ( + fraction, + result, +) { + return this.interpolateUsingSurfaceDistance( + this._distance * fraction, + result, + ); +}; + +/** + * Provides the location of a point at the indicated distance along the geodesic. + * + * @param {number} distance The distance from the initial point to the point of interest along the geodesic + * @param {Cartographic} [result] The object in which to store the result. + * @returns {Cartographic} The location of the point along the geodesic. + * + * @exception {DeveloperError} start and end must be set before calling function interpolateUsingSurfaceDistance + */ +EllipsoidGeodesic.prototype.interpolateUsingSurfaceDistance = function ( + distance, + result, +) { + if (DEBUG) { + Check.defined('distance', this._distance); + } + + const constants = this._constants; + + const s = constants.distanceRatio + distance / constants.b; + + const cosine2S = Math.cos(2.0 * s); + const cosine4S = Math.cos(4.0 * s); + const cosine6S = Math.cos(6.0 * s); + const sine2S = Math.sin(2.0 * s); + const sine4S = Math.sin(4.0 * s); + const sine6S = Math.sin(6.0 * s); + const sine8S = Math.sin(8.0 * s); + + const s2 = s * s; + const s3 = s * s2; + + const u8Over256 = constants.u8Over256; + const u2Over4 = constants.u2Over4; + const u6Over64 = constants.u6Over64; + const u4Over16 = constants.u4Over16; + let sigma = + (2.0 * s3 * u8Over256 * cosine2S) / 3.0 + + s * + (1.0 - + u2Over4 + + (7.0 * u4Over16) / 4.0 - + (15.0 * u6Over64) / 4.0 + + (579.0 * u8Over256) / 64.0 - + (u4Over16 - (15.0 * u6Over64) / 4.0 + (187.0 * u8Over256) / 16.0) * + cosine2S - + ((5.0 * u6Over64) / 4.0 - (115.0 * u8Over256) / 16.0) * cosine4S - + (29.0 * u8Over256 * cosine6S) / 16.0) + + (u2Over4 / 2.0 - + u4Over16 + + (71.0 * u6Over64) / 32.0 - + (85.0 * u8Over256) / 16.0) * + sine2S + + ((5.0 * u4Over16) / 16.0 - + (5.0 * u6Over64) / 4.0 + + (383.0 * u8Over256) / 96.0) * + sine4S - + s2 * + ((u6Over64 - (11.0 * u8Over256) / 2.0) * sine2S + + (5.0 * u8Over256 * sine4S) / 2.0) + + ((29.0 * u6Over64) / 96.0 - (29.0 * u8Over256) / 16.0) * sine6S + + (539.0 * u8Over256 * sine8S) / 1536.0; + + const theta = Math.asin(Math.sin(sigma) * constants.cosineAlpha); + const latitude = Math.atan((constants.a / constants.b) * Math.tan(theta)); + + // Redefine in terms of relative argument of latitude. + sigma = sigma - constants.sigma; + + const cosineTwiceSigmaMidpoint = Math.cos(2.0 * constants.sigma + sigma); + + const sineSigma = Math.sin(sigma); + const cosineSigma = Math.cos(sigma); + + const cc = constants.cosineU * cosineSigma; + const ss = constants.sineU * sineSigma; + + const lambda = Math.atan2( + sineSigma * constants.sineHeading, + cc - ss * constants.cosineHeading, + ); + + const l = + lambda - + computeDeltaLambda( + constants.f, + constants.sineAlpha, + constants.cosineSquaredAlpha, + sigma, + sineSigma, + cosineSigma, + cosineTwiceSigmaMidpoint, + ); + + if (defined(result)) { + result.longitude = this._start.longitude + l; + result.latitude = latitude; + result.height = 0.0; + return result; + } + + return Cartographic.create(this._start.longitude + l, latitude, 0.0); +}; +export function create(start, end, ellipsoid) { + return new EllipsoidGeodesic(start, end, ellipsoid); +} diff --git a/src/geo/ellipsoidGeodesic.test.mjs b/src/geo/ellipsoidGeodesic.test.mjs new file mode 100644 index 000000000..0ab080da1 --- /dev/null +++ b/src/geo/ellipsoidGeodesic.test.mjs @@ -0,0 +1,82 @@ +import test from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as geo from './index.js'; +const tuple = (p) => [p.longitude, p.latitude, p.height]; +test('Vincenty distance, headings and interpolation are exactly the Cesium oracle', () => { + let checked = 0; + for (let i = 0; i < 120; i++) { + const a = geo.cartographic.fromDegrees( + ((i * 43.17) % 340) - 170, + ((i * 17.23) % 150) - 75, + 123, + ); + const b = geo.cartographic.fromDegrees( + ((i * 83.31 + 20) % 340) - 170, + ((i * 7.43 + 9) % 150) - 75, + 567, + ); + let c; + try { + c = new Cesium.EllipsoidGeodesic(a, b); + } catch { + assert.throws(() => geo.ellipsoidGeodesic.create(a, b)); + continue; + } + const g = geo.ellipsoidGeodesic.create(a, b); + checked++; + for (const key of ['surfaceDistance', 'startHeading', 'endHeading']) + assert.equal(g[key], c[key]); + assert.equal(g.start.height, 0); + assert.equal(g.end.height, 0); + for (const f of [-0.2, 0, 0.25, 0.5, 0.75, 1, 1.2]) + assert.deepEqual( + tuple(g.interpolateUsingFraction(f)), + tuple(c.interpolateUsingFraction(f)), + ); + } + assert.ok(checked > 100); +}); +test('endpoint reset, result reuse, custom sphere and dateline/polar paths retain native results', () => { + const g = geo.ellipsoidGeodesic.create( + undefined, + undefined, + geo.ellipsoid.create(1e7, 1e7, 1e7), + ); + const c = new Cesium.EllipsoidGeodesic( + undefined, + undefined, + new Cesium.Ellipsoid(1e7, 1e7, 1e7), + ); + for (const [lon1, lat1, lon2, lat2] of [ + [179.8, 70, -179.8, 69.9], + [0, 89.99, 120, 50], + [12, -50, -80, 70], + ]) { + const a = geo.cartographic.fromDegrees(lon1, lat1), + b = geo.cartographic.fromDegrees(lon2, lat2); + g.setEndPoints(a, b); + c.setEndPoints(a, b); + assert.equal(g.surfaceDistance, c.surfaceDistance); + const result = geo.cartographic.create(); + assert.equal(g.interpolateUsingSurfaceDistance(12345, result), result); + assert.deepEqual( + tuple(result), + tuple(c.interpolateUsingSurfaceDistance(12345)), + ); + } +}); +test('unset and near-antipodal geodesics reject the same debug inputs', () => { + assert.throws(() => geo.ellipsoidGeodesic.create().surfaceDistance, { + name: 'DeveloperError', + }); + const a = geo.cartographic.fromDegrees(0, 0), + b = geo.cartographic.fromDegrees(180, 0); + assert.throws(() => geo.ellipsoidGeodesic.create(a, b), { + name: 'DeveloperError', + }); + assert.throws( + () => geo.ellipsoidGeodesic.create().setEndPoints(undefined, b), + { name: 'DeveloperError' }, + ); +}); diff --git a/src/geo/index.js b/src/geo/index.js index c86b10401..75309996f 100644 --- a/src/geo/index.js +++ b/src/geo/index.js @@ -13,3 +13,5 @@ export * as matrix3 from './matrix3.js'; export * as matrix4 from './matrix4.js'; export * as quaternion from './quaternion.js'; export * as transforms from './transforms.js'; + +export * as ellipsoidGeodesic from './ellipsoidGeodesic.js'; diff --git a/src/layers/launches/paths.js b/src/layers/launches/paths.js index 72819ffb8..e8df3dd0e 100644 --- a/src/layers/launches/paths.js +++ b/src/layers/launches/paths.js @@ -1,4 +1,4 @@ -import * as Cesium from 'cesium'; +import * as geo from '../../geo/index.js'; import { EARTH_ROTATION_RAD_PER_SEC, PROJECTED_ASCENT_ROTATION_SEC, @@ -31,7 +31,7 @@ export function createPaths({ state: layerState, services, parts, source }) { if (!orbitPath?.length) return 5400; const meanRadius = orbitPath.reduce( - (total, point) => total + Cesium.Cartesian3.magnitude(point), + (total, point) => total + geo.cartesian3.magnitude(point), 0, ) / orbitPath.length; return 2 * Math.PI * Math.sqrt(meanRadius ** 3 / 3.986004418e14); @@ -51,20 +51,20 @@ export function createPaths({ state: layerState, services, parts, source }) { : orbitName.includes('medium') ? 20200000 : 550000; - const radius = Cesium.Ellipsoid.WGS84.maximumRadius + altitude; - const longitude = Cesium.Math.toRadians(launch.lon); - const latitude = Cesium.Math.toRadians(launch.lat); - const up = new Cesium.Cartesian3( + const radius = geo.ellipsoid.WGS84.maximumRadius + altitude; + const longitude = geo.math.toRadians(launch.lon); + const latitude = geo.math.toRadians(launch.lat); + const up = geo.cartesian3.create( Math.cos(latitude) * Math.cos(longitude), Math.cos(latitude) * Math.sin(longitude), Math.sin(latitude), ); - const east = new Cesium.Cartesian3( + const east = geo.cartesian3.create( -Math.sin(longitude), Math.cos(longitude), 0, ); - const north = new Cesium.Cartesian3( + const north = geo.cartesian3.create( -Math.sin(latitude) * Math.cos(longitude), -Math.sin(latitude) * Math.sin(longitude), Math.cos(latitude), @@ -82,60 +82,56 @@ export function createPaths({ state: layerState, services, parts, source }) { : isWesternNorthAmerica ? 190 : 90; - const launchAzimuth = Cesium.Math.toRadians(launchAzimuthDeg); - const forward = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.add( - Cesium.Cartesian3.multiplyByScalar( + const launchAzimuth = geo.math.toRadians(launchAzimuthDeg); + const forward = geo.cartesian3.normalize( + geo.cartesian3.add( + geo.cartesian3.multiplyByScalar( north, Math.cos(launchAzimuth), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( east, Math.sin(launchAzimuth), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); // A reconstructed orbit is an estimate, not a claim that the vehicle was // inserted directly above the pad. Offset the orbital plane downrange by a // small launch-to-insertion arc so a top-down view does not draw the entire // orbit on top of the launch site. The ascent still joins this ring at its // propagated insertion reference. - const insertionArc = Cesium.Math.toRadians(isPolar ? 8 : 12); - const orbitAnchor = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.add( - Cesium.Cartesian3.multiplyByScalar( + const insertionArc = geo.math.toRadians(isPolar ? 8 : 12); + const orbitAnchor = geo.cartesian3.normalize( + geo.cartesian3.add( + geo.cartesian3.multiplyByScalar( up, Math.cos(insertionArc), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( forward, Math.sin(insertionArc), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); - const planeNormal = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.cross(orbitAnchor, forward, new Cesium.Cartesian3()), - new Cesium.Cartesian3(), + const planeNormal = geo.cartesian3.normalize( + geo.cartesian3.cross(orbitAnchor, forward, geo.cartesian3.create()), + geo.cartesian3.create(), ); - const crossTrack = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.cross( - planeNormal, - orbitAnchor, - new Cesium.Cartesian3(), - ), - new Cesium.Cartesian3(), + const crossTrack = geo.cartesian3.normalize( + geo.cartesian3.cross(planeNormal, orbitAnchor, geo.cartesian3.create()), + geo.cartesian3.create(), ); return Array.from({ length: 97 }, (_, index) => { const angle = (index / 96) * Math.PI * 2; - return new Cesium.Cartesian3( + return geo.cartesian3.create( radius * (Math.cos(angle) * orbitAnchor.x + Math.sin(angle) * crossTrack.x), radius * @@ -151,21 +147,21 @@ export function createPaths({ state: layerState, services, parts, source }) { // Degenerate returns clone into `result` when provided — handing back a // path vertex would let an in-place caller mutate the path geometry. if (path.length === 1) - return result ? Cesium.Cartesian3.clone(path[0], result) : path[0]; + return result ? geo.cartesian3.clone(path[0], result) : path[0]; let distances = layerState._pathDistanceCache.get(path); if (!distances) { distances = new Float64Array(path.length); for (let index = 1; index < path.length; index++) { distances[index] = distances[index - 1] + - Cesium.Cartesian3.distance(path[index - 1], path[index]); + geo.cartesian3.distance(path[index - 1], path[index]); } layerState._pathDistanceCache.set(path, distances); } const totalDistance = distances.at(-1); if (!(totalDistance > 0)) - return result ? Cesium.Cartesian3.clone(path[0], result) : path[0]; - const targetDistance = Cesium.Math.clamp(progress, 0, 1) * totalDistance; + return result ? geo.cartesian3.clone(path[0], result) : path[0]; + const targetDistance = geo.math.clamp(progress, 0, 1) * totalDistance; let low = 1; let high = distances.length - 1; while (low < high) { @@ -179,11 +175,11 @@ export function createPaths({ state: layerState, services, parts, source }) { segmentDistance > 0 ? (targetDistance - distances[index]) / segmentDistance : 0; - return Cesium.Cartesian3.lerp( + return geo.cartesian3.lerp( path[index], path[index + 1], segmentProgress, - result || new Cesium.Cartesian3(), + result || geo.cartesian3.create(), ); } @@ -198,7 +194,7 @@ export function createPaths({ state: layerState, services, parts, source }) { function orbitProgressAtTime(nowMs, periodSec) { const period = Math.max(1, Number(periodSec) || 1); - return Cesium.Math.mod((Number(nowMs) || 0) / 1000, period) / period; + return geo.math.mod((Number(nowMs) || 0) / 1000, period) / period; } /** @@ -216,13 +212,13 @@ export function createPaths({ state: layerState, services, parts, source }) { insertionPosition, samples = 512, ) { - const ellipsoid = Cesium.Ellipsoid.WGS84; + const ellipsoid = geo.ellipsoid.WGS84; const origin = ellipsoid.cartesianToCartographic(launchPosition); const insertion = ellipsoid.cartesianToCartographic(insertionPosition); if (!origin || !insertion) return [launchPosition, insertionPosition]; - const geodesic = new Cesium.EllipsoidGeodesic( - new Cesium.Cartographic(origin.longitude, origin.latitude), - new Cesium.Cartographic(insertion.longitude, insertion.latitude), + const geodesic = geo.ellipsoidGeodesic.create( + geo.cartographic.create(origin.longitude, origin.latitude), + geo.cartographic.create(insertion.longitude, insertion.latitude), ellipsoid, ); return Array.from({ length: samples + 1 }, (_, index) => { @@ -232,7 +228,7 @@ export function createPaths({ state: layerState, services, parts, source }) { const horizontalProgress = progress ** 4; const cartographic = geodesic.interpolateUsingFraction( horizontalProgress, - new Cesium.Cartographic(), + geo.cartographic.create(), ); // Approximate the inertial eastward lead accumulated during a ten-minute // ascent. The p³ envelope keeps liftoff nearly vertical, peaks during the @@ -241,13 +237,13 @@ export function createPaths({ state: layerState, services, parts, source }) { EARTH_ROTATION_RAD_PER_SEC * PROJECTED_ASCENT_ROTATION_SEC * Math.sin(Math.PI * progress ** 3); - cartographic.longitude = Cesium.Math.negativePiToPi( + cartographic.longitude = geo.math.negativePiToPi( cartographic.longitude + rotationalLead, ); - cartographic.height = Cesium.Math.lerp( + cartographic.height = geo.math.lerp( Math.max(0, origin.height), Math.max(0, insertion.height), - Math.sin(progress * Cesium.Math.PI_OVER_TWO), + Math.sin(progress * geo.math.PI_OVER_TWO), ); return ellipsoid.cartographicToCartesian(cartographic); }); @@ -255,18 +251,18 @@ export function createPaths({ state: layerState, services, parts, source }) { function nearestOrbitIndex(orbitPath, referencePosition) { if (!orbitPath?.length || !referencePosition) return 0; - const referenceDirection = Cesium.Cartesian3.normalize( + const referenceDirection = geo.cartesian3.normalize( referencePosition, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); let bestIndex = 0; let bestDot = -Number.MAX_VALUE; orbitPath.forEach((candidate, index) => { - const direction = Cesium.Cartesian3.normalize( + const direction = geo.cartesian3.normalize( candidate, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); - const dot = Cesium.Cartesian3.dot(referenceDirection, direction); + const dot = geo.cartesian3.dot(referenceDirection, direction); if (dot > bestDot) { bestDot = dot; bestIndex = index; @@ -280,29 +276,29 @@ export function createPaths({ state: layerState, services, parts, source }) { const first = orbitPath[0]; const last = orbitPath.at(-1); const isClosed = - orbitPath.length > 2 && Cesium.Cartesian3.distance(first, last) < 1000; + orbitPath.length > 2 && geo.cartesian3.distance(first, last) < 1000; const core = isClosed ? orbitPath.slice(0, -1) : orbitPath.slice(); if (!core.length) return orbitPath.slice(); - const index = Cesium.Math.mod(insertionIndex, core.length); + const index = geo.math.mod(insertionIndex, core.length); const rotated = [...core.slice(index), ...core.slice(0, index)]; rotated.push(rotated[0]); return rotated; } function surfaceSafeSegment(startPosition, endPosition) { - const ellipsoid = Cesium.Ellipsoid.WGS84; + const ellipsoid = geo.ellipsoid.WGS84; const start = ellipsoid.cartesianToCartographic(startPosition); const end = ellipsoid.cartesianToCartographic(endPosition); if (!start || !end) return [startPosition, endPosition]; - const geodesic = new Cesium.EllipsoidGeodesic( - new Cesium.Cartographic(start.longitude, start.latitude), - new Cesium.Cartographic(end.longitude, end.latitude), + const geodesic = geo.ellipsoidGeodesic.create( + geo.cartographic.create(start.longitude, start.latitude), + geo.cartographic.create(end.longitude, end.latitude), ellipsoid, ); // Replay advances in real time across this path. Keep the geometry // surface-safe, but give the animated marker and chase camera enough // samples that they do not visibly pause at long segment boundaries. - const steps = Cesium.Math.clamp( + const steps = geo.math.clamp( Math.ceil(geodesic.surfaceDistance / 75000), 2, 256, @@ -313,9 +309,9 @@ export function createPaths({ state: layerState, services, parts, source }) { const eased = fraction * fraction * (3 - 2 * fraction); const cartographic = geodesic.interpolateUsingFraction( fraction, - new Cesium.Cartographic(), + geo.cartographic.create(), ); - cartographic.height = Cesium.Math.lerp(start.height, end.height, eased); + cartographic.height = geo.math.lerp(start.height, end.height, eased); positions.push(ellipsoid.cartographicToCartesian(cartographic)); } positions[0] = startPosition; @@ -340,24 +336,14 @@ export function createPaths({ state: layerState, services, parts, source }) { function blendAscentIntoOrbitTangent(ascentPath, orbitPath, insertionIndex) { if (!ascentPath?.length || ascentPath.length < 4 || !orbitPath?.length) return ascentPath; - const ellipsoid = Cesium.Ellipsoid.WGS84; + const ellipsoid = geo.ellipsoid.WGS84; const transferEnd = orbitPath[insertionIndex]; const nextOrbit = orbitPath[(insertionIndex + 1) % orbitPath.length]; - const orbitTangent = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.subtract( - nextOrbit, - transferEnd, - new Cesium.Cartesian3(), - ), - new Cesium.Cartesian3(), + const orbitTangent = geo.cartesian3.normalize( + geo.cartesian3.subtract(nextOrbit, transferEnd, geo.cartesian3.create()), + geo.cartesian3.create(), ); - if ( - Cesium.Cartesian3.equalsEpsilon( - orbitTangent, - Cesium.Cartesian3.ZERO, - 1e-8, - ) - ) + if (geo.cartesian3.equalsEpsilon(orbitTangent, geo.cartesian3.ZERO, 1e-8)) return ascentPath; // Replace a substantial final section with one cubic Bézier transition. @@ -372,62 +358,62 @@ export function createPaths({ state: layerState, services, parts, source }) { .map((position) => Math.max(0, ellipsoid.cartesianToCartographic(position)?.height || 0), ); - const ascentTangent = Cesium.Cartesian3.normalize( - Cesium.Cartesian3.subtract(start, previous, new Cesium.Cartesian3()), - new Cesium.Cartesian3(), + const ascentTangent = geo.cartesian3.normalize( + geo.cartesian3.subtract(start, previous, geo.cartesian3.create()), + geo.cartesian3.create(), ); const chordLength = Math.max( - Cesium.Cartesian3.distance(start, transferEnd), + geo.cartesian3.distance(start, transferEnd), 1000, ); const handleLength = chordLength * 0.28; - const controlA = Cesium.Cartesian3.add( + const controlA = geo.cartesian3.add( start, - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( ascentTangent, handleLength, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); - const controlB = Cesium.Cartesian3.add( + const controlB = geo.cartesian3.add( transferEnd, - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( orbitTangent, -handleLength, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), - new Cesium.Cartesian3(), + geo.cartesian3.create(), ); for (let index = 0; index <= blendCount; index++) { const progress = index / blendCount; const inverse = 1 - progress; - const point = new Cesium.Cartesian3(); - Cesium.Cartesian3.multiplyByScalar(start, inverse ** 3, point); - Cesium.Cartesian3.add( + const point = geo.cartesian3.create(); + geo.cartesian3.multiplyByScalar(start, inverse ** 3, point); + geo.cartesian3.add( point, - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( controlA, 3 * inverse ** 2 * progress, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), point, ); - Cesium.Cartesian3.add( + geo.cartesian3.add( point, - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( controlB, 3 * inverse * progress ** 2, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), point, ); - Cesium.Cartesian3.add( + geo.cartesian3.add( point, - Cesium.Cartesian3.multiplyByScalar( + geo.cartesian3.multiplyByScalar( transferEnd, progress ** 3, - new Cesium.Cartesian3(), + geo.cartesian3.create(), ), point, ); @@ -498,11 +484,11 @@ export function createPaths({ state: layerState, services, parts, source }) { return { lat: launch.lat, lon: launch.lon, accuracy: 'PAD / RTLS' }; } if (!(stage.downrangeKm > 0) || !insertionPosition) return null; - const ellipsoid = Cesium.Ellipsoid.WGS84; - const start = Cesium.Cartographic.fromDegrees(launch.lon, launch.lat); + const ellipsoid = geo.ellipsoid.WGS84; + const start = geo.cartographic.fromDegrees(launch.lon, launch.lat); const insertion = ellipsoid.cartesianToCartographic(insertionPosition); if (!insertion) return null; - const geodesic = new Cesium.EllipsoidGeodesic(start, insertion, ellipsoid); + const geodesic = geo.ellipsoidGeodesic.create(start, insertion, ellipsoid); const bearing = geodesic.startHeading; const angularDistance = (stage.downrangeKm * 1000) / ellipsoid.maximumRadius; @@ -519,8 +505,8 @@ export function createPaths({ state: layerState, services, parts, source }) { Math.cos(angularDistance) - Math.sin(lat1) * Math.sin(lat), ); return { - lat: Cesium.Math.toDegrees(lat), - lon: Cesium.Math.toDegrees(Cesium.Math.negativePiToPi(lon)), + lat: geo.math.toDegrees(lat), + lon: geo.math.toDegrees(geo.math.negativePiToPi(lon)), accuracy: 'EST. DOWNRANGE', }; } @@ -532,13 +518,13 @@ export function createPaths({ state: layerState, services, parts, source }) { stageCount, ) { if (!ascentPath?.length || !endpoint) return []; - const progress = Cesium.Math.clamp( + const progress = geo.math.clamp( 0.28 + (stageIndex / Math.max(stageCount, 1)) * 0.34, 0.28, 0.68, ); const separation = samplePath(ascentPath, progress); - const destination = Cesium.Cartesian3.fromDegrees( + const destination = geo.cartesian3.fromDegrees( endpoint.lon, endpoint.lat, 12, @@ -548,7 +534,7 @@ export function createPaths({ state: layerState, services, parts, source }) { function atmosphericReentryIndex(path) { if (!path?.length) return 0; - const ellipsoid = Cesium.Ellipsoid.WGS84; + const ellipsoid = geo.ellipsoid.WGS84; for (let index = 1; index < path.length; index++) { const previousHeight = ellipsoid.cartesianToCartographic( path[index - 1],