Distance between two points on the sphere.

78 views (last 30 days)
Greetings!
I have to make a script that build a sphere (radius is given by me), then the user inputs two coordinates (x,y,z) ON the sphere, and script shows the closest distance between these two points.
I have no clue how to do that (I was not taught such things on my classes), even though I have to do this...
I wish you help me :)

Accepted Answer

Stephen23
Stephen23 on 18 Jan 2017
Edited: Stephen23 on 18 Jan 2017
I guess you want to find the shortest distance along the surface of the sphere, not just the euclidean distance between the points. The shortest path between two points on a sphere is always located on a great circle, which is thus a "great arc". You can find Roger Staffords mathematically robust method here:
For convenience I repeat it here too:
" P1 = [x1,y1,z1] and P2 = [x2,y2,z2] are two vectors pointing from the center of the sphere to the two given points (x1,y1,z1) and (x2,y2,z2) on the sphere, what is the shortest great circle distance d between them?"
d = radius * atan2(norm(cross(P1,P2)),dot(P1,P2));
If you want an example of how to use this to generate points along a great arc, then see the Mfile colornames_cube in my FEX submission colornames. The nested function cncDemoClBk at the end of the file steps along a great arc, rotating the axes as it goes.
  2 Comments
Piotr Samol
Piotr Samol on 18 Jan 2017
Edited: Piotr Samol on 18 Jan 2017
Great it works. What about estimating 'the closest distance' (in the straight line through the globe) between for example two cities if u know their geographical coordinates? Radius would be equal to radius of the Earth and how to replace spherical coordinates x1,y1,z1 ... to geographical coords?
Stephen23
Stephen23 on 18 Jan 2017
Edited: Stephen23 on 18 Jan 2017
@Piotr Samol: converting between coordinate systems is a totally different topic! As I said, have a look at cncDemoClBk in the Mfile colornames_cube: this also converts between spherical and Cartesian coordinates using inbuilt Cartesian conversions.
You will find other conversions in the mapping toolbox or on FEX.

Sign in to comment.

More Answers (2)

Andrei Bobrov
Andrei Bobrov on 18 Jan 2017
Edited: Andrei Bobrov on 18 Jan 2017
Let c - your two coordinates ( c = [x1 y1 z1; x2 y2 z2] ), r - radius your sphere
distance_out = 2*r*asin(norm(diff(c))/2/r);
Use geographical coordinates:
P1 - coordinates of the first point ( P1 = [Lat1, Lon1] )
P2 - coordinates of the second point ( P2 = [Lat2, Lon2] )
R - the radius of the Earth.
P = [P1;P2];
C = abs(diff(P(:,2)));
a = 90 - P(:,1);
cosa = prod(cosd(a)) + prod(sind(a))*cosd(C);
distance_out = sqrt(2*R^2*(1 - cosa));
  3 Comments
Stephen23
Stephen23 on 18 Jan 2017
Edited: Stephen23 on 18 Jan 2017
Opposite points:
>> P1 = [0,0,-1];
>> P2 = [0,0,+1];
>> c = [P1;P2];
>> r = 1;
>> 2*r*acos(norm(diff(c))/2/r)
ans = 0
My answer:
>> atan2(norm(cross(P1,P2)),dot(P1,P2))
ans = 3.1416
Andrei Bobrov
Andrei Bobrov on 18 Jan 2017
Edited: Andrei Bobrov on 18 Jan 2017
Hi Stephen! Thank you for your comment. Another my mistake. Here should be used asin instead acos.
Fixed 2.

Sign in to comment.


Bruno Luong
Bruno Luong on 17 Jul 2022
Edited: Bruno Luong on 17 Jul 2022
Another formula of angle betwen two (3 x 1) vectors x and y that is also quite accurate is W. Kahan pointed by Jan here
% test vector
x = randn(3,1);
y = randn(3,1);
theta = atan2(norm(cross(x,y)),dot(x,y))
theta = 2.2235
% W. Kahan formula
theta = 2 * atan(norm(x*norm(y) - norm(x)*y) / norm(x * norm(y) + norm(x) * y))
theta = 2.2235
The advantage is for points on spherical surface, i.e., norm(x) = norm(y) = r , such as vectors obtained after thise normalization
r = 6371;
x = r * x / norm(x);
y = r * y / norm(y);
and the formula becomes very simple:
theta = 2 * atan(norm(x-y) / norm(x+y))
theta = 2.2235
or with more statements but less arithmetic operations
d = x-y;
s = x+y;
theta = 2*atan(sqrt((d'*d)/(s'*s))) % This returns correct angle even for s=0
theta = 2.2235
Multiplying theta by the sphere radius r, you obtain then the geodesic distance between x and y.

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!