geospace-code / pymap3d

pure-Python (Numpy optional) 3D coordinate conversions for geospace ecef enu eci
https://geospace-code.github.io/pymap3d
BSD 2-Clause "Simplified" License
394 stars 85 forks source link

loxodrome_direct still returns incorrect longitude when azimuth is close to 90 or 270 #43

Closed noritada closed 3 years ago

noritada commented 3 years ago

Thank you for fixing #42, the lox stability issue, and releasing the fix as v2.7.1. However, loxodrome_direct still returns incorrect longitude values, although is is improved before fixing.

The code is the same as one used in #42.

from pymap3d.lox import loxodrome_direct

clat, clon = 35.0, 140.0
az = [
    89,
    89.9,
    89.99,
    89.999,
    89.9999,
    89.99999,
    89.999999,
    90,
    90.000001,
    90.00001,
    90.0001,
    90.001,
    90.01,
    90.1,
    91,
]
tmp = [az_ + 180 for az_ in az]
az += tmp

for az_ in az:
    lat, lon = loxodrome_direct(clat, clon, 50_000, az_)
    print(f"{az_:.6f}, {lat:.14f}, {lon:.14f}")

This code prints lines like this in PyMap 2.7.1:

89.000000, 35.00786565022930, 140.54765888297607
89.900000, 35.00078660501723, 140.54771788308574
89.990000, 35.00007866054496, 140.54771634402164
89.999000, 35.00000786605369, 140.54771605442008
89.999900, 35.00000078660448, 140.54771541554422
89.999990, 35.00000007865957, 140.54770927042696
89.999999, 35.00000000786506, 140.54764724167964
90.000000, 34.99999999999901, 139.55133723928088
90.000001, 34.99999999213296, 140.54778537068364
90.000010, 34.99999992133847, 140.54772297468475
90.000100, 34.99999921339354, 140.54771678232532
90.001000, 34.99999213394431, 140.54771614079507
90.010000, 34.99992133945204, 140.54771583369782
90.100000, 34.99921339487867, 140.54771264295587
91.000000, 34.99213433955717, 140.54760647858132
269.000000, 34.99213433955717, 139.45239352141868
269.900000, 34.99921339487867, 139.45228735704418
269.990000, 34.99992133945204, 139.45228416630223
269.999000, 34.99999213394431, 139.45228385919867
269.999900, 34.99999921339354, 139.45228321768184
269.999990, 34.99999992133847, 139.45227702538708
269.999999, 34.99999999213296, 139.45221462306554
270.000000, 34.99999999999901, 139.55133723928088
270.000001, 35.00000000786506, 139.45235274366772
270.000010, 35.00000007865956, 139.45229076455408
270.000100, 35.00000078660448, 139.45228458430924
270.001000, 35.00000786605369, 139.45228394556526
270.010000, 35.00007866054496, 139.45228365597760
270.100000, 35.00078660501723, 139.45228211691446
271.000000, 35.00786565022930, 139.45234111702391

Now loxodrome_direct returns the same value 139.55133723928088 for both 90 and 270, and both are incorrect.

pymap3d version: 2.7.1

Many thanks.

scivision commented 3 years ago

I see my comparison with matlab used too short a range. Now I see the issue more clearly.

noritada commented 3 years ago

Thank you very much for fixing! It works for me.