Distance between two points on the sphere.
78 views (last 30 days)
Show older comments
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 :)
0 Comments
Accepted Answer
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
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.
More Answers (2)
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
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.
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))
% W. Kahan formula
theta = 2 * atan(norm(x*norm(y) - norm(x)*y) / norm(x * norm(y) + norm(x) * y))
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))
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
Multiplying theta by the sphere radius r, you obtain then the geodesic distance between x and y.
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!