Match Cesium's fromDegrees exactly and keep toEcef in the stub

fromRadians divides by gamma as Cesium does, so geodetic to ECEF is now
bit for bit Cesium's, and the test checks it exactly. Its tests also
cover three height bands, points inside the Earth and NaN input. toEcef
reads the renderer contract's Position, so it moves to the stub's maths.
The ellipsoid docs pair each conversion with Cesium's.

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;

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),

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);

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;