Distance between two points on the sphere.
61 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Piotr Samol
el 18 de En. de 2017
Editada: Bruno Luong
el 17 de Jul. de 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 comentarios
Respuesta aceptada
Stephen23
el 18 de En. de 2017
Editada: Stephen23
el 18 de En. de 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 comentarios
Stephen23
el 18 de En. de 2017
Editada: Stephen23
el 18 de En. de 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.
Más respuestas (2)
Andrei Bobrov
el 18 de En. de 2017
Editada: Andrei Bobrov
el 18 de En. de 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 comentarios
Andrei Bobrov
el 18 de En. de 2017
Editada: Andrei Bobrov
el 18 de En. de 2017
Hi Stephen! Thank you for your comment. Another my mistake. Here should be used asin instead acos.
Fixed 2.
Bruno Luong
el 17 de Jul. de 2022
Editada: Bruno Luong
el 17 de Jul. de 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 comentarios
Ver también
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!