WIP

test/vitest/examples/example-01.spec.ts

@@ -11,61 +11,63 @@ import {
 /**
  * Example 1: A and B to delta
  *
- * Given two positions A and B. Find the exact vector from A to B in meters
- * north, east and down, and find the direction (azimuth/bearing) to B, relative
- * to north. Use WGS-84 ellipsoid.
- *
- * - Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface.
- * - Use position A to define north, east, and down directions. (Due to the
- *   curvature of Earth and different directions to the North Pole, the north,
- *   east, and down directions will change (relative to Earth) for different
- *   places. Position A must be outside the poles for the north and east
- *   directions to be defined.
- *
  * @see https://www.ffi.no/en/research/n-vector/#example_1
  */
 test("Example 1", () => {
-  // Positions A and B are given in (decimal) degrees and depths:
-
-  // Position A:
-  const aLat = 1;
-  const aLon = 2;
-  const aDepth = 3;
+  // Given two positions, A and B as latitudes, longitudes and depths (relative
+  // to Earth, E):
 
-  // Position B:
-  const bLat = 4;
-  const bLon = 5;
-  const bDepth = 6;
+  const aLat = 1,
+    aLon = 2,
+    aDepth = 3;
+  const bLat = 4,
+    bLon = 5,
+    bDepth = 6;
 
   // Find the exact vector between the two positions, given in meters north,
-  // east, and down, i.e. find p_AB_N.
+  // east, and down, and find the direction (azimuth) to B, relative to north.
+  //
+  // Details:
+  //
+  // - Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface.
+  // - Use position A to define north, east, and down directions. (Due to the
+  //   curvature of Earth and different directions to the North Pole, the north,
+  //   east, and down directions will change (relative to Earth) for different
+  //   places. Position A must be outside the poles for the north and east
+  //   directions to be defined.
 
   // SOLUTION:
 
-  // Step1: Convert to n-vectors (rad() converts to radians):
+  // Step 1: First, the given latitudes and longitudes are converted to
+  //         n-vectors:
   const a = fromGeodeticCoordinates(radians(aLat), radians(aLon));
   const b = fromGeodeticCoordinates(radians(bLat), radians(bLon));
 
-  // Step2: Find p_AB_E (delta decomposed in E). WGS-84 ellipsoid is default:
-  const pe = delta(a, b, aDepth, bDepth);
+  // Step 2: When the positions are given as n-vectors (and depths), it is easy
+  //         to find the delta vector decomposed in E. No ellipsoid is specified
+  //         when calling the function, thus WGS-84 (default) is used:
+  const abE = delta(a, b, aDepth, bDepth);
 
-  // Step3: Find R_EN for position A:
-  const r = toRotationMatrix(a);
+  // Step 3: We now have the delta vector from A to B, but the three coordinates
+  //         of the vector are along the Earth coordinate frame E, while we need
+  //         the coordinates to be north, east and down. To get this, we define
+  //         a North-East-Down coordinate frame called N, and then we need the
+  //         rotation matrix (direction cosine matrix) rEN to go between E and
+  //         N. We have a simple function that calculates rEN from an n-vector,
+  //         and we use this function (using the n-vector at position A):
+  const rEN = toRotationMatrix(a);
 
-  // Step4: Find p_AB_N
-  const pn = transform(transpose(r), pe);
-  // (Note the transpose of R_EN: The "closest-rule" says that when decomposing,
-  // the frame in the subscript of the rotation matrix that is closest to the
-  // vector, should equal the frame where the vector is decomposed. Thus the
-  // calculation R_NE*p_AB_E is correct, since the vector is decomposed in E,
-  // and E is closest to the vector. In the above example we only had R_EN, and
-  // thus we must transpose it: R_EN' = R_NE)
+  // Step 4: Now the delta vector is easily decomposed in N. Since the vector is
+  //         decomposed in E, we must use rNE (rNE is the transpose of rEN):
+  const abN = transform(transpose(rEN), abE);
 
-  // Step5: Also find the direction (azimuth) to B, relative to north:
-  const azimuth = Math.atan2(pn[1], pn[0]);
+  // Step 5: The three components of abN are the north, east and down
+  //         displacements from A to B in meters. The azimuth is simply found
+  //         from element 1 and 2 of the vector (the north and east components):
+  const azimuth = Math.atan2(abN[1], abN[0]);
 
-  expect(pn[0]).toBeCloseTo(331730.2347808944, 8); // meters
-  expect(pn[1]).toBeCloseTo(332997.8749892695, 8); // meters
-  expect(pn[2]).toBeCloseTo(17404.27136193635, 8); // meters
+  expect(abN[0]).toBeCloseTo(331730.2347808944, 8); // meters
+  expect(abN[1]).toBeCloseTo(332997.8749892695, 8); // meters
+  expect(abN[2]).toBeCloseTo(17404.27136193635, 8); // meters
   expect(azimuth).toBeCloseTo(radians(45.10926323826139), 15);
 });