refactor(launches): geodesics on src/geo

Port native Vincenty geodesics and launch-path math with exact Cesium oracles.
Keep all existing launch path-test expectations unchanged.
This commit is contained in:
German Khokhlov committed 2026-10-22 23:48:50 +04:00
1 parent 6bdb995fe3
commit 1839d62160
7 files changed
+751 -113

No files matched your search

+2 -1
View File
@@ -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"
]
+24
View File
@@ -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",
+16
View File
@@ -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);
}
+527
View File
@@ -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);
}
+82
View File
@@ -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' },
);
});
+2
View File
@@ -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';
+98 -112
View File
@@ -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],