Distance between two points on the sphere.
61 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Stephen23
2017 年 1 月 18 日
編集済み: Stephen23
2017 年 1 月 18 日
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 件のコメント
Stephen23
2017 年 1 月 18 日
編集済み: Stephen23
2017 年 1 月 18 日
@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.
その他の回答 (2 件)
Andrei Bobrov
2017 年 1 月 18 日
編集済み: Andrei Bobrov
2017 年 1 月 18 日
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 件のコメント
Andrei Bobrov
2017 年 1 月 18 日
編集済み: Andrei Bobrov
2017 年 1 月 18 日
Hi Stephen! Thank you for your comment. Another my mistake. Here should be used asin instead acos.
Fixed 2.
Bruno Luong
2022 年 7 月 17 日
編集済み: Bruno Luong
2022 年 7 月 17 日
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 件のコメント
参考
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!