diff --git a/src/geo/ellipsoid.js b/src/geo/ellipsoid.js index 45d67da0d..dbff51e3f 100644 --- a/src/geo/ellipsoid.js +++ b/src/geo/ellipsoid.js @@ -4,10 +4,20 @@ * (ECEF) conversions, surface normals and projection onto the surface. The * formulas follow the default renderer's step for step, so they agree with * it to floating-point noise. + * + * Into ECEF, positions are arguments (`fromDegrees(lon, lat, height)`, + * `fromRadians` likewise); out of it, objects (`toLonLatH` gives + * `{lon, lat, height}` in degrees, `toCartographicRadians` the same in + * radians). Cesium's `Cartographic.fromDegrees(lon, lat, h)` is plain + * `{lon: toRadians(lon), lat: toRadians(lat), height: h}`; for its ECEF + * point, use `fromDegrees`. + * + * NaN in gives NaN fields out, never a throw (Cesium's debug build throws). + * The Earth's centre, and points so near it that their squares round to + * zero, answer null, as Cesium answers undefined. */ import { add, - clone, dot, magnitude, normalize, @@ -58,7 +68,7 @@ export function fromRadians(lon, lat, height = 0) { const r = WGS84.radiiSquared; let k = vec(r.x * n.x, r.y * n.y, r.z * n.z); const gamma = Math.sqrt(dot(n, k)); - k = scale(k, 1 / gamma); + k = vec(k.x / gamma, k.y / gamma, k.z / gamma); return add(k, scale(n, height)); } @@ -66,13 +76,6 @@ export function fromRadians(lon, lat, height = 0) { export const fromDegrees = (lon, lat, height = 0) => fromRadians(toRadians(lon), toRadians(lat), height); -/** A position in either contract form, as ECEF. */ -export function toEcef(position) { - if ('lon' in position && 'lat' in position) - return fromDegrees(position.lon, position.lat, position.height ?? 0); - return clone(position); -} - /** The point on the ellipsoid along the geodetic normal, or null at the centre. */ export function scaleToGeodeticSurface(p) { const { x: positionX, y: positionY, z: positionZ } = p; diff --git a/src/geo/geo.test.mjs b/src/geo/geo.test.mjs index 9e0ad1a58..e0fcde931 100644 --- a/src/geo/geo.test.mjs +++ b/src/geo/geo.test.mjs @@ -29,12 +29,20 @@ function points(count = 500) { { lon: 180, lat: 0, height: -400 }, { lon: -97.7431, lat: 30.2672, height: 150 }, ]; - for (let i = 0; i < count; i++) + // Each height band gets its share: near the ground, low orbit, far out. + const bands = [ + [-500, 1e4], + [1e4, 1e6], + [1e6, 4e7], + ]; + for (let i = 0; i < count; i++) { + const [lo, hi] = bands[i % bands.length]; out.push({ lon: next() * 360 - 180, lat: next() * 180 - 90, - height: next() * 4e7 - 1e4, + height: lo + next() * (hi - lo), }); + } return out; } @@ -60,7 +68,7 @@ test('geodetic to ECEF matches the default renderer', () => { for (const p of points()) { const ours = fromDegrees(p.lon, p.lat, p.height); const theirs = Cesium.Cartesian3.fromDegrees(p.lon, p.lat, p.height); - closeVec(ours, theirs, 1e-6, `${p.lon},${p.lat},${p.height}`); + closeVec(ours, theirs, 0, `${p.lon},${p.lat},${p.height}`); } }); @@ -76,6 +84,37 @@ test('ECEF to geodetic matches the default renderer', () => { close(ours.lat, theirs.latitude, 1e-12, `${where} lat`); close(ours.height, theirs.height, 1e-6, `${where} height`); } + // Inside the Earth: within the centre cutoff (the first four) and between + // it and the surface (the last). + for (const p of [ + { x: 1e6, y: 0, z: 0 }, + { x: 1, y: 2, z: 3 }, + { x: -2e5, y: 3e5, z: 1e5 }, + { x: 1e6, y: 1e6, z: 5e5 }, + { x: 2e6, y: 1.5e6, z: 1.5e6 }, + ]) { + const ours = toCartographicRadians(p); + const theirs = Cesium.Cartographic.fromCartesian( + new Cesium.Cartesian3(p.x, p.y, p.z), + ); + const where = `${p.x},${p.y},${p.z}`; + close(ours.lon, theirs.longitude, 1e-12, `${where} lon`); + close(ours.lat, theirs.latitude, 1e-12, `${where} lat`); + close(ours.height, theirs.height, 1e-6, `${where} height`); + } + for (const p of [ + { x: NaN, y: 0, z: 0 }, + { x: 7e6, y: NaN, z: 0 }, + { x: 7e6, y: 0, z: NaN }, + ]) { + const nan = toLonLatH(p); + assert.ok( + Number.isNaN(nan.lon) && + Number.isNaN(nan.lat) && + Number.isNaN(nan.height), + JSON.stringify(nan), + ); + } assert.equal(toCartographicRadians({ x: 0, y: 0, z: 0 }), null); assert.equal( Cesium.Cartographic.fromCartesian(Cesium.Cartesian3.ZERO), diff --git a/src/geo/vector.js b/src/geo/vector.js index e020dd5cc..a5b7400fc 100644 --- a/src/geo/vector.js +++ b/src/geo/vector.js @@ -18,6 +18,7 @@ export const dot = (a, b) => a.x * b.x + a.y * b.y + a.z * b.z; export const cross = (a, b) => vec(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); export const magnitude = (a) => Math.sqrt(dot(a, a)); +/** Straight-line (chord) distance, not along the surface. */ export const distance = (a, b) => magnitude(sub(a, b)); export function normalize(a) { const m = magnitude(a); diff --git a/src/render/stub/math.js b/src/render/stub/math.js index 4fa6665b1..896309eca 100644 --- a/src/render/stub/math.js +++ b/src/render/stub/math.js @@ -14,6 +14,8 @@ import { clamp, equalsEpsilon, negate, + clone, + fromDegrees, rotate, scale, vec, @@ -22,6 +24,13 @@ import { export * from '../../geo/index.js'; +/** A contract `Position` (`Ecef` or `LonLatH`, degrees), as ECEF. */ +export function toEcef(position) { + if ('lon' in position && 'lat' in position) + return fromDegrees(position.lon, position.lat, position.height ?? 0); + return clone(position); +} + const EPSILON3 = 1e-3; const EPSILON10 = 1e-10;