Skip to content

Commit 2a20c3c

Browse files
felixpalmerfelixpalmer
authored
Authalic conversion functions (#21)
1 parent b35bb44 commit 2a20c3c

3 files changed

Lines changed: 258 additions & 0 deletions

File tree

‎authalic_constants.py‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
# Fundamental flattening constant from WGS84 definition
2+
INVERSE_FLATTENING = 298.257223563
3+
4+
# Third flattening is given by:
5+
# n = (a - b) / (a + b), where a is the semi-major axis and b is the semi-minor axis of the ellipsoid
6+
# The first flattening (f) is given by: 1 / f = (a - b) / a, thus rewrite n as:
7+
THIRD_FLATTENING = 1 / (2 * INVERSE_FLATTENING - 1)
8+
9+
# Constants for computing coefficients from: https://arxiv.org/pdf/2212.05818
10+
C_ξφ = [
11+
-4/3.0, -4/45.0, 88/315.0, 538/4725.0, 20824/467775.0, -44732/2837835.0,
12+
34/45.0, 8/105.0, -2482/14175.0, -37192/467775.0, -12467764/212837625.0,
13+
-1532/2835.0, -898/14175.0, 54968/467775.0, 100320856/1915538625.0,
14+
6007/14175.0, 24496/467775.0, -5884124/70945875.0,
15+
-23356/66825.0, -839792/19348875.0,
16+
570284222/1915538625.0
17+
]
18+
19+
C_φξ = [
20+
4/3.0, 4/45.0, -16/35.0, -2582/14175.0, 60136/467775.0, 28112932/212837625.0,
21+
46/45.0, 152/945.0, -11966/14175.0, -21016/51975.0, 251310128/638512875.0,
22+
3044/2835.0, 3802/14175.0, -94388/66825.0, -8797648/10945935.0,
23+
6059/4725.0, 41072/93555.0, -1472637812/638512875.0,
24+
768272/467775.0, -455935736/638512875.0,
25+
4210684958/1915538625.0
26+
]
27+
28+
def compute_coefficients(n, C):
29+
"""Compute series expansion coefficients based on Horner's method.
30+
31+
Args:
32+
n: Third flattening of the ellipsoid
33+
C: Coefficients for generation of the series expansion
34+
"""
35+
cp = [0.0] * 6
36+
37+
# Using Horner's method to compute coefficients
38+
d = n
39+
cp[0] = (((((C[5] * n + C[4]) * n + C[3]) * n + C[2]) * n + C[1]) * n + C[0]) * d
40+
d *= n
41+
cp[1] = ((((C[10] * n + C[9]) * n + C[8]) * n + C[7]) * n + C[6]) * d
42+
d *= n
43+
cp[2] = (((C[14] * n + C[13]) * n + C[12]) * n + C[11]) * d
44+
d *= n
45+
cp[3] = ((C[17] * n + C[16]) * n + C[15]) * d
46+
d *= n
47+
cp[4] = (C[19] * n + C[18]) * d
48+
d *= n
49+
cp[5] = C[20] * d
50+
51+
return cp
52+
53+
def main():
54+
# Compute coefficients
55+
geodetic_to_authalic = [float(f"{x:.16e}") for x in compute_coefficients(THIRD_FLATTENING, C_ξφ)]
56+
authalic_to_geodetic = [float(f"{x:.16e}") for x in compute_coefficients(THIRD_FLATTENING, C_φξ)]
57+
58+
# Print coefficients
59+
print("geodetic_to_authalic:")
60+
for x in geodetic_to_authalic:
61+
print(f" {x:.16e},")
62+
print("\nauthalic_to_geodetic:")
63+
for x in authalic_to_geodetic:
64+
print(f" {x:.16e},")
65+
66+
if __name__ == "__main__":
67+
main()

‎modules/core/authalic.ts‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
import { Radians } from "./coordinate-systems";
2+
3+
// Authalic conversion coefficients obtained from: https://arxiv.org/pdf/2212.05818
4+
// See: authalic_constants.py for the derivation of the coefficients
5+
const GEODETIC_TO_AUTHALIC = new Float64Array([
6+
-2.2392098386786394e-03,
7+
2.1308606513250217e-06,
8+
-2.5592576864212742e-09,
9+
3.3701965267802837e-12,
10+
-4.6675453126112487e-15,
11+
6.6749287038481596e-18
12+
]);
13+
14+
const AUTHALIC_TO_GEODETIC = new Float64Array([
15+
2.2392089963541657e-03,
16+
2.8831978048607556e-06,
17+
5.0862207399726603e-09,
18+
1.0201812377816100e-11,
19+
2.1912872306767718e-14,
20+
4.9284235482523806e-17
21+
]);
22+
23+
// Adaptation of applyCoefficients from DGGAL project: authalic.ec
24+
//
25+
// BSD 3-Clause License
26+
//
27+
// Copyright (c) 2014-2025, Ecere Corporation
28+
//
29+
// Redistribution and use in source and binary forms, with or without
30+
// modification, are permitted provided that the following conditions are met:
31+
//
32+
// 1. Redistributions of source code must retain the above copyright notice, this
33+
// list of conditions and the following disclaimer.
34+
//
35+
// 2. Redistributions in binary form must reproduce the above copyright notice,
36+
// this list of conditions and the following disclaimer in the documentation
37+
// and/or other materials provided with the distribution.
38+
//
39+
// 3. Neither the name of the copyright holder nor the names of its
40+
// contributors may be used to endorse or promote products derived from
41+
// this software without specific prior written permission.
42+
//
43+
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
44+
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
45+
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
46+
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
47+
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
48+
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
49+
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
50+
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
51+
// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
52+
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
53+
54+
/**
55+
* Applies coefficients using Clenshaw summation algorithm (order 6)
56+
* @param phi Angle in radians
57+
* @param C Array of coefficients
58+
* @returns Transformed angle in radians
59+
*/
60+
function applyCoefficients(phi: Radians, C: Float64Array): Radians {
61+
const sinPhi = Math.sin(phi);
62+
const cosPhi = Math.cos(phi);
63+
const X = 2 * (cosPhi - sinPhi) * (cosPhi + sinPhi);
64+
let u0, u1;
65+
66+
u0 = X * C[5] + C[4];
67+
u1 = X * u0 + C[3];
68+
u0 = X * u1 - u0 + C[2];
69+
u1 = X * u0 - u1 + C[1];
70+
u0 = X * u1 - u0 + C[0];
71+
72+
return phi + 2 * sinPhi * cosPhi * u0 as Radians;
73+
}
74+
75+
export function geodeticToAuthalic(phi: Radians): Radians {
76+
return applyCoefficients(phi, GEODETIC_TO_AUTHALIC);
77+
}
78+
79+
export function authalicToGeodetic(phi: Radians): Radians {
80+
return applyCoefficients(phi, AUTHALIC_TO_GEODETIC);
81+
}

‎tests/authalic.test.ts‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
import { describe, it, expect } from 'vitest'
2+
import { geodeticToAuthalic, authalicToGeodetic } from 'a5/core/authalic'
3+
import type { Radians } from 'a5/core/coordinate-systems'
4+
5+
describe('Authalic Conversion', () => {
6+
describe('geodeticToAuthalic', () => {
7+
it('converts zero latitude', () => {
8+
const result = geodeticToAuthalic(0 as Radians)
9+
expect(result).toBeCloseTo(0, 10)
10+
})
11+
12+
it('converts small latitudes', () => {
13+
const input = Math.PI / 180 as Radians // 1 degree
14+
const result = geodeticToAuthalic(input)
15+
expect(result).toBeCloseTo(input, 2)
16+
})
17+
18+
it('converts medium latitudes', () => {
19+
const input = Math.PI / 4 as Radians // 45 degrees
20+
const result = geodeticToAuthalic(input)
21+
expect(result).toBeCloseTo(input, 2)
22+
})
23+
24+
it('converts large latitudes', () => {
25+
const input = Math.PI / 2 as Radians // 90 degrees
26+
const result = geodeticToAuthalic(input)
27+
expect(result).toBeCloseTo(input, 2)
28+
})
29+
30+
it('converts negative latitudes', () => {
31+
const input = -Math.PI / 4 as Radians // -45 degrees
32+
const result = geodeticToAuthalic(input)
33+
expect(result).toBeCloseTo(input, 2)
34+
})
35+
})
36+
37+
describe('authalicToGeodetic', () => {
38+
it('converts zero latitude', () => {
39+
const result = authalicToGeodetic(0 as Radians)
40+
expect(result).toBeCloseTo(0, 10)
41+
})
42+
43+
it('converts small latitudes', () => {
44+
const input = Math.PI / 180 as Radians // 1 degree
45+
const result = authalicToGeodetic(input)
46+
expect(result).toBeCloseTo(input, 2)
47+
})
48+
49+
it('converts medium latitudes', () => {
50+
const input = Math.PI / 4 as Radians // 45 degrees
51+
const result = authalicToGeodetic(input)
52+
expect(result).toBeCloseTo(input, 2)
53+
})
54+
55+
it('converts large latitudes', () => {
56+
const input = Math.PI / 2 as Radians // 90 degrees
57+
const result = authalicToGeodetic(input)
58+
expect(result).toBeCloseTo(input, 2)
59+
})
60+
61+
it('converts negative latitudes', () => {
62+
const input = -Math.PI / 4 as Radians // -45 degrees
63+
const result = authalicToGeodetic(input)
64+
expect(result).toBeCloseTo(input, 2)
65+
})
66+
})
67+
68+
describe('round-trip conversion', () => {
69+
it('preserves latitude through geodetic->authalic->geodetic conversion for all latitudes', () => {
70+
for (let deg = -90; deg <= 90; deg++) {
71+
const lat = (deg * Math.PI / 180) as Radians
72+
const authalic = geodeticToAuthalic(lat)
73+
const geodetic = authalicToGeodetic(authalic)
74+
expect(geodetic).toBeCloseTo(lat, 15)
75+
}
76+
})
77+
78+
it('preserves latitude through authalic->geodetic->authalic conversion for all latitudes', () => {
79+
for (let deg = -90; deg <= 90; deg++) {
80+
const lat = (deg * Math.PI / 180) as Radians
81+
const geodetic = authalicToGeodetic(lat)
82+
const authalic = geodeticToAuthalic(geodetic)
83+
expect(authalic).toBeCloseTo(lat, 15)
84+
}
85+
})
86+
})
87+
88+
describe('specific conversion values', () => {
89+
it('matches reference conversion values', () => {
90+
const testCases = [
91+
{ geodetic: -90, authalic: -90.0000 },
92+
{ geodetic: -67.5, authalic: -67.4092 },
93+
{ geodetic: -45, authalic: -44.8717 },
94+
{ geodetic: -22.5, authalic: -22.4094 },
95+
{ geodetic: 0, authalic: 0 },
96+
{ geodetic: 22.5, authalic: 22.4094 },
97+
{ geodetic: 45, authalic: 44.8717 },
98+
{ geodetic: 67.5, authalic: 67.4092 },
99+
{ geodetic: 90, authalic: 90.0000 }
100+
]
101+
102+
testCases.forEach(({ geodetic, authalic }) => {
103+
const geodeticRad = (geodetic * Math.PI / 180) as Radians
104+
const authalicRad = (authalic * Math.PI / 180) as Radians
105+
const result = geodeticToAuthalic(geodeticRad)
106+
expect(result).toBeCloseTo(authalicRad, 5)
107+
})
108+
})
109+
})
110+
})

0 commit comments

Comments
 (0)