# use of GeoDistance on WGS84 spheroid

The reference WGS84 describes the Earth as oblate spheroid with semi-minor axis (radius to poles) of 6356.752 km and semi-major axis (radius to equator) of 6378.137 km.

I am considering the distance along the 50th parallel (50 degrees north) from 8° E to 11° E. The radius of the circle parallel to the equatorial plane at 50° N, calculated by is x = 4091.704 km, and 3/360 of the circumference is 214.241 km, which is the lenght of the line under consideration.

Calculating the geodetic distance between the two end-points of this line, which is assumed to be shorter than the line along the 50th parallel, by GeoDistance[{50,8},{50,11}] results in 215.073 km, which is longer.

What went wrong?

Let's use the formula for an ellipse and equation for the geocentric latitude, $\psi$, in terms of the geographic latitude, $\phi$, $$\frac{x^2}{r_x^2}+\frac{y^2}{r_y^2} = 1 \hspace{3em} \tan\psi = \frac{x}{y} = \frac{r_y^2}{r_x^2}\tan\phi$$

See the sketch at geocentric latitude. We solve for the $x$-coordinate of a position at latitude 50 degrees and use the value as the radius of our arc.

eqns = {x^2/rx^2 + y^2/ry^2 == 1,
y/x == ry^2/rx^2 Tan[ϕ]
};
soln = NSolve[eqns, {x, y}] // Last;
{x, y} = ({x, y} /. soln);

rx = 6378.1370~Quantity~"Kilometers";
ry = 6356.7523~Quantity~"Kilometers";
ϕ = 50 π/180;
long1 = 8 π/180;
long2 = 11 π/180;
arc = x (long2 - long1)

(*   215.087 km   *)


This distance is slightly greater than the geodetic distance given by

pos1 = GeoPosition[{50, 8}, "ITRF00"];
pos2 = GeoPosition[{50, 11}, "ITRF00"];

GeoDistance[pos1, pos2]~UnitConvert~"Kilometers"

(*   215.073   *)


We can also calculate the distance along the east-west path by

pts = Table[{50, λ}, {λ, 8, 11, .3}];
GeoLength[GeoPath[pts]]~UnitConvert~"Kilometers"

(*   215.087   *)


Another way to calculate the geodetic distance is to consider the great ellipse through the two points and in the plane of the center of the earth. We calculate the radius r from the center to the points, the colatitude θ, unit vectors from the center to the two points, the angle α between the two points in the plane of the great ellipse, and the semi-minor axis of the great ellipse as

r = Sqrt[x^2 + y^2]~UnitConvert~"Kilometers"
θ = ArcTan[y, x];
vec1 = {Cos[long1] Sin[θ], Sin[long1] Sin[θ], Cos[θ]};
vec2 = {Cos[long2] Sin[θ], Sin[long2] Sin[θ], Cos[θ]};
α = ArcCos[vec1.vec2] // N;
ryy = rx r Cos[α/2]/Sqrt[rx^2 - r^2 Sin[α/2]^2]

(*   6365.63 km
6365.63 km   *)


Since the semi-minor axis and the radius to either point is approximately the same, we use a great circle through the points to calculate the gedetic distance as

α ryy

(*   215.073 km   *)