From 666eac3262f4bbdf2b0534866e63d10c58ce11bc Mon Sep 17 00:00:00 2001 From: Erin Millard Date: Tue, 28 May 2024 17:36:59 +1000 Subject: [PATCH] WIP --- src/coords.ts | 12 ++++++------ src/ecef.ts | 6 +++--- src/ellipsoid.ts | 22 +++++++++++++--------- src/euler.ts | 6 +++--- src/index.ts | 2 +- src/rotation-matrix.ts | 10 +++++----- test/arbitrary.ts | 5 ++--- test/vitest/ecef.spec.ts | 2 +- test/vitest/ellipsoid.spec.ts | 13 +++++++++---- 9 files changed, 43 insertions(+), 35 deletions(-) diff --git a/src/coords.ts b/src/coords.ts index 14d5200..54e7d79 100644 --- a/src/coords.ts +++ b/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]; } diff --git a/src/ecef.ts b/src/ecef.ts index bcb4f3f..45a7f9e 100644 --- a/src/ecef.ts +++ b/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); diff --git a/src/ellipsoid.ts b/src/ellipsoid.ts index b03f2b1..4ca701e 100644 --- a/src/ellipsoid.ts +++ b/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; diff --git a/src/euler.ts b/src/euler.ts index 27de7da..e834ac7 100644 --- a/src/euler.ts +++ b/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); } diff --git a/src/index.ts b/src/index.ts index 9c84d8a..849ecf4 100644 --- a/src/index.ts +++ b/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, diff --git a/src/rotation-matrix.ts b/src/rotation-matrix.ts index b4617ca..2e45e6f 100644 --- a/src/rotation-matrix.ts +++ b/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( diff --git a/test/arbitrary.ts b/test/arbitrary.ts index 7a25bc6..bdb4be9 100644 --- a/test/arbitrary.ts +++ b/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 { export function arbitraryEllipsoid(): fc.Arbitrary { return fc.oneof( fc.constant(WGS_84), - fc.constant(WGS_84_SPHERE), - fc.constant(WGS_72), fc.constant(GRS_80), + fc.constant(WGS_72), ); } diff --git a/test/vitest/ecef.spec.ts b/test/vitest/ecef.spec.ts index 91a789f..5a6dba2 100644 --- a/test/vitest/ecef.spec.ts +++ b/test/vitest/ecef.spec.ts @@ -99,7 +99,7 @@ describe("fromECEF()", () => { ); }); -describe("n_EB_E2p_EB_E()", () => { +describe("toECEF()", () => { let nvectorTestClient: NvectorTestClient; beforeAll(async () => { diff --git a/test/vitest/ellipsoid.spec.ts b/test/vitest/ellipsoid.spec.ts index d5fdae0..f216738 100644 --- a/test/vitest/ellipsoid.spec.ts +++ b/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); }); });