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

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?

• Calculating the vector from New York to New Delhi is described by MvG. What do you want to measure the angle from? – Ross Millikan Mar 12 '13 at 14:24
• Given any two points on the Earth, $A$ and $B,$ the straight line to $A$ from $B$, going straight through the Earth as if it were not an obstacle, intersects the great circle through $A$ and $B$ at two points, $A$ and $B$. In order to look in that exact direction, you may first look in the direction of the great-circle route from $B$ to $A,$ and then direct your gaze downward until it is directed toward $A.$ If $B$ is at exactly the opposite side of the Earth from $A$, the direction to $A$ is nearly straight down. (I write "nearly" because the Earth is not a perfect sphere.) – David K Nov 29 '17 at 15:58

# 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.

• If you don't like the spherical assumption, this approach still works if you use a more refined shape of the earth. – Ross Millikan Mar 12 '13 at 14:23
• I would like to know this direction: "Face the QIBLAH (the direction of Mecca): In the U.S. A. the correct direction is generally SOUTHEAST. DO NOT follow the "magnetic" direction suggested by magnetic compasses. They are wrong. Please call us for an explanation. Briefly, the soul cannot be stopped by walls or mountains. Thus, the shortest way to Mecca is right through the earth, and not around the North Pole. During our Salat we bow and prostrate and we face a semi-circle, a wide range of angles." – user2117566 Mar 12 '13 at 19:29
• currently I am using Rhumb line heading , but I am not sure if its giving me the correct result although it still gives me ESE. – user2117566 Mar 12 '13 at 19:33
• The angle should be relative to true north on a x,y plane through the Earth – user2117566 Mar 12 '13 at 20:15
• 'Mecca' is a city in the Middle East. 'Salat' is a form of Contact Prayer with the Creator of the Universe made in that direction as a focal point from any place on Earth – user2117566 Mar 12 '13 at 20:19

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) }