How can I calculate the Euclidian displacement of two places on a sphere (earth in this case ) and calculate the

Question

I would like to get the formula on how to calculate the distance between two geographical co-ordinates on earth and heading angle relative to True North. Say from New York to New Dehli , I draw a straight line THROUGH THE EARTH - as they were two points in space. How can I calculate that angle from say New York to New Dehli if I was to draw a straight line through the surface of the earth . What kind of mathematical calculation/formula would be involved in order to do that?

Answer 1

Distance

Convert from spherical coordinates to cartesian ones. If you have latitude $\varphi$ and longitude $\lambda$, you can compute

\begin{align*} x &= r\cos\varphi\cos\lambda \\ y &= r\cos\varphi\sin\lambda \\ z &= r\sin\varphi \end{align*}

This assumes that the earth is a sphere of radius $r$, and has the $z$ axis pointing north and the $x$ axis pointing twoards the intersection of the equator and the prime meridian.

Do this conversion for both points, and you have two positions in $\mathbb R^3$. Compute their element-wise difference and take the length of that, by squaring the elements, adding them up and computing the square root.

$$\ell = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}$$

Direction

In order to compute a kind of bearing from the difference between two locations, you'd most likely want to orthogonally project the displacement vector onto a plane tangent to the surface of the earth at your current location. To do this, let us compute two tangent directions, as derivatives with respect to latitude and longitude. The first points north, the second east.

\begin{align*} n = \frac{\partial(x,y,z)^T}{\partial\varphi} &= \begin{pmatrix} -\sin\varphi\cos\lambda \\ -\sin\varphi\sin\lambda \\ \cos\varphi \end{pmatrix} & e = \frac{\partial(x,y,z)^T}{\partial\lambda} &= \begin{pmatrix} -\cos\varphi\sin\lambda \\ \cos\varphi\cos\lambda \\ 0 \end{pmatrix} \end{align*}

Next you have to normalize these to unit length. $n$ already has unit length, but $e$ has to be rescaled.

\begin{align*} \hat n &= \frac{n}{\lVert n\rVert} = n = \begin{pmatrix} -\sin\varphi\cos\lambda \\ -\sin\varphi\sin\lambda \\ \cos\varphi \end{pmatrix} & \hat e &= \frac{e}{\lVert e\rVert} = \frac1{\cos\varphi}e = \begin{pmatrix} -\sin\lambda \\ \cos\lambda \\ 0 \end{pmatrix} \end{align*}

Then you can use dot products to project your displacement vector onto these directions. This gives you two components, the ratio of which corresponds to the tangens of the azimuth. With zero at north and positive angle going east, you want the north-pointing component in the denominator and the east-pointing component in the numerator.

\begin{align*} \tan\alpha &= \frac{\hat e\cdot(p_2-p_1)}{\hat n\cdot(p_2-p_1)} \\&= \frac {-\sin\lambda_1(x_2-x_1) +\cos\lambda_1(y_2-y_1)} {-\sin\varphi_1\cos\lambda_1(x_2-x_1) -\sin\varphi_1\sin\lambda_1(y_2-y_1) +\cos\varphi_1(z_2-z_1)} \end{align*}

Note that there are two solutions to this, so you either have to check the signs or use an atan2 function.

If you wanted to, you could perform additional computations to determine the elevation, i.e. the angle between the displacement vector and the tangent plane.

Properties of this bearing

I would assume the difference between magnetic north and true north to be rather small in most places, so I'm not sure how much more detail you can expect from a computation like this, with simplifications like a spherical earth and so on.

The initial bearing chosen in this fashion is identical to the one obtained for a greatcircle course. The reason is because in both of these cases you obtain the bearing by intersecting the tangent plane with the plane spanned by origin, destination and center of the earth. From the USA to Mecca that bearing would be mostly NE, according to Google Earth:

I hope this matches your question. If the answer didn't match the answer you expected, then you should probably confer with the authorities which caused that expectation, to see whether you asked the right question in the first place. The rhumb line you mention as an alternative has certain benefits for simple compass navigation, but is not related to an orthogonal projection of the direct displacement vector, or to the shortest distance either through the earth or along its surface.

Answer 2

Rhumb Line Navigation

Rhumb lines or loxodromes are tracks of constant true course. With the exception of meridians and the equator, they are not the same as great circles. They are not very useful approaching either pole, where they become tightly wound spirals. The formulae below fail if any point actually is a pole.

East-West rhumb lines are special. They follow the latitude parallels and form a closed curve. Other rhumb lines extend from pole-to-pole, encircling each pole an infinite number of times. Despite this, they have a finite length given by pi/abs(cos(tc)) (in our angular units, multiply by the radius of the earth to get it in distance units).

When two points (lat1,lon1), (lat2,lon2) are connected by a rhumb line with true course tc :

lon2-lon1=-tan(tc)(log((1+sin(lat2))/cos(lat2))- log((1+sin(lat1))/cos(lat1))) =-tan(tc)(log((1+tan(lat2/2))/(1-tan(lat2/2)))- log((1+tan(lat1/2))/(1-tan(lat1/2)))) =-tan(tc)*(log(tan(lat2/2+pi/4)/tan(lat1/2+pi/4))) (logs are "natural" logarithms to the base e.)

The true course between the points is given by:

tc= mod(atan2(lon1-lon2,log(tan(lat2/2+pi/4)/tan(lat1/2+pi/4))),2*pi) The dist, d between the points is given by:

     if (abs(lat2-lat1) < sqrt(TOL)){
         q=cos(lat1)
     } else {
         q= (lat2-lat1)/log(tan(lat2/2+pi/4)/tan(lat1/2+pi/4))
     }
     d=sqrt((lat2-lat1)^2+ q^2*(lon2-lon1)^2)

This formula fails if the rhumb line in question crosses the 180 E/W meridian. Allowing this as a possibility, the true course tc, and distance d, for the shortest rhumb line connecting two points is given by:

dlon_W=mod(lon2-lon1,2*pi) dlon_E=mod(lon1-lon2,2*pi) dphi=log(tan(lat2/2+pi/4)/tan(lat1/2+pi/4)) if (abs(lat2-lat1) < sqrt(TOL)){ q=cos(lat1) } else { q= (lat2-lat1)/dphi } if (dlon_W < dlon_E){// Westerly rhumb line is the shortest tc=mod(atan2(-dlon_W,dphi),2*pi) d= sqrt(q^2*dlon_W^2 + (lat2-lat1)^2) } else{ tc=mod(atan2(dlon_E,dphi),2*pi) d= sqrt(q^2*dlon_E^2 + (lat2-lat1)^2) }