diff --git a/README.md b/README.md
index 76a5140..5f68e99 100644
--- a/README.md
+++ b/README.md
@@ -103,8 +103,8 @@ test("Example 1", () => {
   // 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));
+  const a = fromGeodeticCoordinates(radians(aLon), radians(aLat));
+  const b = fromGeodeticCoordinates(radians(bLon), radians(bLat));
 
   // Step 2
   //
@@ -228,7 +228,7 @@ test("Example 2", () => {
   const [c, cDepth] = destination(b, bcE, bDepth, e);
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(c);
+  const [lon, lat] = toGeodeticCoordinates(c);
   const height = -cDepth;
 
   expect(degrees(lat)).toBeCloseTo(53.32637826433107, 13);
@@ -284,7 +284,7 @@ test("Example 3", () => {
   // Step 2
   //
   // Find latitude, longitude and height:
-  const [lat, lon] = toGeodeticCoordinates(b);
+  const [lon, lat] = toGeodeticCoordinates(b);
   const height = -bDepth;
 
   expect(degrees(lat)).toBeCloseTo(5.685075734513181, 14);
@@ -327,7 +327,7 @@ test("Example 4", () => {
   // SOLUTION:
 
   // Step 1: First, the given latitude and longitude are converted to n-vector:
-  const b = fromGeodeticCoordinates(radians(bLat), radians(bLon));
+  const b = fromGeodeticCoordinates(radians(bLon), radians(bLat));
 
   // Step 2: Convert to an ECEF-vector:
   const pb = toECEF(b, -bHeight);
@@ -370,8 +370,8 @@ test("Example 5", () => {
   // PROBLEM:
 
   // Given two positions A and B as n-vectors:
-  const a = fromGeodeticCoordinates(radians(88), radians(0));
-  const b = fromGeodeticCoordinates(radians(89), radians(-170));
+  const a = fromGeodeticCoordinates(radians(0), radians(88));
+  const b = fromGeodeticCoordinates(radians(-170), radians(89));
 
   // Find the surface distance (i.e. great circle distance). The heights of A
   // and B are not relevant (i.e. if they do not have zero height, we seek the
@@ -430,8 +430,8 @@ test("Example 6", () => {
   const t0 = 10,
     t1 = 20,
     ti = 16;
-  const pt0 = fromGeodeticCoordinates(radians(89.9), radians(-150));
-  const pt1 = fromGeodeticCoordinates(radians(89.9), radians(150));
+  const pt0 = fromGeodeticCoordinates(radians(-150), radians(89.9));
+  const pt1 = fromGeodeticCoordinates(radians(150), radians(89.9));
 
   // Find an interpolated position at time ti, pti. All positions are given as
   // n-vectors.
@@ -444,7 +444,7 @@ test("Example 6", () => {
   );
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(pti);
+  const [lon, lat] = toGeodeticCoordinates(pti);
 
   expect(degrees(lat)).toBeCloseTo(89.91282199988446, 12);
   expect(degrees(lon)).toBeCloseTo(173.4132244463705, 12);
@@ -480,9 +480,9 @@ test("Example 7", () => {
   // PROBLEM:
 
   // Three positions A, B, and C are given as n-vectors:
-  const a = fromGeodeticCoordinates(radians(90), radians(0));
-  const b = fromGeodeticCoordinates(radians(60), radians(10));
-  const c = fromGeodeticCoordinates(radians(50), radians(-20));
+  const a = fromGeodeticCoordinates(radians(0), radians(90));
+  const b = fromGeodeticCoordinates(radians(10), radians(60));
+  const c = fromGeodeticCoordinates(radians(-20), radians(50));
 
   // Find the mean position, M. Note that the calculation is independent of the
   // heights/depths of the positions.
@@ -534,7 +534,7 @@ test("Example 8", () => {
   // PROBLEM:
 
   // Position A is given as n-vector:
-  const a = fromGeodeticCoordinates(radians(80), radians(-90));
+  const a = fromGeodeticCoordinates(radians(-90), radians(80));
 
   // We also have an initial direction of travel given as an azimuth (bearing)
   // relative to north (clockwise), and finally the distance to travel along a
@@ -587,7 +587,7 @@ test("Example 8", () => {
   );
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(b);
+  const [lon, lat] = toGeodeticCoordinates(b);
 
   expect(degrees(lat)).toBeCloseTo(79.99154867339445, 13);
   expect(degrees(lon)).toBeCloseTo(-90.01769837291397, 13);
@@ -633,12 +633,12 @@ test("Example 9", () => {
   // the two positions are not antipodal).
 
   // Path A is given by a1 and a2:
-  const a1 = fromGeodeticCoordinates(radians(50), radians(180));
-  const a2 = fromGeodeticCoordinates(radians(90), radians(180));
+  const a1 = fromGeodeticCoordinates(radians(180), radians(50));
+  const a2 = fromGeodeticCoordinates(radians(180), radians(90));
 
   // While path B is given by b1 and b2:
-  const b1 = fromGeodeticCoordinates(radians(60), radians(160));
-  const b2 = fromGeodeticCoordinates(radians(80), radians(-140));
+  const b1 = fromGeodeticCoordinates(radians(160), radians(60));
+  const b2 = fromGeodeticCoordinates(radians(-140), radians(80));
 
   // Find the position C where the two paths intersect.
 
@@ -660,7 +660,7 @@ test("Example 9", () => {
   const c = apply((n) => Math.sign(dot(cTmp, a1)) * n, cTmp);
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(c);
+  const [lon, lat] = toGeodeticCoordinates(c);
 
   expect(degrees(lat)).toBeCloseTo(74.16344802135536, 16);
   expect(degrees(lon)).toBeCloseTo(180, 16);
@@ -701,10 +701,10 @@ test("Example 10", () => {
   // Path A is given by the two n-vectors a1 and a2 (as in the previous
   // example):
   const a1 = fromGeodeticCoordinates(radians(0), radians(0));
-  const a2 = fromGeodeticCoordinates(radians(10), radians(0));
+  const a2 = fromGeodeticCoordinates(radians(0), radians(10));
 
   // And a position B is given by b:
-  const b = fromGeodeticCoordinates(radians(1), radians(0.1));
+  const b = fromGeodeticCoordinates(radians(0.1), radians(1));
 
   // Find the cross track distance between the path A (i.e. the great circle
   // through a1 and a2) and the position B (i.e. the shortest distance at the
diff --git a/src/coords.ts b/src/coords.ts
index 54e7d79..3eeef9e 100644
--- a/src/coords.ts
+++ b/src/coords.ts
@@ -9,15 +9,15 @@ import { transform } from "./vector.js";
  *
  * @see https://github.com/FFI-no/n-vector/blob/82d749a67cc9f332f48c51aa969cdc277b4199f2/nvector/lat_long2n_E.m
  *
- * @param latitude - Geodetic latitude in radians.
  * @param longitude - Geodetic longitude in radians.
+ * @param latitude - Geodetic latitude in radians.
  * @param frame - Coordinate frame in which the n-vector is decomposed.
  *
  * @returns An n-vector.
  */
 export function fromGeodeticCoordinates(
-  latitude: number,
   longitude: number,
+  latitude: number,
   frame: Matrix = Z_AXIS_NORTH,
 ): Vector {
   // Equation (3) from Gade (2010):
@@ -39,12 +39,12 @@ export function fromGeodeticCoordinates(
  * @param vector - An n-vector.
  * @param frame - Coordinate frame in which the n-vector is decomposed.
  *
- * @returns Geodetic latitude and longitude in radians.
+ * @returns Geodetic longitude and latitude in radians.
  */
 export function toGeodeticCoordinates(
   vector: Vector,
   frame: Matrix = Z_AXIS_NORTH,
-): [latitude: number, longitude: number] {
+): [longitude: number, latitude: number] {
   // Equation (5) in Gade (2010):
   const [x, y, z] = transform(frame, vector);
   const longitude = Math.atan2(y, -z);
@@ -59,5 +59,5 @@ export function toGeodeticCoordinates(
   // ill-conditioned which may lead to numerical inaccuracies (and it will give
   // imaginary results for norm(vector)>1)
 
-  return [latitude, longitude];
+  return [longitude, latitude];
 }
diff --git a/src/rotation-matrix.ts b/src/rotation-matrix.ts
index 2e45e6f..39512a7 100644
--- a/src/rotation-matrix.ts
+++ b/src/rotation-matrix.ts
@@ -94,7 +94,7 @@ export function toRotationMatrixUsingWanderAzimuth(
   wanderAzimuth: number,
   frame: Matrix = Z_AXIS_NORTH,
 ): Matrix {
-  const [latitude, longitude] = toGeodeticCoordinates(vector, frame);
+  const [longitude, latitude] = toGeodeticCoordinates(vector, frame);
 
   // Longitude, -latitude, and wander azimuth are the x-y-z Euler angles (about
   // new axes) for rotation. See also the second paragraph of Section 5.2 in
diff --git a/test/arbitrary.ts b/test/arbitrary.ts
index 2c300d8..407f53b 100644
--- a/test/arbitrary.ts
+++ b/test/arbitrary.ts
@@ -70,16 +70,16 @@ export function arbitraryEllipsoidECEFVector({
   return arbitrary3dVector({ min: a - b, max: a + b, noNaN: true });
 }
 
-export function arbitraryLatLon(): fc.Arbitrary<[number, number]> {
+export function arbitraryGeodeticCoordinates(): fc.Arbitrary<[number, number]> {
   return fc.tuple(
     fc.double({
-      min: -90 * RADIAN,
-      max: 90 * RADIAN,
+      min: -180 * RADIAN,
+      max: 180 * RADIAN,
       noNaN: true,
     }),
     fc.double({
-      min: -180 * RADIAN,
-      max: 180 * RADIAN,
+      min: -90 * RADIAN,
+      max: 90 * RADIAN,
       noNaN: true,
     }),
   );
diff --git a/test/nvector-test-api.ts b/test/nvector-test-api.ts
index 15736f9..062c785 100644
--- a/test/nvector-test-api.ts
+++ b/test/nvector-test-api.ts
@@ -47,11 +47,11 @@ export async function createNvectorTestClient(): Promise<NvectorTestClient> {
   });
 
   return {
-    async fromGeodeticCoordinates(latitude, longitude, frame) {
+    async fromGeodeticCoordinates(longitude, latitude, frame) {
       return unwrapVector3(
         await call<WrappedVector3>("lat_lon2n_E", {
-          latitude,
           longitude,
+          latitude,
           R_Ee: frame,
         }),
       );
@@ -66,15 +66,15 @@ export async function createNvectorTestClient(): Promise<NvectorTestClient> {
     },
 
     async toGeodeticCoordinates(vector, frame) {
-      const { latitude, longitude } = await call<{
-        latitude: number;
+      const { longitude, latitude } = await call<{
         longitude: number;
+        latitude: number;
       }>("n_E2lat_lon", {
         n_E: wrapVector3(vector),
         R_Ee: frame,
       });
 
-      return [latitude, longitude];
+      return [longitude, latitude];
     },
 
     async toRotationMatrix(vector, frame) {
diff --git a/test/vitest/coords.spec.ts b/test/vitest/coords.spec.ts
index 9118965..ce140a4 100644
--- a/test/vitest/coords.spec.ts
+++ b/test/vitest/coords.spec.ts
@@ -7,7 +7,7 @@ import { afterAll, beforeAll, describe, expect } from "vitest";
 import {
   arbitrary3dRotationMatrix,
   arbitrary3dUnitVector,
-  arbitraryLatLon,
+  arbitraryGeodeticCoordinates,
 } from "../arbitrary.js";
 import type { NvectorTestClient } from "../nvector-test-api.js";
 import { createNvectorTestClient } from "../nvector-test-api.js";
@@ -27,16 +27,16 @@ describe("fromGeodeticCoordinates()", () => {
 
   it.prop(
     [
-      arbitraryLatLon(),
+      arbitraryGeodeticCoordinates(),
       fc.option(arbitrary3dRotationMatrix(), { nil: undefined }),
     ],
     { interruptAfterTimeLimit: TEST_DURATION, numRuns: Infinity },
   )(
     "matches the reference implementation",
-    async ([latitude, longitude], frame) => {
+    async ([longitude, latitude], frame) => {
       const expected = await nvectorTestClient.fromGeodeticCoordinates(
-        latitude,
         longitude,
+        latitude,
         frame,
       );
 
@@ -46,7 +46,7 @@ describe("fromGeodeticCoordinates()", () => {
         expect.any(Number),
       ]);
 
-      const actual = fromGeodeticCoordinates(latitude, longitude, frame);
+      const actual = fromGeodeticCoordinates(longitude, latitude, frame);
 
       expect(actual).toMatchObject([
         expect.any(Number),
diff --git a/test/vitest/examples/example-01.spec.ts b/test/vitest/examples/example-01.spec.ts
index a5a6dc4..43736ed 100644
--- a/test/vitest/examples/example-01.spec.ts
+++ b/test/vitest/examples/example-01.spec.ts
@@ -46,8 +46,8 @@ test("Example 1", () => {
   // 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));
+  const a = fromGeodeticCoordinates(radians(aLon), radians(aLat));
+  const b = fromGeodeticCoordinates(radians(bLon), radians(bLat));
 
   // Step 2
   //
diff --git a/test/vitest/examples/example-02.spec.ts b/test/vitest/examples/example-02.spec.ts
index 5b60055..8254558 100644
--- a/test/vitest/examples/example-02.spec.ts
+++ b/test/vitest/examples/example-02.spec.ts
@@ -72,7 +72,7 @@ test("Example 2", () => {
   const [c, cDepth] = destination(b, bcE, bDepth, e);
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(c);
+  const [lon, lat] = toGeodeticCoordinates(c);
   const height = -cDepth;
 
   expect(degrees(lat)).toBeCloseTo(53.32637826433107, 13);
diff --git a/test/vitest/examples/example-03.spec.ts b/test/vitest/examples/example-03.spec.ts
index fd2791d..b8d6894 100644
--- a/test/vitest/examples/example-03.spec.ts
+++ b/test/vitest/examples/example-03.spec.ts
@@ -35,7 +35,7 @@ test("Example 3", () => {
   // Step 2
   //
   // Find latitude, longitude and height:
-  const [lat, lon] = toGeodeticCoordinates(b);
+  const [lon, lat] = toGeodeticCoordinates(b);
   const height = -bDepth;
 
   expect(degrees(lat)).toBeCloseTo(5.685075734513181, 14);
diff --git a/test/vitest/examples/example-04.spec.ts b/test/vitest/examples/example-04.spec.ts
index c86a8cc..05418aa 100644
--- a/test/vitest/examples/example-04.spec.ts
+++ b/test/vitest/examples/example-04.spec.ts
@@ -22,7 +22,7 @@ test("Example 4", () => {
   // SOLUTION:
 
   // Step 1: First, the given latitude and longitude are converted to n-vector:
-  const b = fromGeodeticCoordinates(radians(bLat), radians(bLon));
+  const b = fromGeodeticCoordinates(radians(bLon), radians(bLat));
 
   // Step 2: Convert to an ECEF-vector:
   const pb = toECEF(b, -bHeight);
diff --git a/test/vitest/examples/example-05.spec.ts b/test/vitest/examples/example-05.spec.ts
index 2d85615..5911eb7 100644
--- a/test/vitest/examples/example-05.spec.ts
+++ b/test/vitest/examples/example-05.spec.ts
@@ -20,8 +20,8 @@ test("Example 5", () => {
   // PROBLEM:
 
   // Given two positions A and B as n-vectors:
-  const a = fromGeodeticCoordinates(radians(88), radians(0));
-  const b = fromGeodeticCoordinates(radians(89), radians(-170));
+  const a = fromGeodeticCoordinates(radians(0), radians(88));
+  const b = fromGeodeticCoordinates(radians(-170), radians(89));
 
   // Find the surface distance (i.e. great circle distance). The heights of A
   // and B are not relevant (i.e. if they do not have zero height, we seek the
diff --git a/test/vitest/examples/example-06.spec.ts b/test/vitest/examples/example-06.spec.ts
index a7fb772..ea80082 100644
--- a/test/vitest/examples/example-06.spec.ts
+++ b/test/vitest/examples/example-06.spec.ts
@@ -23,8 +23,8 @@ test("Example 6", () => {
   const t0 = 10,
     t1 = 20,
     ti = 16;
-  const pt0 = fromGeodeticCoordinates(radians(89.9), radians(-150));
-  const pt1 = fromGeodeticCoordinates(radians(89.9), radians(150));
+  const pt0 = fromGeodeticCoordinates(radians(-150), radians(89.9));
+  const pt1 = fromGeodeticCoordinates(radians(150), radians(89.9));
 
   // Find an interpolated position at time ti, pti. All positions are given as
   // n-vectors.
@@ -37,7 +37,7 @@ test("Example 6", () => {
   );
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(pti);
+  const [lon, lat] = toGeodeticCoordinates(pti);
 
   expect(degrees(lat)).toBeCloseTo(89.91282199988446, 12);
   expect(degrees(lon)).toBeCloseTo(173.4132244463705, 12);
diff --git a/test/vitest/examples/example-07.spec.ts b/test/vitest/examples/example-07.spec.ts
index 3282d20..fddcfb3 100644
--- a/test/vitest/examples/example-07.spec.ts
+++ b/test/vitest/examples/example-07.spec.ts
@@ -17,9 +17,9 @@ test("Example 7", () => {
   // PROBLEM:
 
   // Three positions A, B, and C are given as n-vectors:
-  const a = fromGeodeticCoordinates(radians(90), radians(0));
-  const b = fromGeodeticCoordinates(radians(60), radians(10));
-  const c = fromGeodeticCoordinates(radians(50), radians(-20));
+  const a = fromGeodeticCoordinates(radians(0), radians(90));
+  const b = fromGeodeticCoordinates(radians(10), radians(60));
+  const c = fromGeodeticCoordinates(radians(-20), radians(50));
 
   // Find the mean position, M. Note that the calculation is independent of the
   // heights/depths of the positions.
diff --git a/test/vitest/examples/example-08.spec.ts b/test/vitest/examples/example-08.spec.ts
index 9d38661..ece123f 100644
--- a/test/vitest/examples/example-08.spec.ts
+++ b/test/vitest/examples/example-08.spec.ts
@@ -24,7 +24,7 @@ test("Example 8", () => {
   // PROBLEM:
 
   // Position A is given as n-vector:
-  const a = fromGeodeticCoordinates(radians(80), radians(-90));
+  const a = fromGeodeticCoordinates(radians(-90), radians(80));
 
   // We also have an initial direction of travel given as an azimuth (bearing)
   // relative to north (clockwise), and finally the distance to travel along a
@@ -77,7 +77,7 @@ test("Example 8", () => {
   );
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(b);
+  const [lon, lat] = toGeodeticCoordinates(b);
 
   expect(degrees(lat)).toBeCloseTo(79.99154867339445, 13);
   expect(degrees(lon)).toBeCloseTo(-90.01769837291397, 13);
diff --git a/test/vitest/examples/example-09.spec.ts b/test/vitest/examples/example-09.spec.ts
index f1467cf..0b19bcf 100644
--- a/test/vitest/examples/example-09.spec.ts
+++ b/test/vitest/examples/example-09.spec.ts
@@ -26,12 +26,12 @@ test("Example 9", () => {
   // the two positions are not antipodal).
 
   // Path A is given by a1 and a2:
-  const a1 = fromGeodeticCoordinates(radians(50), radians(180));
-  const a2 = fromGeodeticCoordinates(radians(90), radians(180));
+  const a1 = fromGeodeticCoordinates(radians(180), radians(50));
+  const a2 = fromGeodeticCoordinates(radians(180), radians(90));
 
   // While path B is given by b1 and b2:
-  const b1 = fromGeodeticCoordinates(radians(60), radians(160));
-  const b2 = fromGeodeticCoordinates(radians(80), radians(-140));
+  const b1 = fromGeodeticCoordinates(radians(160), radians(60));
+  const b2 = fromGeodeticCoordinates(radians(-140), radians(80));
 
   // Find the position C where the two paths intersect.
 
@@ -53,7 +53,7 @@ test("Example 9", () => {
   const c = apply((n) => Math.sign(dot(cTmp, a1)) * n, cTmp);
 
   // Use human-friendly outputs:
-  const [lat, lon] = toGeodeticCoordinates(c);
+  const [lon, lat] = toGeodeticCoordinates(c);
 
   expect(degrees(lat)).toBeCloseTo(74.16344802135536, 16);
   expect(degrees(lon)).toBeCloseTo(180, 16);
diff --git a/test/vitest/examples/example-10.spec.ts b/test/vitest/examples/example-10.spec.ts
index 2627fc6..07e78d4 100644
--- a/test/vitest/examples/example-10.spec.ts
+++ b/test/vitest/examples/example-10.spec.ts
@@ -21,10 +21,10 @@ test("Example 10", () => {
   // Path A is given by the two n-vectors a1 and a2 (as in the previous
   // example):
   const a1 = fromGeodeticCoordinates(radians(0), radians(0));
-  const a2 = fromGeodeticCoordinates(radians(10), radians(0));
+  const a2 = fromGeodeticCoordinates(radians(0), radians(10));
 
   // And a position B is given by b:
-  const b = fromGeodeticCoordinates(radians(1), radians(0.1));
+  const b = fromGeodeticCoordinates(radians(0.1), radians(1));
 
   // Find the cross track distance between the path A (i.e. the great circle
   // through a1 and a2) and the position B (i.e. the shortest distance at the
diff --git a/test/vitest/rotation-matrix.spec.ts b/test/vitest/rotation-matrix.spec.ts
index c7acc57..2418df4 100644
--- a/test/vitest/rotation-matrix.spec.ts
+++ b/test/vitest/rotation-matrix.spec.ts
@@ -158,7 +158,7 @@ describe("toRotationMatrixUsingWanderAzimuth()", () => {
           // Avoid situations where components of the xyz2R matrix are close
           // to zero. The Python implementation rounds to zero in these cases,
           // which produces very different results.
-          const [latitude, longitude] = toGeodeticCoordinates(vector, frame);
+          const [longitude, latitude] = toGeodeticCoordinates(vector, frame);
           const rotation = eulerXYZToRotationMatrix(
             longitude,
             -latitude,