Skip to content

refactor: use non-iterative Bowring algorithm for ecef2geodetic#1126

Merged
olemke merged 6 commits into
atmtools:mainfrom
olemke:ecef2geodetic-non-iterative
May 22, 2026
Merged

refactor: use non-iterative Bowring algorithm for ecef2geodetic#1126
olemke merged 6 commits into
atmtools:mainfrom
olemke:ecef2geodetic-non-iterative

Conversation

@olemke
Copy link
Copy Markdown
Member

@olemke olemke commented May 12, 2026

Replace the iterative ECEF-to-geodetic coordinate conversion with Bowring's closed-form algorithm.

Key changes:

  • Removes the iterative loop and its convergence check, eliminating the max_ecef2geodetic_iter limit and associated error handling.
  • Removes special-casing for spherical ellipsoids and near-poles; the closed-form approach handles these robustly.
  • Uses a robust height formula that avoids division by cos(lat).

This simplifies the implementation while maintaining accuracy across all latitudes and ellipsoid shapes. The new implementation is about 5x faster on average.

test_path_point makes heavy use of ecef2geodetic and compares results against calculated reference values. This test continues to pass with the new algorithm.

Accuracy and performance comparison between the iterative (old) and the non-iterative (new) implementation:

=== Accuracy comparison (old vs new) ===
WGS-84 overall:
  max_dh   = 5.587935e-09 m
  max_dlat = 7.663132e-10 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 7.036659e-10 m
  avg_dlat = 4.240752e-11 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 810
Near-pole:
  max_dh   = 1.862645e-09 m
  max_dlat = 0.000000e+00 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 1.086543e-09 m
  avg_dlat = 0.000000e+00 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 216
Mid-latitude:
  max_dh   = 1.862645e-09 m
  max_dlat = 7.663132e-10 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.069606e-10 m
  avg_dlat = 6.953180e-11 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 486
Sphere:
  max_dh   = 9.313226e-10 m
  max_dlat = 1.908518e-11 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.161588e-10 m
  avg_dlat = 2.548511e-12 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 810
Runtime:
  ecef2geodetic_old: 145929 μs
  ecef2geodetic_new: 28419 μs

Replace the iterative ECEF-to-geodetic conversion with Bowring's
closed-form algorithm. This removes special-casing for spherical
ellipsoids and near-poles, and eliminates the convergence loop
and its associated iteration limit.
Comment thread src/core/geodesy/geodetic.cpp Outdated
@riclarsson
Copy link
Copy Markdown
Contributor

I searched for the algorithm. Can you check the difference between these two algorithms far from Earth. e.g., the moon, Venus, Jupiter, Pluto. My google result told me it starts losing accuracy perhaps as distances reach 10,000 km.

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 13, 2026

That's a very good point. I'll do comparisons for larger distances. Give me some time, I'll be off for a couple of days (back next Wednesday).

olemke added 3 commits May 20, 2026 11:28
Replace direct computation of e2 and ep2 with algebraically
equivalent expressions that avoid catastrophic cancellation.

Using (a-b)(a+b) = a² - b² identity prevents precision loss
when a and b are close in value (as is the case for Earth
reference ellipsoids where flattening is small).

- Compute a_minus_b = a - b and a_plus_b = a + b
- Express e2 and ep2 using these intermediates
@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 20, 2026

The comparison to the old method has been added as a test case. Accuracy seems to be fine (at least in comparison to our old method). Here is the test output. It shows separate stats for low (<=100km) and high (>100km) altitudes. High altitudes include distances to several objects from the Moon to Pluto:

=== Accuracy comparison (old vs new) ===
WGS-84 overall (low altitude <=100km):
  max_dh   = 6.519258e-09 m
  max_dlat = 7.663061e-10 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 6.530756e-10 m
  avg_dlat = 4.240701e-11 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 810
WGS-84 overall (high altitude >100km):
  max_dh   = 2.075195e-03 m
  max_dlat = 5.463190e-08 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 3.243281e-04 m
  avg_dlat = 2.779372e-09 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 945
Near-pole (low altitude <=100km):
  max_dh   = 2.793968e-09 m
  max_dlat = 0.000000e+00 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 6.208817e-10 m
  avg_dlat = 0.000000e+00 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 216
Near-pole (high altitude >100km):
  max_dh   = 1.953125e-03 m
  max_dlat = 1.421085e-14 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.236694e-04 m
  avg_dlat = 1.015061e-15 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 252
Mid-latitude (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 7.663061e-10 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.951104e-10 m
  avg_dlat = 6.953095e-11 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 486
Mid-latitude (high altitude >100km):
  max_dh   = 1.953125e-03 m
  max_dlat = 5.463190e-08 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.745006e-04 m
  avg_dlat = 4.557258e-09 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 567
Sphere (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 1.908518e-11 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 2.161588e-10 m
  avg_dlat = 2.548511e-12 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 810
Sphere (high altitude >100km):
  max_dh   = 9.765625e-04 m
  max_dlat = 1.908518e-11 deg
  max_dlon = 0.000000e+00 deg
  avg_dh   = 7.737809e-05 m
  avg_dlat = 2.351449e-12 deg
  avg_dlon = 0.000000e+00 deg
  samples  = 945
ecef2geodetic (old): 105717 μs
ecef2geodetic_new: 28240 μs

@riclarsson
Copy link
Copy Markdown
Contributor

So errors/differences are below 3mm until Pluto? If so, please merge! Great results!

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 20, 2026

That's at least what my test is telling. And of course only in comparison to our old method, no idea about absolute error. ;-)

@erikssonpatrick
Copy link
Copy Markdown
Contributor

To get a grip on the absolute accuracy, start with a geodetic position, convert to ECEF and back, and compare the result to the original value. This doesn't tell you how much each step contributes to the error, but presumably the combined error is small enough to be ignored.

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 21, 2026

@erikssonpatrick Good call, should've thought of that myself since I already start with geodetic positions in the test anyway. Here are the round-trip comparison values. The deviations seem to be acceptable (in the same order of magnitude or smaller as for the comparison between old and new):

=== Round-trip error (original geodetic vs ecef2geodetic) ===
WGS-84 round-trip (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 7.662990e-10 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 4.599124e-10 m
  avg_dlat = 4.240510e-11 deg
  avg_dlon = 7.719477e-16 deg
  samples  = 810
WGS-84 round-trip (high altitude >100km):
  max_dh   = 1.953125e-03 m
  max_dlat = 5.463189e-08 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 1.149681e-04 m
  avg_dlat = 2.779374e-09 deg
  avg_dlon = 5.413659e-16 deg
  samples  = 945
Near-pole round-trip (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 0.000000e+00 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 6.208817e-10 m
  avg_dlat = 0.000000e+00 deg
  avg_dlon = 3.947460e-16 deg
  samples  = 216
Near-pole round-trip (high altitude >100km):
  max_dh   = 1.953125e-03 m
  max_dlat = 0.000000e+00 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 1.787459e-04 m
  avg_dlat = 0.000000e+00 deg
  avg_dlon = 6.767074e-16 deg
  samples  = 252
Mid-latitude round-trip (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 7.662990e-10 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 3.181061e-10 m
  avg_dlat = 6.952830e-11 deg
  avg_dlon = 1.052656e-15 deg
  samples  = 486
Mid-latitude round-trip (high altitude >100km):
  max_dh   = 4.882812e-04 m
  max_dlat = 5.463189e-08 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 6.566588e-05 m
  avg_dlat = 4.557260e-09 deg
  avg_dlon = 5.012647e-16 deg
  samples  = 567
Sphere round-trip (low altitude <=100km):
  max_dh   = 9.313226e-10 m
  max_dlat = 7.105427e-15 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 1.678680e-10 m
  avg_dlat = 2.228122e-15 deg
  avg_dlon = 5.789607e-16 deg
  samples  = 810
Sphere round-trip (high altitude >100km):
  max_dh   = 9.765625e-04 m
  max_dlat = 7.105427e-15 deg
  max_dlon = 7.105427e-15 deg
  avg_dh   = 8.099502e-05 m
  avg_dlat = 2.172983e-15 deg
  avg_dlon = 6.917453e-16 deg
  samples  = 945

@erikssonpatrick
Copy link
Copy Markdown
Contributor

We can live with max 2 mm error in altitude! I have done tests like this, way back. And according to my memory, those results were not better. So this looks perfectly fine.
But when at it, why not include za and aa in the testing? To be complete among the five variables specifying position and direction.
And there are dedicated tests for the poles? For exactly lat=90, and values where close to it? Or are these handled by special expressions (like in my old code)?

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 21, 2026

The algorithm doesn't require any special treatment for the poles. Latitudes ±89.99 and ±90 are included in the "Near-pole" test category.
I'll have a look at including aa and za into the test.

Add comprehensive LOS tests to test_geodetic.cpp:
- Compare ecef2geodetic_los internally using old/new ecef2geodetic
  implementations
- Measure round-trip errors for geodetic_los2ecef -> ecef2geodetic_los
@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 21, 2026

I've added statistics comparing za/aa using ecef2geodetic_los using the old and new implementation of ecef2geodetic internally. Round-trip compares the result of geodetic_los2ecef -> ecef2geodetic_los to the original input los. Zenith cases have been separated out because the azimuth angles naturally produce large differences under this condition.

=== LOS accuracy comparison (ecef2geodetic_los_old vs ecef2geodetic_los) ===
LOS (low altitude <=100km):
  max_dza  = 8.018224e-10 deg
  max_daa  = 4.390622e-06 deg
  avg_dza  = 2.465304e-11 deg
  avg_daa  = 2.114028e-08 deg
  samples  = 90882
LOS (high altitude >100km):
  max_dza  = 5.466956e-08 deg
  max_daa  = 3.130177e-04 deg
  avg_dza  = 1.588562e-09 deg
  avg_daa  = 1.385428e-06 deg
  samples  = 106029
Near-pole LOS (low altitude <=100km):
  max_dza  = 3.552714e-14 deg
  max_daa  = 1.457892e-10 deg
  avg_dza  = 3.515499e-15 deg
  avg_daa  = 4.791701e-12 deg
  samples  = 32076
Near-pole LOS (high altitude >100km):
  max_dza  = 5.151435e-14 deg
  max_daa  = 1.457892e-10 deg
  avg_dza  = 7.934318e-16 deg
  avg_daa  = 1.026818e-12 deg
  samples  = 37422
Mid-latitude LOS (low altitude <=100km):
  max_dza  = 8.018224e-10 deg
  max_daa  = 4.390622e-06 deg
  avg_dza  = 4.581863e-11 deg
  avg_daa  = 3.928004e-08 deg
  samples  = 48114
Mid-latitude LOS (high altitude >100km):
  max_dza  = 5.466956e-08 deg
  max_daa  = 3.130177e-04 deg
  avg_dza  = 2.951996e-09 deg
  avg_daa  = 2.574533e-06 deg
  samples  = 56133
LOS at zenith/nadir (low altitude <=100km):
  max_dza  = 0.000000e+00 deg
  max_daa  = 1.800000e+02 deg
  avg_dza  = 0.000000e+00 deg
  avg_daa  = 5.998796e+01 deg
  samples  = 16524
LOS at zenith/nadir (high altitude >100km):
  max_dza  = 0.000000e+00 deg
  max_daa  = 1.800000e+02 deg
  avg_dza  = 0.000000e+00 deg
  avg_daa  = 5.119203e+01 deg
  samples  = 19278
=== LOS round-trip error (geodetic_los2ecef -> ecef2geodetic_los) ===
LOS round-trip (low altitude <=100km):
  max_dza  = 8.191691e-10 deg
  max_daa  = 4.390624e-06 deg
  avg_dza  = 2.820474e-11 deg
  avg_daa  = 1.056923e-08 deg
  samples  = 90882
LOS round-trip (high altitude >100km):
  max_dza  = 5.468720e-08 deg
  max_daa  = 3.130177e-04 deg
  avg_dza  = 1.591584e-09 deg
  avg_daa  = 6.927150e-07 deg
  samples  = 106029
Near-pole LOS round-trip (low altitude <=100km):
  max_dza  = 5.379320e-11 deg
  max_daa  = 1.421085e-13 deg
  avg_dza  = 5.034575e-12 deg
  avg_daa  = 1.407435e-15 deg
  samples  = 32076
Near-pole LOS round-trip (high altitude >100km):
  max_dza  = 5.379320e-11 deg
  max_daa  = 1.421085e-13 deg
  avg_dza  = 5.034540e-12 deg
  avg_daa  = 1.468029e-15 deg
  samples  = 37422
Mid-latitude LOS round-trip (low altitude <=100km):
  max_dza  = 8.191691e-10 deg
  max_daa  = 4.390624e-06 deg
  avg_dza  = 4.856053e-11 deg
  avg_daa  = 1.963998e-08 deg
  samples  = 48114
Mid-latitude LOS round-trip (high altitude >100km):
  max_dza  = 5.468720e-08 deg
  max_daa  = 3.130177e-04 deg
  avg_dza  = 2.953801e-09 deg
  avg_daa  = 1.287268e-06 deg
  samples  = 56133
LOS at zenith/nadir round-trip (low altitude <=100km):
  max_dza  = 0.000000e+00 deg
  max_daa  = 1.800000e+02 deg
  avg_dza  = 0.000000e+00 deg
  avg_daa  = 4.236673e+01 deg
  samples  = 16524
LOS at zenith/nadir round-trip (high altitude >100km):
  max_dza  = 0.000000e+00 deg
  max_daa  = 1.800000e+02 deg
  avg_dza  = 0.000000e+00 deg
  avg_daa  = 4.694658e+01 deg
  samples  = 19278

@riclarsson
Copy link
Copy Markdown
Contributor

Nadir matters. We must define it to be consistent. For polarized calculations, it must be consistently defined at every level, at the very least. If that's different from the old method, fine - apply a rotation. But I'm confused by the comment.

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 21, 2026

All I can currently say is that for the round-trip test, the resulting azimuth deviates on average from the input azimuth by ~45 degrees in the new algorithm and by ~39 degrees with the old algorithm. I'll add specific tests that only change the altitude and keep everything else constant. If I understand you correctly, the returned azimuth should also stay constant for those cases. I can see that this consistency is important.

@erikssonpatrick
Copy link
Copy Markdown
Contributor

Richard is right. The azimuth angle must be maintained for nadir/zenith, to keep track of polarisation. That said, we need to come up with a definition of the polarisation directions for nadir/zenith that makes sense.
Maybe we can take this tomorrow? (It seems that I can join from start anyhow)

@olemke
Copy link
Copy Markdown
Member Author

olemke commented May 22, 2026

As discussed, aa treatment will be handled in a separate PR.

@olemke olemke merged commit becca91 into atmtools:main May 22, 2026
17 of 18 checks passed
@olemke olemke deleted the ecef2geodetic-non-iterative branch May 22, 2026 08:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants