Distance from Lat/Lng point to Minor Arc segment

Question

I need to calculate the shortest distance from a lat/lng GPS point P to a line segment described by 2 other lat/lng GPS points A and B.

'Cross-track distance' helps me to calculate the shortest distance between P and the great circle described by A and B.

However, this is not what I want. I need need the distance between P and the line segment of A-B, not the entire great circle.

I have used the following implementation from http://www.movable-type.co.uk/scripts/latlong.html

Formula:    dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R where: δ13 is (angular) distance from start point to third point θ13 is (initial) bearing from start point to third point θ12 is (initial) bearing from start point to end point R is the earth’s radius 

The following images hopefully demonstrate the problem I am trying to solve:

In the first image the Cross-Track distance, indicated by the green line is correct and indeed the shortest distance to the line segment AB.

In the second image the problem with cross-track distance is shown, In this case I would want the shortest distance to be the simple distance AP, but Cross-Track distance gives me the distance indicated by the red line.

How do I change my algoritm to take this into account, or check whether or not point X is within AB. Is it possible to do this computationally? Or is iterative the only possible (expensive) solution? (take N points along AB and calculate the min distance from P to all these points)

For simplicity purposes all lines in the images are straight. In reality, these are minor arcs on a great circle

Answer 1

First, some nomenclature:
Our arc is drawn from p1 to p2.
Our third point is p3.
The imaginary point that intersects the great circle is p4.
p1 is defined by lat1,lon1; p2 by lat2,lon2; etc.
dis12 is the distance from p1 to p2; etc.
bear12 is the bearing from p1 to p2; etc.
dxt is cross-track distance.
dxa is cross-arc distance, our goal!

Notice that the cross-track formula relies on the relative bearing, bear13-bear12

We have 3 cases to deal with.

Case 1: The relative bearing is obtuse. So, dxa=dis13.

Case 2.1: The relative bearing is acute, AND p4 falls on our arc. So, dxa=dxt.

Case 2.2: The relative bearing is acute,AND p4 falls beyond our arc. So, dxa=dis23

The algorithm:

Step 1: If relative bearing is obtuse, dxa=dis13
Done!
Step 2: If relative bearing is acute:
2.1: Find dxt.
2.3: Find dis12.
2.4: Find dis14.
2.4: If dis14>dis12, dxa=dis23.
Done!
2.5: If we reach here, dxa=abs(dxt)

MATLAB code:

function [ dxa ] = crossarc( lat1,lon1,lat2,lon2,lat3,lon3 ) %// CROSSARC Calculates the shortest distance in meters  %// between an arc (defined by p1 and p2) and a third point, p3. %// Input lat1,lon1,lat2,lon2,lat3,lon3 in degrees.     lat1=deg2rad(lat1); lat2=deg2rad(lat2); lat3=deg2rad(lat3);     lon1=deg2rad(lon1); lon2=deg2rad(lon2); lon3=deg2rad(lon3);      R=6371000; %// Earth's radius in meters     %// Prerequisites for the formulas     bear12 = bear(lat1,lon1,lat2,lon2);     bear13 = bear(lat1,lon1,lat3,lon3);     dis13 = dis(lat1,lon1,lat3,lon3);      diff = abs(bear13-bear12);     if diff > pi         diff = 2 * pi - diff;     end     %// Is relative bearing obtuse?     if diff>(pi/2)         dxa=dis13;     else         %// Find the cross-track distance.         dxt = asin( sin(dis13/R)* sin(bear13 - bear12) ) * R;          %// Is p4 beyond the arc?         dis12 = dis(lat1,lon1,lat2,lon2);         dis14 = acos( cos(dis13/R) / cos(dxt/R) ) * R;         if dis14>dis12             dxa=dis(lat2,lon2,lat3,lon3);         else             dxa=abs(dxt);         end        end end  function [ d ] = dis( latA, lonA, latB, lonB ) %DIS Finds the distance between two lat/lon points. R=6371000; d = acos( sin(latA)*sin(latB) + cos(latA)*cos(latB)*cos(lonB-lonA) ) * R; end  function [ b ] = bear( latA,lonA,latB,lonB ) %BEAR Finds the bearing from one lat/lon point to another. b=atan2( sin(lonB-lonA)*cos(latB) , ...     cos(latA)*sin(latB) - sin(latA)*cos(latB)*cos(lonB-lonA) ); end 

Sample outputs: Demonstrate all cases. See maps below.

>> crossarc(-10.1,-55.5,-15.2,-45.1,-10.5,-62.5) ans =    7.6709e+05 >> crossarc(40.5,60.5,50.5,80.5,51,69) ans =    4.7961e+05 >> crossarc(21.72,35.61,23.65,40.7,25,42) ans =    1.9971e+05 

Those same outputs on the map!:

Demonstrates case 1:

Demonstrates case 2.1:

Demonstrates case 2.2:

Credit to: http://www.movable-type.co.uk/scripts/latlong.html
for the formulas
and: http://www.darrinward.com/lat-long/?id=1788764
for generating the map images.