refactor: use non-iterative Bowring algorithm for ecef2geodetic#1126
Conversation
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.
|
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. |
|
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). |
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
|
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: |
|
So errors/differences are below 3mm until Pluto? If so, please merge! Great results! |
|
That's at least what my test is telling. And of course only in comparison to our old method, no idea about absolute error. ;-) |
|
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. |
|
@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): |
|
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. |
|
The algorithm doesn't require any special treatment for the poles. Latitudes ±89.99 and ±90 are included in the "Near-pole" test category. |
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
|
I've added statistics comparing za/aa using |
|
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. |
|
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. |
|
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. |
|
As discussed, aa treatment will be handled in a separate PR. |
Replace the iterative ECEF-to-geodetic coordinate conversion with Bowring's closed-form algorithm.
Key changes:
This simplifies the implementation while maintaining accuracy across all latitudes and ellipsoid shapes. The new implementation is about 5x faster on average.
test_path_pointmakes heavy use ofecef2geodeticand 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: