The n-vector page

nvector klode 2.png

Position calculations – simple and exact solutions by means of n-vector

 

n-vector is used to replace latitude and longitude, making the calculations simple and non-singular. Full accuracy is achieved for any global position (and for any distance).

 

This web page contains solutions to ten common geographical position calculations.

Background

Written by Kenneth Gade (PhD, Principal Scientist in the Navigation Group at FFI, Norwegian Defence Research Establishment).

For references, please cite: 

Gade, K. (2010). A Non-singular Horizontal Position Representation, The Journal of Navigation, Volume 63, Issue 03, pp 395-417, July 2010.

Used worldwide for critical applications:

n-vector is used worldwide, e.g. for navigation, traffic control, and surveillance, of aircraft, cars, ships and submarines/submersibles. One example is Thales, who is world-leading at air traffic management (ATM). Their newest generation of ATM system is based on n-vector (replacing previous alternatives).

Solving ten common problems

We will now illustrate how ten common geographical position calculations can be solved using n-vector. The solutions are also available as Matlab code (and as C++, C#, Python, JavaScript, and other languages).

Color coding for figures and variables:

  • Input (red): The initial input that is given
  • Output (green): The solution to be found
  • Function (blue): A function that is available from the downloadable code.

Example 1: A and B to delta

Given two positions A and B. Find the exact vector from A to B in meters north, east and down, and find the direction (azimuth/bearing) to B, relative to north. Use WGS-84 ellipsoid.

Problem

See vector symbols explained for more details about the mathematical notation.
 
Given two positions, A and B as latitudes, longitudes and depths (relative to Earth, E):
– 
la t EA
lon g EA
 and 
z EA
 (depth = –height)
la t EB
lon g EB
 and 
z EB
 
Find the exact vector between the two positions, given in meters north, east, and down, and find the direction (azimuth) to B, relative to north.
 
Details:
 
  • Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface.
  • Use position A to define north, east, and down directions. (Due to the curvature of Earth and different directions to the North Pole, the north, east, and down directions will change (relative to Earth) for different places. Position A must be outside the poles for the north and east directions to be defined.

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
1. First, the given latitudes and longitudes are converted to n-vectors:

n EA E =lat_long2n_E(la t EA ,lon g EA ) n EB E =lat_long2n_E(la t EB ,lon g EB ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGceaqabeaacaWHUb Waa0baaSqaaiaadweacaWGbbaabaGaamyraaaakiabg2da9GWabKaa abbbaaaaaWRaefMCRjeB1vgapeGaa8hBaiaa=fgacaWF0bGaa83xai aa=XgacaWFVbGaa8NBaiaa=DgacaWFYaGaa8NBaiaa=9facaWFfbGc paGaa8hkaabbOpaaaaaasvgza8GacaWGSbGaamyyaiaadshadaWgaa WcbaGaamyraiaadgeaaeqaaOWdaiaa=XcapiGaamiBaiaad+gacaWG UbGaam4zamaaBaaaleaacaWGfbGaamyqaaqabaGcpaGaa8xkaaqaai aah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOGaeyypa0tc aa0dbiaa=XgacaWFHbGaa8hDaiaa=9facaWFSbGaa83Baiaa=5gaca WFNbGaa8Nmaiaa=5gacaWFFbGaa8xraOWdaiaa=HcapiGaamiBaiaa dggacaWG0bWaaSbaaSqaaiaadweacaWGcbaabeaak8aacaWFSaWdci aadYgacaWGVbGaamOBaiaadEgadaWgaaWcbaGaamyraiaadkeaaeqa [email protected]@

 

2. When the positions are given as n-vectors (and depths), it is easy to find the delta vector decomposed in E (using the first function mentioned at the end of n-vector explained):

p AB E =n_EA_E_and_n_EB_E2p_AB_E( n EA E , n EB E , z EA , z EB ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCiCamaaDa aaleaacaWGbbGaamOqaaqaaiaadweaaaGccqGH9aqpimqajaaqqqaa aaaa8karHj3AcXwDLbWdbiaa=5gacaWFFbGaa8xraiaa=feacaWFFb Gaa8xraiaa=9facaWFHbGaa8NBaiaa=rgacaWFFbGaa8NBaiaa=9fa caWFfbGaa8Nqaiaa=9facaWFfbGaa8Nmaiaa=bhacaWFFbGaa8xqai aa=jeacaWFFbGaa8xraOWdaiaa=HcacaWHUbWaa0baaSqaaiaadwea caWGbbaabaGaamyraaaakiaa=XcacaaMe8UaaCOBamaaDaaaleaaca WGfbGaamOqaaqaaiaadweaaaGccaWFSaGaaGjbVdbbOpaaaaaasvgz a8GacaWG6bWaaSbaaSqaaiaadweacaWGbbaabeaak8aacaWFSaGaaG jbV=GacaWG6bWaaSbaaSqaaiaadweacaWGcbaabeaak8aacaWFPaaa [email protected]@

No ellipsoid is specified when calling the function, thus WGS-84 (default) is used.

 

3. We now have the delta vector from A to B, but (as seen from the superscript) the three coordinates of the vector are along the Earth coordinate frame E, while we need the coordinates to be north, east and down. To get this, we define a North-East-Down coordinate frame called N, and then we need the rotation matrix (direction cosine matrix) 
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 to go between E and N (as explained here). In our math library, we have a simple function that calculates 
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 from n-vector, and we use this function (using n-vector at A-position):

R EN =n_E2R_EN( n EA E ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa aaleaacaWGfbGaamOtaaqaaaaakiabg2da9GWabKaaabbbaaaaaWRa efMCRjeB1vgapeGaa8NBaiaa=9facaWFfbGaa8Nmaiaa=jfacaWFFb Gaa8xraiaa=5eak8aacaWFOaGaaCOBamaaDaaaleaacaWGfbGaamyq [email protected]@

 

4. Now the delta vector is easily decomposed in N. When deciding if we should premultiply 
p AB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCiCamaaDa [email protected]@
 with 
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
or 
R NE [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
, we use the "closest-rule" saying that the subscript that is closest to the vector should be the frame where the vector is decomposed (also explained here). Since the vector is decomposed in E, we must use 
R NE [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 (
R NE [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 is the transpose of
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
):

p AB N = R NE p AB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHWbWaa0baaSqaaiaadgeacaWGcbaabaGaamOtaaaak8aacqGH 9aqpcaWHsbWaa0baaSqaaiaad6eacaWGfbaabaaaaOGaaCiCamaaDa [email protected]@

 

5. The three components of 
p AB N [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
 are the north, east and down displacements from A to B in meters. The azimuth is simply found from element 1 and 2 of the vector (the north and east components):

azimuth=atan2 p AB N (2), p AB N (1) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWGHbGaamOEaiaadMgacaWGTbGaamyDaiaadshacaWGObWdaiab g2da9iaabggacaqG0bGaaeyyaiaab6gacaqGYaWaaeWaaeaacaWHWb Waa0baaSqaaiaadgeacaWGcbaabaGaamOtaaaakiaacIcacaaIYaGa aiykaiaacYcacaWHWbWaa0baaSqaaiaadgeacaWGcbaabaGaamOtaa [email protected]@

Example 2: B and delta to C

Given the position of vehicle B and a bearing and distance to an object C. Find the exact position of C. Use WGS-72 ellipsoid.

Problem

See vector symbols explained for more details about the mathematical notation.
 
A radar or sonar attached to a vehicle B (Body coordinate frame) measures the distance and direction to an object C. We assume that the distance and two angles measured by the sensor (typically bearing and elevation relative to B) are already converted (by converting from spherical to Cartesian coordinates) to the vector 
p BC B [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaahchadaqhaaWcbaGaamOqaiaadoeaaeaacaWGcbaaaaaa @[email protected]
 (i.e. the vector from B to C, decomposed in B). The position of B is given as 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaaaa @[email protected]
 and 
z EB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai [email protected]@
, and the orientation (attitude) of B is given as 
R NB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai [email protected]@
 (this rotation matrix can be found from roll/pitch/yaw by using zyx2R.m).
 
Find the exact position of object C as n-vector and depth (
n EC E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
 and 
z EC [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
), assuming Earth ellipsoid with semi-major axis a and flattening f. For WGS-72, use a = 6 378 135 m and f = 1/298.26.

 

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
1. The delta vector is given in B. It should be decomposed in E before using it, and thus we need 
R EB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
. This matrix is found from the matrices 
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 and 
R NB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai [email protected]@
, and we need to find 
R EN [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
, as in Example 1:

R EN =n_E2R_EN( n EB E ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa aaleaacaWGfbGaamOtaaqaaaaakiabg2da9GWabKaaabbbaaaaaWRa efMCRjeB1vgapeGaa8NBaiaa=9facaWFfbGaa8Nmaiaa=jfacaWFFb Gaa8xraiaa=5eak8aacaWFOaaeeG+aaaaaaivzKbWdciaah6gadaqh [email protected]@

 

2. Now, we can find 
R EB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa [email protected]@
 by using that the closest frames cancel when multiplying two rotation matrices (i.e. N is cancelled here):

R EB = R EN R NB [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOuamaaDa aaleaacaWGfbGaamOqaaqaaaaakiabg2da9iaahkfadaqhaaWcbaGa amyraiaad6eaaeaaaaGcqqa6daaaaaGuLrgapeGaaCOuamaaDaaale [email protected]@

 

3. The delta vector is now decomposed in E (details about decomposing are found here):

p BC E = R EB p BC B [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCiCamaaDa aaleaacaWGcbGaam4qaaqaaiaadweaaaGccqGH9aqpcaWHsbWaa0ba aSqaaiaadweacaWGcbaabaaaaOaeeG+aaaaaaivzKbWdbiaahchada [email protected]@

 

4. It is now easy to find the position of C using the second function mentioned here (with custom ellipsoid overriding the default WGS-84).

[ n EC E , z EC ]=n_EA_E_and_p_AB_E2n_EB_E( n EB E , p BC E , z EB ,a,f) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaai4waabaaa q9baaapeGaaCOBamaaDaaaleaacaWGfbGaam4qaaqaaiaadweaaaGc paGaaiila8qacaWG6bWaaSbaaSqaaiaadweacaWGdbaabeaak8aaca GGDbGaeyypa0dcdeqcaaeeeaaaaaaVcquyYTMqSvxza8GacaWFUbGa a83xaiaa=veacaWFbbGaa83xaiaa=veacaWFFbGaa8xyaiaa=5gaca WFKbGaa83xaiaa=bhacaWFFbGaa8xqaiaa=jeacaWFFbGaa8xraiaa =jdacaWFUbGaa83xaiaa=veacaWFcbGaa83xaiaa=veak8aacaWFOa aeeG+aaaaaaivzKbWddiaah6gadaqhaaWcbaGaamyraiaadkeaaeaa caWGfbaaaOaeeaaaaaa6dieB1vgapqGaaGjcVlaayIW7paGaa8hlai aaysW7caWHWbWaa0baaSqaaiaadkeacaWGdbaabaGaamyraaaakiaa yIW7caaMi8Uaa8hlaiaaysW7pmGaamOEamaaBaaaleaacaWGfbGaam OqaaqabaGcpqGaaGPaV=aacaWFSaWddiaadggapaGaa8hlaiaaykW7 [email protected]@

 

Details: For simplicity, the orientation of the vehicle in this example was given relative to N. If the navigation system of the vehicle was using the non-singular coordinate frame L instead, the solution would be the same, but replacing the function n_E2R_EN.m with n_E_and_wa2R_EL.m. See Table 2 in Gade (2010) and the help text in n_E_and_wa2R_EL.m for details.

Example 3: ECEF-vector to geodetic latitude

Given an ECEF-vector of a position. Find geodetic latitude, longitude and height (using WGS-84 ellipsoid).

Problem

See vector symbols explained for more details about the mathematical notation.
 
Position B is given as an “ECEF-vector” 
p EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaahchadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaaaa @[email protected]
 (i.e. a vector from E, the center of the Earth, to B, decomposed in E).
 
Find the geodetic latitude, longitude and height, 
la t EB
, 
lon g EB
 and 
h EB
, assuming WGS-84 ellipsoid.

 

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
1. We have a function that converts p-vector to n-vector, see Appendix B in Gade (2010):

n EB E , z EB =p_EB_E2n_EB_E( p EB E ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaWaamWaaeaaca WHUbWaa0baaSqaaiaadweacaWGcbaabaGaamyraaaaimqakiaa=Xca caaMe8UaamOEamaaBaaaleaacaWGfbGaamOqaaqabaaakiaawUfaca GLDbaacqGH9aqpjaaqqqaaaaaa8karHj3AcXwDLbWdbiaa=bhacaWF FbGaa8xraiaa=jeacaWFFbGaa8xraiaa=jdacaWFUbGaa83xaiaa=v eacaWFcbGaa83xaiaa=veak8aacaWFOaaeeG+aaaaaaivzKbWdciaa hchadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOWdaiaa=Lcaaa [email protected]@

 

2. Find latitude, longitude and height from 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOBamaaDa [email protected]@
 and 
z EB

la t EB ,lon g EB =n_E2lat_long( n EB E ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaWaamWaaeaaqa aaauFaaaWdbiaadYgacaWGHbGaamiDamaaBaaaleaacaWGfbGaamOq aaqabaacdeGcpaGaa8hla8qacaWGSbGaam4Baiaad6gacaWGNbWaaS baaSqaaiaadweacaWGcbaabeaaaOWdaiaawUfacaGLDbaacqGH9aqp jaaqqqaaaaaa8karHj3AcXwDLbWdciaa=5gacaWFFbGaa8xraiaa=j dacaWFSbGaa8xyaiaa=rhacaWFFbGaa8hBaiaa=9gacaWFUbGaa83z aOWdaiaa=HcacaWHUbWaa0baaSqaaiaadweacaWGcbaabaGaamyraa [email protected]@

h EB
 =  – 
z EB

Example 4: Geodetic latitude to ECEF-vector

Given geodetic latitude, longitude and height. Find the ECEF-vector (using WGS-84 ellipsoid).

Problem

See vector symbols explained for more details about the mathematical notation.
 
Geodetic latitude, longitude and height are given for position B as 
la t EB
, 
lon g EB
 and 
h EB
, find the ECEF-vector for this position, 
p EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
.
 

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
1. First, the given latitude and longitude are converted to n-vector:

n EB E =lat_long2n_E(la t EB ,lon g EB ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOBamaaDa aaleaacaWGfbGaamOqaaqaaiaadweaaaGccqGH9aqpimqajaaqqqaa aaaa8karHj3AcXwDLbWdbiaa=XgacaWFHbGaa8hDaiaa=9facaWFSb Gaa83Baiaa=5gacaWFNbGaa8Nmaiaa=5gacaWFFbGaa8xraOWdaiaa =Hcaqqa6daaaaaGuLrgapiGaamiBaiaadggacaWG0bWaaSbaaSqaai aadweacaWGcbaabeaak8aacaWFSaWdciaadYgacaWGVbGaamOBaiaa [email protected]@

 

2. We want
p EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCiCamaaDa [email protected]@
from
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCOBamaaDa [email protected]@
, and have functions that convert between these two. Now we need n_EB_E2p_EB_E.m. See also Appendix B in Gade (2010).

p EB E =n_EB_E2p_EB_E( n EB E , h EB ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHWbWaa0baaSqaaiaadweacaWGcbaabaGaamyraaaak8aacqGH 9aqpimqajaaqqqaaaaaa8karHj3AcXwDLbWdciaa=5gacaWFFbGaa8 xraiaa=jeacaWFFbGaa8xraiaa=jdacaWFWbGaa83xaiaa=veacaWF cbGaa83xaiaa=veak8aacaWFOaGaaCOBamaaDaaaleaacaWGfbGaam OqaaqaaiaadweaaaGccaWFSaGaaGjbVlabgkHiTabbOpaaaaaasvgz [email protected]@

In the remaining examples, we assume that a spherical Earth model is sufficiently accurate, and then the position functions used above are not needed. The examples demonstrate the simplicity of these calculations when using n-vector instead of latitude and longitude. The n-vector solutions are also non-singular for all Earth-positions and work across the Date Line (180th meridian).

Example 5: Surface distance

Given position A and B. Find the surface distance (i.e. great circle distance) and the Euclidean distance.

Problem

See vector symbols explained for more details about the mathematical notation.

Given two positions A and B as 

n EA E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaaaa @[email protected]
.

 
Find the surface distance 
s AB
 (i.e. great circle distance). The heights of A and B are not relevant (i.e. if they do not have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B). The Euclidean distance (chord length) 
d AB
 should also be found (for 
d AB
, nonzero heights can be used). Use Earth radius 
r Earth
.

 

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
The surface distance can be found in three different ways; by utilizing the dot product of the vectors, the cross product, or a combination:

 

s AB =acos n EA E n EB E r Earth =asin n EA E × n EB E r Earth =atan2 n EA E × n EB E , n EA E n EB E r Earth [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWGZbWaa0baaSqaaiaadgeacaWGcbaabaaaaOWdaiabg2da9GWa bKaaabbbaaaaaWRaefMCRjeB1vgapiGaa8xyaiaa=ngacaWFVbGaa8 3CaOWdamaabmaabaaeeG+aaaaaaivzKbWddiaah6gadaqhaaWcbaGa amyraiaadgeaaeaacaWGfbaaaOWdaiabgwSix=WacaWHUbWaa0baaS qaaiaadweacaWGcbaabaGaamyraaaaaOWdaiaawIcacaGLPaaacqGH flY1pmGaamOCamaaBaaaleaacaWGfbGaamyyaiaadkhacaWG0bGaam iAaaqabaGcpaGaeyypa0tcaa0dciaa=fgacaWFZbGaa8xAaiaa=5ga k8aadaqadaqaamaaemaabaWddiaah6gadaqhaaWcbaGaamyraiaadg eaaeaacaWGfbaaaOWdaiabgEna0+WacaWHUbWaa0baaSqaaiaadwea caWGcbaabaGaamyraaaaaOWdaiaawEa7caGLiWoaaiaawIcacaGLPa aacqGHflY1pmGaamOCamaaBaaaleaacaWGfbGaamyyaiaadkhacaWG 0bGaamiAaaqabaGcpaGaeyypa0tcaa0dciaa=fgacaWF0bGaa8xyai aa=5gacaWFYaGcpaWaaeWaaeaadaabdaqaa8WacaWHUbWaa0baaSqa aiaadweacaWGbbaabaGaamyraaaak8aacqGHxdaTpmGaaCOBamaaDa aaleaacaWGfbGaamOqaaqaaiaadweaaaaak8aacaGLhWUaayjcSdGa aiila8WacaWHUbWaa0baaSqaaiaadweacaWGbbaabaGaamyraaaak8 aacqGHflY1pmGaaCOBamaaDaaaleaacaWGfbGaamOqaaqaaiaadwea aaaak8aacaGLOaGaayzkaaGaeyyXIC9ddiaadkhadaWgaaWcbaGaam [email protected]@

 

The atan2 expression is numerically well-conditioned for all distances, see also equation (16) in Gade (2010).

 

The Euclidean distance is found by

d AB = n EA E n EB E r Earth [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWGKbWaa0baaSqaaiaadgeacaWGcbaabaaaaOWdaiabg2da9maa emaabaaeeG+aaaaaaivzKbWdciaah6gadaqhaaWcbaGaamyraiaadg eaaeaacaWGfbaaaOWdaiabgkHiT8GacaWHUbWaa0baaSqaaiaadwea caWGcbaabaGaamyraaaaaOWdaiaawEa7caGLiWoacqGHflY1piGaam OCamaaBaaaleaacaWGfbGaamyyaiaadkhacaWG0bGaamiAaaqabaaa [email protected]@
   

If nonzero heights for A and B are wanted for the Euclidean distance, each n-vector should be multiplied with 
r Earth
h i
, where 
h i
 is the corresponding height of each n-vector. If ellipsoidal Earth model is wanted, use n_EA_E_and_n_EB_E2p_AB_E.m from Example 1, and take the norm of the output vector.

Example 6: Interpolated position

Given the position of B at time
t 0
 and
t 1
. Find an interpolated position at time
t i
 .

Problem

See vector symbols explained for more details about the mathematical notation.
 
Given the position of B at time 
t 0
 and
t 1
n EB E ( t 0 ) MathTy[email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOGa [email protected]@
 and 
n EB E ( t 1 ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOGa [email protected]@
.
 
Find an interpolated position at time 
t i
n EB E ( t i ) [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOGa [email protected]@
. All positions are given as n-vectors.

 

 

Solution

Note: The solutions are also available as downloadable code in several different programming languages.
 
Standard interpolation can be used directly with n-vector:

n EB E ( t i )=unit n EB E ( t 0 )+ t i t 0 n EB E ( t 1 ) n EB E ( t 0 ) t 1 t 0 [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbiqaaqzeqaaaau FaaaWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOGa aiikaiaadshadaWgaaWcbaGaamyAaaqabaGccaGGPaWdaiabg2da9G WabKaaabbbaaaaaWRaefMCRjeB1vgapiGaa8xDaiaa=5gacaWFPbGa a8hDaOWdamaabmaabaaeeG+aaaaaaivzKbWddiaah6gadaqhaaWcba GaamyraiaadkeaaeaacaWGfbaaaOGaaiikaiaadshadaWgaaWcbaGa aGimaaqabaGccaGGPaWdaiabgUcaRmaabmaabaWddiaadshadaWgaa WcbaGaamyAaaqabaGcpaGaeyOeI0YddiaadshadaWgaaWcbaGaaGim aaqabaaak8aacaGLOaGaayzkaaWaaSaaaeaapmGaaCOBamaaDaaale aacaWGfbGaamOqaaqaaiaadweaaaGccaGGOaGaamiDamaaBaaaleaa caaIXaaabeaakiaacMcapaGaeyOeI0Yddiaah6gadaqhaaWcbaGaam yraiaadkeaaeaacaWGfbaaaOGaaiikaiaadshadaWgaaWcbaGaaGim aaqabaGccaGGPaaapaqaa8WacaWG0bWaaSbaaSqaaiaaigdaaeqaaO WdaiabgkHiT8WacaWG0bWaaSbaaSqaaiaaicdaaeqaaaaaaOWdaiaa [email protected]@

Example 7: Mean position/center

Given three positions A, B, and C. Find the mean position (center/midpoint).

Problem 7

See vector symbols explained for more details about the mathematical notation.
 
Three positions AB, and C are given as n-vectors 
n EA E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
, and 
n EC E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
. Find the mean position, M, given as 
n EM E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
. Note that the calculation is independent of the heights/depths of the positions.

 

Solution 7

Note: The solutions are also available as downloadable code in several different programming languages.
 
The mean position is simply given by the mean n-vector:

n EM E =unit n EA E + n EB E + n EC E = n EA E + n EB E + n EC E n EA E + n EB E + n EC E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHUbWaa0baaSqaaiaadweacaWGnbaabaGaamyraaaak8aacqGH 9aqpimqajaaqqqaaaaaa8karHj3AcXwDLbWdciaa=vhacaWFUbGaa8 xAaiaa=rhak8aadaqadaqaaabbOpaaaaaasvgza8WacaWHUbWaa0ba aSqaaiaadweacaWGbbaabaGaamyraaaak8aacqGHRaWkpmGaaCOBam aaDaaaleaacaWGfbGaamOqaaqaaiaadweaaaGcpaGaey4kaSYddiaa h6gadaqhaaWcbaGaamyraiaadoeaaeaacaWGfbaaaaGcpaGaayjkai aawMcaaiabg2da9maalaaabaWddiaah6gadaqhaaWcbaGaamyraiaa dgeaaeaacaWGfbaaaOWdaiabgUcaR8WacaWHUbWaa0baaSqaaiaadw eacaWGcbaabaGaamyraaaak8aacqGHRaWkpmGaaCOBamaaDaaaleaa caWGfbGaam4qaaqaaiaadweaaaaak8aabaWaaqWaaeaapmGaaCOBam aaDaaaleaacaWGfbGaamyqaaqaaiaadweaaaGcpaGaey4kaSYddiaa h6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaOWdaiabgUcaR8 WacaWHUbWaa0baaSqaaiaadweacaWGdbaabaGaamyraaaaaOWdaiaa [email protected]@

See also equation (17) in Gade (2010).
 
n-vector Red: n-vectors for positions A, B, and C. Green: n-vector for the mean position.
Red: n-vectors for positions A, B, and C. Green: n-vector for the mean position.

Details about the definition of the horizontal geographical mean position

Several definitions of a geographical mean/center/midpoint are possible. The solution given here can be described as the “center of gravity” of the given positions, and can be compared to the centroid of a geometrical shape. Informally, the found mean position can be described by the following: Assume a sphere that is a model of the Earth, and let the sphere roll freely on a horizontal plane (in a uniform gravitational field). Place weights at the given positions of the sphere, and the weights will pull one side of the sphere down due to the gravitation (assuming the positions are not antipodal). When the sphere is in equilibrium, the “center of gravity”-position is the tangent point between the sphere and the plane.

Another possible way to define a midpoint is the position where the sum of surface distances (great circle distances) from the original positions is at a minimum. This midpoint is undefined when only two positions are given. If three positions are given in one dimension as 0, 0, and 3, this midpoint is at 0, while the arithmetic mean is 1. Iterations are probably needed to calculate this midpoint, and code for this is not included at this web site.

Example 8: A and azimuth/distance to B

Given position A and an azimuth/bearing and a (great circle) distance. Find the destination point B.
Problem 8
See vector symbols explained for more details about the mathematical notation.
 
Position A is given as n-vector 
n EA E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
. We also have an initial direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle is given as 
s AB
. Use Earth radius 
r Earth
.

 

Find the destination point B, given as 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaaaa @[email protected]
.

 

In geodesy, this is known as "The first geodetic problem" or "The direct geodetic problem" for a sphere, and we see that this is similar to Example 2, but now the delta is given as an azimuth and a great circle distance. "The second/inverse geodetic problem" for a sphere is already solved in Examples 1 and 5.

 

Solution 8

Note: The solutions are also available as downloadable code in several different programming languages.
 
The azimuth (relative to north) is a singular quantity (undefined at the Poles), but from this angle we can find a (non-singular) quantity that is more convenient when working with vector algebra: a vector 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 that points in the initial direction. We find this from azimuth by first finding the north and east vectors at the start point, with unit lengths (see also equations (9) and (10) in Gade (2010), and the figure below):
 
k east E =unit 0 0 1 × n EA E k north E = n EA E × k east E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGceaqabeaacaWHRb Waa0baaSqaaiaadwgacaWGHbGaam4CaiaadshaaeaacaWGfbaaaOGa eyypa0dcdeqcaaeeeaaaaaaVcquyYTMqSvxza8qacaWF1bGaa8NBai aa=LgacaWF0bGcpaWaaeWaaeaadaWadaqaauaabeqadeaaaeaacaaI WaaabaGaaGimaaqaaiaaigdaaaaacaGLBbGaayzxaaGaey41aqleeG +aaaaaaivzKbWdciaah6gadaqhaaWcbaGaamyraiaadgeaaeaacaWG fbaaaaGcpaGaayjkaiaawMcaaaqaaiaahUgadaqhaaWcbaGaamOBai aad+gacaWGYbGaamiDaiaadIgaaeaacaWGfbaaaOGaeyypa0Zdciaa h6gadaqhaaWcbaGaamyraiaadgeaaeaacaWGfbaaaOWdaiabgEna0k aahUgadaqhaaWcbaGaamyzaiaadggacaWGZbGaamiDaaqaaiaadwea [email protected]@
 
Here we have assumed that our coordinate frame E has its z-axis along the rotational axis of the Earth, pointing towards the North Pole. Hence, this axis is given by 
0 0 1 [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaWaamWaaeaafa qabeWabaaabaGaaGimaaqaaiaaicdaaeaacaaIXaaaaaGaay5waiaa [email protected]@
.
 
nvector globe example 8
The north and east vectors are simple to find from the n-vector. The gray Earth-rotation-axis vector is also drawn at position A for convenience (vectors can be drawn at any position since they only have direction and magnitude).
 
The two vectors 
k north E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaC4AamaaDa aaleaacaWGUbGaam4BaiaadkhacaWG0bGaamiAaaqaaiaadweaaaaa [email protected]@
 and 
k east E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaC4AamaaDa [email protected]@
 are horizontal, orthogonal, and span the tangent plane at the initial position. A unit vector 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 in the direction of the azimuth is now given by:

d E = k north E cos azimuth + k east E sin azimuth [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa aaleqabaGaamyraaaakiabg2da9iaahUgadaqhaaWcbaGaamOBaiaa d+gacaWGYbGaamiDaiaadIgaaeaacaWGfbaaaOGaeyyXICTaci4yai aac+gacaGGZbWaaeWaaabbOpaaaaaasvgza8qabaGaamyyaiaadQha caWGPbGaamyBaiaadwhacaWG0bGaamiAaaWdaiaawIcacaGLPaaacq GHRaWkcaWHRbWaa0baaSqaaiaadwgacaWGHbGaam4Caiaadshaaeaa caWGfbaaaOGaeyyXICTaci4CaiaacMgacaGGUbWaaeWaa8qabaGaam yyaiaadQhacaWGPbGaamyBaiaadwhacaWG0bGaamiAaaWdaiaawIca [email protected]@

Math graph

With the initial direction given as 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 instead of azimuth, it is now quite simple to find 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
. We know that 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 and 
n EA E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
 are orthogonal, and they will span the plane where 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
 will lie. Thus, we can use sin and cos in the same manner as above, with the angle travelled given by 
s AB
/
r Earth
:

n EB E = n EA E cos s AB r Earth + d E sin s AB r Earth [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHUbWaa0baaSqaaiaadweacaWGcbaabaGaamyraaaak8aacqGH 9aqpqqa6daaaaaGuLrgapiGaaCOBamaaDaaaleaacaWGfbGaamyqaa qaaiaadweaaaGcpaGaeyyXICTaci4yaiaac+gacaGGZbWaaeWaa8Ga baWdamaalaaapiqaaiaadohadaWgaaWcbaGaamyqaiaadkeaaeqaaa GcbaGaamOCamaaBaaaleaacaWGfbGaamyyaiaadkhacaWG0bGaamiA aaqabaaaaaGcpaGaayjkaiaawMcaaiabgUcaRiaahsgadaqhaaWcba aabaGaamyraaaakiabgwSixlGacohacaGGPbGaaiOBamaabmaapiqa a8aadaWcaaWdceaacaWGZbWaaSbaaSqaaiaadgeacaWGcbaabeaaaO qaaiaadkhadaWgaaWcbaGaamyraiaadggacaWGYbGaamiDaiaadIga [email protected]@

Math graph

 

Example 9: Intersection of two paths / triangulation

Given path A going through
A 1
 and
A 2
, and path B going through
B 1
 and
B 2
. Find the intersection of the two paths.
 
Variant: Each of the two paths is specified by start point and azimuth instead of two points, which corresponds to triangulation.
 
Credits: Some readers of the paper Gade (2010) have contacted the author and reported new examples (other than those included in the paper) they have solved by means of n-vector. The following example has been reported twice, and is thus included here. Thanks to Even Børhaug at Kongsberg Maritime and Chris Veness at Movable Type Ltd, for this example. Børhaug reported that the solution he found that was based on latitude and longitude used multiple pages of code, while Veness says that such a solution is "horribly complex". Using n-vector, the intersection is very simple to find.

 

Problem 9

See vector symbols explained for more details about the mathematical notation.
 
Define a path from two given positions (at the surface of a spherical Earth), as the great circle that goes through the two points (assuming that the two positions are not antipodal).
 
Path A is given by 
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 and 
n EA,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
, while path B is given by 
n EB,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 and 
n EB,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
.
 
Find the position C where the two paths intersect, as 
n EC E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHUbWaa0baaSqaaiaadweacaWGdbaabaGaamyr[email protected]@
.

 

Solution 9

Note: The solutions are also available as downloadable code in several different programming languages.
 
A convenient way to represent a great circle is by its normal vector (i.e. the normal vector to the plane containing the great circle). This normal vector is simply found by taking the cross product of the two n-vectors defining the great circle (path). Having the normal vectors to both paths, the intersection is now simply found by taking the cross product of the two normal vectors:

n EC E =±unit n EA,1 E × n EA,2 E × n EB,1 E × n EB,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHUbWaa0baaSqaaiaadweacaWGdbaabaGaamyraaaak8aacqGH 9aqpcqGHXcqSimqajaaqqqaaaaaa8karHj3AcXwDLbWdciaa=vhaca WFUbGaa8xAaiaa=rhak8aadaqadaqaamaabmaabaaeeG+aaaaaaivz KbWddiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqaai aadweaaaGcpaGaey41aq7ddiaah6gadaqhaaWcbaGaamyraiaadgea caGGSaGaaGOmaaqaaiaadweaaaaak8aacaGLOaGaayzkaaGaey41aq 7aaeWaaeaapmGaaCOBamaaDaaaleaacaWGfbGaamOqaiaacYcacaaI XaaabaGaamyraaaak8aacqGHxdaTpmGaaCOBamaaDaaaleaacaWGfb GaamOqaiaacYcacaaIYaaabaGaamyraaaaaOWdaiaawIcacaGLPaaa [email protected]@

Note that there will be two places where the great circles intersect, and thus two solutions are found. Selecting the solution that is closest to e.g. 
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 can be achieved by selecting the solution that has a positive dot product with 
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 (or the mean position from Example 7 could be used instead of 
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
). In program code, this can be implemented by selecting any of the two solutions above (+ or –), called C', and calculate:
 
n EC E =sign n EC' E n EA,1 E n EC' E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWHUbWaa0baaSqaaiaadweacaWGdbaabaGaamyraaaak8aacqGH 9aqpimqajaaqqqaaaaaa8karHj3AcXwDLbWdciaa=nhacaWFPbGaa8 3zaiaa=5gak8aadaqadaqaa8qacaWHUbWaa0baaSqaaiaadweacaWG dbGaai4jaaqaaiaadweaaaGcpaGaeyyXICneeG+aaaaaaivzKbWddi aah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqaaiaadwea aaaak8aacaGLOaGaayzkaaWdbiaah6gadaqhaaWcbaGaamyraiaado [email protected]@

 

Variant / triangulation: If each of the two paths is specified by start point and azimuth instead of two points, the normal vectors to the great circles can be found by the following:
1. Find the direction vector for each path in the same manner as 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 in Example 8.
2. Take the cross product of the direction vector and the n-vector for the initial position.

Example 10: Cross track distance (cross track error)

Given path A going through
A 1
 and
A 2
, and a point B. Find the cross track distance/cross track error between B and the path.
 
Variant: The path is specified by start point and azimuth instead of two points.
 
(Along track distance (and the corresponding point on the path) can also be found quite easily)
 
Credits: Thanks to Chris Veness at Movable Type Ltd who first reported this example.

 

Problem 10

See vector symbols explained for more details about the mathematical notation.
 
Path A is given by the two n-vectors 
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 and 
n EA,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 (as in the previous example), and a position B is given by 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 [email protected]@
.
 
Find the cross track distance 
s xt
between the path A (i.e. the great circle through
n EA,1 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
 and 
n EA,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadgeacaGGSaGaaGymaaqa [email protected]@
) and the position B (i.e. the shortest distance at the surface, between the great circle and B). Also, find the Euclidean distance 
d xt
 between B and the plane defined by the great circle. Use Earth radius 
r Earth
.

 

Solution 10

Note: The solutions are also available as downloadable code in several different programming languages.
 
First, find the normal 
c E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 to the great circle, with direction given by the right hand rule and the direction of travel:
 
c E =unit n EA,1 E × n EA,2 E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaC4yamaaDa aaleaaaeaacaWGfbaaaOGaeyypa0dcdeqcaaeeeaaaaaaVcquyYTMq Svxza8qacaWF1bGaa8NBaiaa=LgacaWF0bGcpaWaaeWaaeaaqqa6da aaaaGuLrgapiGaaCOBamaaDaaaleaacaWGfbGaamyqaiaacYcacaaI XaaabaGaamyraaaak8aacqGHxdaTpiGaaCOBamaaDaaaleaacaWGfb GaamyqaiaacYcacaaIYaaabaGaamyraaaaaOWdaiaawIcacaGLPaaa [email protected]@

 

Math graph

We now see that the cross track distance to B is the great circle distance from 
c E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 minus 90° (π/2 radians):
 
s xt = acos c E n EB E π 2 r Earth =asin c E nEBErEarth[email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWGZbWaa0baaSqaaiaadIhacaWG0baabaaaaOWdaiabg2da9maa bmaabaacdeqcaaeeeaaaaaaVcquyYTMqSvxza8GacaWFHbGaa83yai aa=9gacaWFZbGcpaWaaeWaaeaacaWHJbWaa0baaSqaaaqaaiaadwea aaGccqGHflY1qqa6daaaaaGuLrgapmGaaCOBamaaDaaaleaacaWGfb GaamOqaaqaaiaadweaaaaak8aacaGLOaGaayzkaaGaeyOeI0YaaSaa aeaacqaHapaCaeaacaaIYaaaaaGaayjkaiaawMcaaiabgwSix=Waca WGYbWaaSbaaSqaaiaadweacaWGHbGaamOCaiaadshacaWGObaabeaa [email protected]@
 
Finding the Euclidean distance is even simpler, since it is the projection of 
n EB E [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeeG+aaaaaai vzKbWdbiaah6gadaqhaaWcbaGaamyraiaadkeaaeaacaWGfbaaaaaa @[email protected]
 onto 
c E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
, thus simply the dot product:

d xt = c E n EB E r Earth [email protected]@[email protected]@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaaeaaaa1haaa8 qacaWGKbWaa0baaSqaaiaadIhacaWG0baabaaaaOWdaiabg2da9iab gkHiTiaahogadaqhaaWcbaaabaGaamyraaaakiabgwSixdbbOpaaaa aasvgza8GacaWHUbWaa0baaSqaaiaadweacaWGcbaabaGaamyraaaa k8aacqGHflY1piGaamOCamaaBaaaleaacaWGfbGaamyyaiaadkhaca [email protected]@

For both 
s xt
 and 
d xt
, positive answers means that B is to the right of the track.
 
 
Variant: If the path is specified by start point and azimuth instead of two points, the normal vector to the great circle can be found by the following:
1. Find the direction vector for the path in the same manner as 
d E [email protected]@[email protected]@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFv0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaaCizamaaCa [email protected]@
 in Example 8.
2. Take the cross product of the direction vector and the n-vector for the initial position.
 
 
Along track distance to the closest point on path A can be found by using the four vectors in this example, and then calculating the intersection point as in Example 9 (when the n-vector for the intersection point is found, the along track distance from 
A 1
 is found as in Example 5). Thanks to Enrico Spinielli at EUROCONTROL for reporting this.
 

About the solutions

The four first examples are simply solved by using the functions provided for download. The functions are exact, and can use ellipsoidal or spherical Earth model. They are without iterations (i.e. on closed form), and WGS-84 reference ellipsoid is used as default. Custom ellipsoids/spheres can be specified. The functions are based on n-vector instead of latitude and longitude.

The last six examples are solved for spherical Earth, and then no functions are needed. Using n-vector instead of latitude and longitude, these examples are quite straightforward to solve by simple vector algebra.

The solutions to all 10 examples are non-singular, and they work equally well for all global positions (i.e. they work at the Poles, close to the Poles, and across the Date Line). The solutions also have full numerical accuracy for both short and long distances (they are numerically well-conditioned).