why is there no spherical solution to the direct geodesy problem? There is!

[Claudio L.]

(12-05-2016 08:29 AM)StephenG1CMZ Wrote:  I have attempted an implementation of the Vincenty direct formula, but it is not correct - so I was looking for a way of generating some test data for comparison.

Here's my simple version of how to solve the direct problem on a sphere:

Given a point P on a sphere R. Point P is given in spherical coordinates by [R longitude latitude]. You want to find coordinate for point Q, traveling a distance S in the bearing B (measured from North, clockwise).

The entire travel will happen in a plane passing through the center of the sphere, that contains P and Q, so it's perfectly defined by those 3 points.
Assume a coordinate system inside that plane, with 3 unit vectors [i j k], where i coincides with the direction of P, j coincides with the direction of bearing and k is just perpendicular to the others.

That plane intersects the sphere in a circle. The arc traveled on that circle will cover an angle S/R (in radians).

Local coordinates of P = [ R 0 0 ]
Local coordinates of Q are trivial: [ R*cos(S/R) , R*sin(S/R) , 0 ]

And the problem is basically solved. We need to express Q in global coordinates, for which we need to determine the actual directions for i,j,k in our global system x,y,z:

The direction of the first vector, i is the same as P:

P0 = [ R*cos(Long.)*cos(Lat.), R*sin(Long.)*cos(Lat.), R*sin(Lat.) ]
then the unit vector is:
i = [ cos(Long.)*cos(Lat.), sin(Long.)*cos(Lat.), sin(Lat.) ]

To determine j is a bit more involved, we need to create a vector perpendicular to i first, in the i-z plane:

k0 = [ -sin(Long.), cos(Long.), 0 ] (this is perpendicular to P in the xy plane)

j0 = k0 ^ i = [ -cos(Long.)*sin(Lat.) , -sin(Long.)*sin(Lat.), cos^2(Long.)*cos(Lat.)+sin^2(Long.)*cos(Lat.) ]

j0 = [ -cos(Long.)*sin(Lat.) , -sin(Long.)*sin(Lat.), cos(Lat.) ]

This is pointing North, so it needs to be rotated about the 'i' axis an angle B to account for the bearing.

Here there's some detail in how to get to it, but basically the rotation matrix on a vector u is given by this.

Just assemble the rotation matrix Rx and do:

j = Rx * j0

Finally, we need k, which we can easily compute as:

k = i ^ j (cross product)

The coordinate transformation matrix A is assembled by putting [ i j k ] in its columns.

Then you have:

Q'=A*Q

And you have Q in global coordinates (Q').

There's one more simple step to determine latitude and longitude from Q'.

EDIT: I should clarify this is only for sufficiently short distances S. For long distances, you need to break it into steps, making P=Q(i-1), computing Qi for each step. To use for comparison with more advanced methods, you could use this computing a new sphere radius Ri at each step (from the definition of the ellipsoid you are trying to compare to), and scaling Qi to that new length Ri (in other words, projecting it onto the ellipsoid at each step), it should give you relatively ok approximations.