fortyninemaps / karta

A tidy Python package for geospatial computation
https://karta.fortyninemaps.com
MIT License
100 stars 11 forks source link

convergence failure for newton iteration on inverse problem #26

Closed njwilson23 closed 9 years ago

njwilson23 commented 9 years ago

Solving the inverse problem (0.000000, -57.244630) -> (179.795459, -12.169046) on the WGS84 ellipsoid results in a convergence failure warning from the Newton solver and an inaccurate result.

In [3]: LonLatWGS84.inverse(0.000000, -57.244630, 179.795459, -12.169046)
RuntimeWarning: Convergence failure (0.000000, -57.244630) -> (179.795459, -12.169046)
Out [3]: (-168.75633860993707, -173.7906929918106, 12246936.66410405)

libproj result:

In [6]: import pyproj
In [7]: g = pyproj.Geod("+ellps=WGS84")
In [8]: g.inv(0.000000, -57.244630, 179.795459, -12.169046)
Out[8]: (179.786987431945, -179.88183727918067, 12310967.277411604)

geographiclib result:

    ellipsoid (a f)     = 6378137 1/298.257223563 (WGS84)
    status              = OK

    lat1 lon1 fazi1 (°) = -57.24463000 0.00000000 179.78698743
    lat2 lon2 fazi2 (°) = -12.16904600 179.79545900 0.11816272
    s12 (m)             = 12310967.277

    a12 (°)             = 110.71332679
    m12 (m)             = 5987750.843
    M12 M21             = -0.3469428985 -0.3523063207
    S12 (m^2)           = -127281236984754

matlab result:

>> azimuth( -57.244630, 0, -12.169046, 179.795459, referenceEllipsoid('wgs84'))
ans =
  179.7870
>> azimuth(-12.169046, 179.795459, -57.244630, 0, referenceEllipsoid('wgs84'))
ans =
  180.1182
njwilson23 commented 9 years ago

Fixed by 4c3f6aad9345d3c01596704c17acee4d2488f08b