WIP

src/coords.ts

@@ -23,7 +23,7 @@ export function fromGeodeticCoordinates(
   // Equation (3) from Gade (2010):
   const cosLat = Math.cos(latitude);
 
-  // R_Ee selects correct E-axes
+  // frame selects correct E-axes
   return transform(transpose(frame), [
     Math.sin(latitude),
     Math.sin(longitude) * cosLat,
@@ -51,13 +51,13 @@ export function toGeodeticCoordinates(
 
   // Equation (6) in Gade (2010) (Robust numerical solution)
   // vector component in the equatorial plane
-  const equatorial_component = Math.hypot(y, z);
+  const ec = Math.hypot(y, z);
   // atan() could also be used since latitude is within [-pi/2,pi/2]
-  const latitude = Math.atan2(x, equatorial_component);
+  const latitude = Math.atan2(x, ec);
 
-  // latitude = asin(n_E(1)) is a theoretical solution, but close to the Poles
-  // it is ill-conditioned which may lead to numerical inaccuracies (and it will
-  // give imaginary results for norm(n_E)>1)
+  // latitude = asin(x) is a theoretical solution, but close to the Poles it is
+  // ill-conditioned which may lead to numerical inaccuracies (and it will give
+  // imaginary results for norm(vector)>1)
 
   return [latitude, longitude];
 }

src/ecef.ts

@@ -21,7 +21,7 @@ export function fromECEF(
   vector: Vector,
   { a, f }: Ellipsoid = WGS_84,
   frame: Matrix = Z_AXIS_NORTH,
-): [n_EB_E: Vector, depth: number] {
+): [vector: Vector, depth: number] {
   // frame selects correct E-axes
   const [x, y, z] = transform(frame, vector);
 
@@ -94,8 +94,8 @@ export function toECEF(
   const denom = Math.hypot(x, y / (1 - f), z / (1 - f));
 
   // We first calculate the position at the origin of coordinate system L, which
-  // has the same n-vector as B (n_EL_E = n_EB_E), but lies at the surface of
-  // the Earth (depth = 0).
+  // has the same n-vector as B (bEL = bEB), but lies at the surface of the
+  // Earth (depth = 0).
 
   const ex = (b / denom) * x;
   const ey = (b / denom) * (y / (1 - f) ** 2);

src/ellipsoid.ts

@@ -32,27 +32,31 @@ export const WGS_84: Ellipsoid = {
 } as const;
 
 /**
- * A sphere with the same semi-major axis as the WGS-84 ellipsoid.
+ * Create a spherical ellipsoid with the given radius.
  *
- * @see https://github.com/chrisveness/geodesy/blob/761587cd748bd9f7c9825195eba4a9fc5891b859/latlon-ellipsoidal-datum.js#L39
+ * @param radius - A radius in meters.
+ *
+ * @returns A spherical ellipsoid.
  */
-export const WGS_84_SPHERE: Ellipsoid = {
-  a: WGS_84.a,
-  b: WGS_84.a,
-  f: 0,
-} as const;
+export function sphere(radius: number): Ellipsoid {
+  return {
+    a: radius,
+    b: radius,
+    f: 0,
+  };
+}
 
 /**
  * An ellipsoid.
  */
 export type Ellipsoid = {
   /**
-   * The semi-major axis of the ellipsoid.
+   * The semi-major axis of the ellipsoid in meters.
    */
   a: number;
 
   /**
-   * The semi-minor axis of the ellipsoid.
+   * The semi-minor axis of the ellipsoid in meters.
    */
   b: number;
 

src/euler.ts

@@ -102,7 +102,7 @@ export function rotationMatrixToEulerXYZ([
     // Only the sum/difference of x and z is now given, choosing x = 0:
     x = 0;
 
-    // Lower left 2x2 elements of R_AB now only consists of sin_z and cos_z.
+    // Lower left 2x2 elements of rotation now only consists of sin z and cos z.
     // Using the two whose signs are the same for both singularities:
     z = Math.atan2(r21, r22);
   }
@@ -148,8 +148,8 @@ export function rotationMatrixToEulerZYX([
     // Only the sum/difference of x and z is now given, choosing x = 0:
     x = 0;
 
-    // Upper right 2x2 elements of R_AB now only consists of sin_z and cos_z.
-    // Using the two whose signs are the same for both singularities:
+    // Upper right 2x2 elements of rotation now only consists of sin z and cos
+    // z. Using the two whose signs are the same for both singularities:
     z = Math.atan2(-r12, r22);
   }
 

src/index.ts

@@ -3,7 +3,7 @@ export { X_AXIS_NORTH, Z_AXIS_NORTH } from "./coord-frame.js";
 export { fromGeodeticCoordinates, toGeodeticCoordinates } from "./coords.js";
 export { delta, destination } from "./delta.js";
 export { fromECEF, toECEF } from "./ecef.js";
-export { GRS_80, WGS_72, WGS_84, WGS_84_SPHERE } from "./ellipsoid.js";
+export { GRS_80, WGS_72, WGS_84, sphere } from "./ellipsoid.js";
 export type { Ellipsoid } from "./ellipsoid.js";
 export {
   eulerXYZToRotationMatrix,

src/rotation-matrix.ts

@@ -38,7 +38,7 @@ export function toRotationMatrix(
 
   // N coordinate frame (North-East-Down) is defined in Table 2 in Gade (2010)
 
-  // R_EN is constructed by the following three column vectors: The x, y and z
+  // rEN is constructed by the following three column vectors: The x, y and z
   // basis vectors (axes) of N, each decomposed in E.
 
   // Find z-axis of N (Nz):
@@ -55,7 +55,7 @@ export function toRotationMatrix(
   const yzDir = y;
   const yDirNorm = Math.hypot(yyDir, yzDir);
   const onPoles = Math.hypot(yyDir, yzDir) === 0;
-  // Ny_E_x is always 0, so it's factored out in the following equations
+  // yx is always 0, so it's factored out in the following equations
   const yy = onPoles
     ? 1 // Pole position: selected y-axis direction
     : yyDir / yDirNorm; // outside Poles:
@@ -69,7 +69,7 @@ export function toRotationMatrix(
   const xy = yz * zx;
   const xz = -yy * zx;
 
-  // Form R_EN from the unit vectors:
+  // Form rotation from the unit vectors:
   // frame selects correct E-axes
   return multiply(transpose(frame), [
     [xx, 0, zx],
@@ -97,8 +97,8 @@ export function toRotationMatrixUsingWanderAzimuth(
   const [latitude, longitude] = toGeodeticCoordinates(vector, frame);
 
   // Longitude, -latitude, and wander azimuth are the x-y-z Euler angles (about
-  // new axes) for R_EL. See also the second paragraph of Section 5.2 in Gade
-  // (2010):
+  // new axes) for rotation. See also the second paragraph of Section 5.2 in
+  // Gade (2010):
 
   // frame selects correct E-axes
   return multiply(

test/arbitrary.ts

@@ -1,6 +1,6 @@
 import { fc } from "@fast-check/vitest";
 import type { Ellipsoid, Matrix, Vector } from "../src/index.js";
-import { GRS_80, WGS_72, WGS_84, WGS_84_SPHERE } from "../src/index.js";
+import { GRS_80, WGS_72, WGS_84 } from "../src/index.js";
 
 const RADIAN = Math.PI / 180;
 
@@ -52,9 +52,8 @@ export function arbitrary3dUnitVector(): fc.Arbitrary<Vector> {
 export function arbitraryEllipsoid(): fc.Arbitrary<Ellipsoid> {
   return fc.oneof(
     fc.constant(WGS_84),
-    fc.constant(WGS_84_SPHERE),
-    fc.constant(WGS_72),
     fc.constant(GRS_80),
+    fc.constant(WGS_72),
   );
 }
 

test/vitest/ecef.spec.ts

@@ -99,7 +99,7 @@ describe("fromECEF()", () => {
   );
 });
 
-describe("n_EB_E2p_EB_E()", () => {
+describe("toECEF()", () => {
   let nvectorTestClient: NvectorTestClient;
 
   beforeAll(async () => {

test/vitest/ellipsoid.spec.ts

@@ -1,6 +1,6 @@
 import { it } from "@fast-check/vitest";
 import { describe, expect } from "vitest";
-import { GRS_80, WGS_72, WGS_84, WGS_84_SPHERE } from "../../src/index.js";
+import { GRS_80, WGS_72, WGS_84, sphere } from "../../src/index.js";
 
 describe("GRS_80", () => {
   it("has the correct semi-minor axis", () => {
@@ -20,8 +20,13 @@ describe("WGS_84", () => {
   });
 });
 
-describe("WGS_84_SPHERE", () => {
-  it("has the correct semi-minor axis", () => {
-    expect(WGS_84_SPHERE.b).toEqual(WGS_84_SPHERE.a);
+describe("sphere", () => {
+  it("produces a sphere with the given radius", () => {
+    const radius = 6371e3;
+    const ellipsoid = sphere(radius);
+
+    expect(ellipsoid.a).toEqual(radius);
+    expect(ellipsoid.b).toEqual(radius);
+    expect(ellipsoid.f).toEqual(0);
   });
 });