# HP Forums

You're currently viewing a stripped down version of our content. View the full version with proper formatting.
The direct geodesy problem: given Lat Lon distance and bearing find Lat Lon 2.
The inverse: Given two points find the distance.

All the references I have seen, such as Wikipedia, start by giving a spherical solution for the inverse problem, then progress to the Vincenty inverse and direct solutions for ellipsoidal accuracy.

I'm curious - perhaps someone can explain...

Does this mean that the direct solution was not solvable until Vincenty? As far as on-line references go, there is no mention of how it was solved.

Or is it meant to be blindingly obvious that you can just run the spherical code in reverse...surely if it can be done, someone would have implemented it?

The reason I am so interested...
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.
Perhaps because an exact solution doesn't exist ????

Pauli
(12-05-2016 08:29 AM)StephenG1CMZ Wrote: [ -> ]The reason I am so interested ... I was looking for a way of generating some test data for comparison.
Your inquiry is interesting ... however, for the moment I offer a partial solution to your multifaceted posting: generating TEST data - the ONLINE GEODETIC CALCULATOR may facilitate your need.
I have numerous texts and algorithms for avigation, navigation, geodetic computation, Great Circle calculations, etc and may have something close to your desired end product. It may take a little while to research.

after a little research utilizing ... [attachment=4254] as well as ... [attachment=4255] ... for a better understanding ...

that eventually leads to a program at ...
Conversions between Geographical & Transverse Mercator (UTM, GAUSS.KRUGER) Grid Coordinates ... for the HP-67 Calculator.

BEST!
SlideRule

ps: my best reference

[attachment=4259]

from "Solving the Direct & Inverse Geodetic Problems on the Ellipsoid by Numerical Integration"
We have demonstrated how to solve the direct and inverse problems on the ellipsoid by adding the strict solution for the sphere and an ellipsoidal correction determined by numerical integration ... The tests show that both the direct and inverse problem solutions practically agree with Vincenty’s solutions. Although Vincenty’s iterative method apparently is more practical, our method, to our knowledge, is the only independent method for validating Vincenty’s inverse method. (pg 16).
(12-05-2016 04:03 PM)SlideRule Wrote: [ -> ]BEST!
SlideRule

ps: my best reference

<image of cover of "Geodesic Math and How to Use it">

Ahhh! That's Hugh Kenner, PPC member 103, and author of a couple of dozen articles over the first 7 or so years of the club. Also wrote the famous "Calcu-nuts" article for the Wall Street Journal Magazine, a lengthy article about HP calculator fanatics in general and highlighting Bill Wickes and HP-41 synthetic programming. Those were the days.

Jake
(12-05-2016 08:29 AM)StephenG1CMZ Wrote: [ -> ]Does this mean that the direct solution was not solvable until Vincenty? As far as on-line references go, there is no mention of how it was solved.

On a sphere this is quite simple to solve, they probably don't mention it because it's trivial. The main reason it's hard to solve with ellipsoid is because the length of a curve doesn't have an analytic solution for that shape, but it does for spheres (EDIT: it should say "it does for circles" which is what you need to solve the distance on a sphere), just using simple trig. functions.
(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.
I assume 'we' have at least tickled your curiosity w our responses, any alibi's (just curious)? I have one final & hopefully definitive reference from SPON Press

PLANE AND GEODETIC SURVEYING, SECOND EDITION
Aylmer Johnson
CRC Press
© 2014 by Taylor & Francis Group, LLC

8 Geoids, Ellipsoids and Co-ordinate Transforms 101
8.1 Definition of the Geoid 101
8.2 The Need for an Ellipsoid 102
8.3 Orthometric Heights and Bench Marks 106
8.4 Geometry of the Ellipse 108
8.4.1 Defining the Shape of an Ellipse 108
8.4.2 Curvature on an Ellipse 110
8.5 Transformations between Ellipsoids 112
8.5.1 Using a Transform Matrix 112
8.5.2 Making a Transform Matrix 117
8.6 ETRS89 and the International Terrestrial Reference System 118
8.7 Further Properties of Ellipsoids 119
8.7.1 Curvature on an Ellipsoid 119
8.7.1.1 Principal Curvatures 119
8.7.1.2 Curvature in Other Directions 122
8.7.1.3 Average Curvature at a Point 124
8.7.1.4 Mean Curvature of an Ellipsoid 125
8.7.2 Geodesics 125

&
10 Reduction of Distance Measurements 155
10.1 Correction for the Curvature of the Ellipsoid 155
10.2 Correction for Light Curvature 158
10.3 Corrections to Slope Distance Measurements 163
10.4 Final Calculation of Reduced Distance 165
10.5 Slope Distances 167
10.6 Summary 169

&
Appendix C: Worked Example in Transforming between Ellipsoids 211

ALL as apropos

BEST!
SlideRule
My 2 cents: http://proj4.org/geodesic.html seems to have plenty of source code references to delve into.
Yes Sliderule, that has tickled my interest. Thank you all for your helpful replies.
I would love to read those references to help understand and visualise the transformations, rather than just follow them - but my local library rarely has any books.

I also asked the question on Yahoo Answers and Morningfox pointed me to a solution on Movable- types website that I had missed: