Distance between two points on the sphere.
61 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Piotr Samol
am 18 Jan. 2017
Bearbeitet: Bruno Luong
am 17 Jul. 2022
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 Kommentare
Akzeptierte Antwort
Stephen23
am 18 Jan. 2017
Bearbeitet: Stephen23
am 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 Kommentare
Stephen23
am 18 Jan. 2017
Bearbeitet: Stephen23
am 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.
Weitere Antworten (2)
Andrei Bobrov
am 18 Jan. 2017
Bearbeitet: Andrei Bobrov
am 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 Kommentare
Andrei Bobrov
am 18 Jan. 2017
Bearbeitet: Andrei Bobrov
am 18 Jan. 2017
Hi Stephen! Thank you for your comment. Another my mistake. Here should be used asin instead acos.
Fixed 2.
Bruno Luong
am 17 Jul. 2022
Bearbeitet: Bruno Luong
am 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 Kommentare
Siehe auch
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!