I seem to have found an error in the "reckon" function.

I've been comparing solutions to the "direct" geodesic problem (find the ending (lat, lon) and azimuth, given a starting (lat, lon) azimuth and distance). "reckon" provides different ending longitudes between true geodesic using WGS-84 ellipsoid, and a sphere with radius = WGS-84 semi-major axis. Here's the test case:

1) WGS-84 geodesic

>> ellipsoid = referenceEllipsoid('wgs84');

>> start_lat_deg = 90; start_lon_deg = 57; arclen_m = 100e3; az_CWfromN_deg = 30;

>> [lat_reckon, lon_reckon] = reckon(start_lat_deg, start_lon_deg, arclen_m, az_CWfromN_deg, ellipsoid)

lat_reckon =

89.1047

lon_reckon =

-153.0000

2) Great circle

>> earth_radius_km = 6378.137;

>> ellipsoid = referenceEllipsoid;

ellipsoid.LengthUnit = 'm';

sphere_radius_m = 1.e3 * earth_radius_km;

ellipsoid.SemimajorAxis = sphere_radius_m;

ellipsoid.Eccentricity = 0;

>> [lat_reckon, lon_reckon] = reckon(start_lat_deg, start_lon_deg, arclen_m, az_CWfromN_deg, ellipsoid)

lat_reckon =

89.1017

lon_reckon =

57.0000

It's understandable that the two latitudes differ, but the great cicle ending longitude is 210 degrees from the geodesic ending longitude! The "greatcirclefwd" function is used by "reckon" to do the great circle forward problem, and "greatcirclefwd" changes the azimuth in this case from 30 deg to 180 deg (see line 19).

I checked against Charles Karney's "GeodSolve" utility (https://geographiclib.sourceforge.io/cgi-bin/GeodSolve?type=D&input=90+57+30+100e3&format=g&azi2=f&unroll=r&prec=3&radius=6378137&flattening=1%2F10000000&option=Submit), and that gives the final longitude as -153 deg.

What am I missing?