# How to create a point along a line given a distance from another point

I'm trying to figure out a point(D) on B-C that is a specific distance from A.
I really don't have any idea on how to proceed with the ellipsoidal math so I figured I'd ask here. How would I figure out the longitude, latitude of D which is 5mi(26400ft) away from A, on the line B-C? The bearing of A-D would be just as useful as i could figure out the point 5mi away from A with the bearing.

I wouldn't mind if the formula lined up with the Vincenty model of the earth.

Using the Vincenty formula's I've already calculated the following.

Points (longitude, latitude)
A: -86.32899412316736,41.12507719625815
B: -86.29237110757369,41.16845596588609
C: -86.19689279712584,41.16923568742915

A-B Distance: 18750ft
B-C Distance: 26289ft
A-B Initial Bearing: 32.5331889647
B-C Initial Bearing: 89.3493884487

Thanks for the help!

This is most easily solve using the equidistant azimuthal projection. Guess some point for D, e.g., A. Transform A, B, C, to equidistant azimuthal projection using D as a center. In projected space, solve for D' (i.e., D' lies on BC and is a distance s from A). Project D' back to lat,lon and update D with this position. Repeat. This converges quadratically.

To see this in action, grab http://www.mathworks.com/matlabcentral/fileexchange/39108 and http://www.mathworks.com/matlabcentral/fileexchange/39366 from MathWorks File Exchange.

Write a .m file to solve the planar problem, planesolve.m. Here is a not-very-elegant version

``````function D = planesolve(pts, s)
% pts is 3x2 matrix of [A;B;C] lying in a plane.  Find point D which lies
% on B-C and which is a distance s from A.  There may be 0, 1, or 2
% solutions.  Here assume there are two and pick the one most in the
% direction B-C.

xa=pts(1,1); ya=pts(1,2);
xb=pts(2,1); yb=pts(2,2);
xc=pts(3,1); yc=pts(3,2);
% write D = t*B + (1-t)*C, solve for t
t = ((yc-ya)*(yc-yb)+(xc-xa)*(xc-xb) + ...
sqrt((-xb^2+2*xa*xb-xa^2+s^2)*yc^2+...
(((2*xb-2*xa)*xc-2*xa*xb+2*xa^2-2*s^2)*yb+...
((2*xa-2*xb)*xc+2*xb^2-2*xa*xb)*ya)*yc+...
(-xc^2+2*xa*xc-xa^2+s^2)*yb^2+...
(2*xc^2+(-2*xb-2*xa)*xc+2*xa*xb)*ya*yb+...
(-xc^2+2*xb*xc-xb^2)*ya^2+s^2*xc^2-2*s^2*xb*xc+s^2*xb^2))/...
((yc-yb)^2+(xc-xb)^2);
xd=t*xb+(1-t)*xc; yd=t*yb+(1-t)*yc;
D=[xd,yd];
end
``````

Then run the following script (in either matlab or octave)

``````g=[41.12507719625815,-86.32899412316736;...
41.16845596588609,-86.29237110757369;...
41.16923568742915,-86.19689279712584];
s = 5*1760*36*0.0254;
gd=g(1,:);
for i=1:5,
[x,y] = eqdazim_fwd(gd(end,1),gd(end,2),g(:,1),g(:,2));
d = planesolve([x,y],s);
[lat,lon] = eqdazim_inv(gd(end,1),gd(end,2),d(1),d(2));
gd = [gd;lat,lon];
end
format long;
gd
``````

This prints out

``````gd =

41.125077196258147 -86.328994123167362
41.167417462129485 -86.406778034888063
41.167417401850848 -86.406778092275701
41.167417401850948 -86.406778092275616
41.167417401850948 -86.406778092275601
41.167417401850948 -86.406778092275601
``````

which gives the converging sequence of approximations (lat,lon) for D. (Here the WGS84 ellipsoid is assumed. eqdazi_fwd and eqdazi_inv take an optional argument which lets you specify the ellipsoid.)

• This ends up being really close, the distance is correct but It doesn't land on line segment B-C? – Phil May 26 '13 at 17:41
• Yes! The distance is right because the projection is equidistant (with D as the center). D lies on B-C because the projection is azimuthal (with D as the center). I should add that you can pretty much choose any point as a starting point: A, B, C, (B+C)/2, the North Pole, the Eiffel Tower all work fine. – cffk May 28 '13 at 12:42
• the distance is right.. but how would you go about getting a point between B & C. instead of just on the line B-C. – Phil May 31 '13 at 4:05
• My apologies, I got you the solution outside the range B-C. Change the sign of the sqrt in planesolve.m and rerun and you will get the solution you want 41.1687907834032 -86.2525504601801. In general you should modify the code to try both solutions in planesolve and then select the one (if any) which gives you opposite azimuths D->B and D->C. – cffk May 31 '13 at 9:36

Intersection of two paths given start points and bearings could be used. You would need to check the resulting point to make sure it is between B and C. Since there's no spheroid flattening, I'm not sure how close it would match up with Vincenty's formula.

``````Formula:

d12 = 2.asin( √(sin²(Δφ/2) + cos(φ1).cos(φ2).sin²(Δλ/2)) )
φ1 = acos( sin(φ2) − sin(φ1).cos(d12) / sin(d12).cos(φ1) )
φ2 = acos( sin(φ1) − sin(φ2).cos(d12) / sin(d12).cos(φ2) )

if sin(λ2−λ1) > 0
θ12 = φ1, θ21 = 2.π − φ2
else
θ12 = 2.π − φ1, θ21 = φ2

α1 = (θ1 − θ12 + π) % 2.π − π
α2 = (θ21 − θ2 + π) % 2.π − π

α3 = acos( −cos(α1).cos(α2) + sin(α1).sin(α2).cos(d12) )
d13 = atan2( sin(d12).sin(α1).sin(α2), cos(α2)+cos(α1).cos(α3) )
φ3 = asin( sin(φ1).cos(d13) + cos(φ1).sin(d13).cos(θ1) )
Δλ13 = atan2( sin(θ1).sin(d13).cos(φ1), cos(d13)−sin(φ1).sin(φ3) )
λ3 = (λ1+Δλ13+π) % 2.π − π
where

φ1, λ1, θ1 : 1st point & bearing
φ2, λ2, θ2 : 2nd point & bearing
φ3, λ3 : intersection point

% = mod
note –  if sin(α1)=0 and sin(α2)=0: infinite solutions
if sin(α1).sin(α2) < 0: ambiguous solution
this formulation is not always well-conditioned for meridional or equatorial lines
``````
• wouldn't that just give point B – Phil May 23 '13 at 20:07