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

[noritada]

When azimuth is equal or close to 90 or 270, loxodrome_direct will return an incorrect longitude value.

The following code will reproduce the problem:

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:

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, -19494.22711642675131
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, -6404.74237214225013
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

At the very least, the values of longitude when the azimuth approaching 90 or 270 to less than 1e-5 look obviously wrong. On the other hand, the latitude values look fine, including when the azimuth is close to 0 or 180.

pymap3d version: 2.7.0

Many thanks.

[noritada]

tan(a12) changes significantly with respect to a12 around a12=pi/2, and is more susceptible to rounding errors. I think the equation below needs to be replaced with some other code around a12=pi/2.

    newiso = geodetic2isometric(lat2, ell, deg=False)
    iso = geodetic2isometric(lat1, ell, deg=False)

    dlon = tan(a12) * (newiso - iso)