Distance between two points on the sphere.

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 :)

Answer 1

[Stephen23]

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:

https://www.mathworks.com/matlabcentral/newsreader/view_thread/278872

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.

[Piotr Samol]

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]

@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.

Answer 2

[Andrei Bobrov]

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));

[Andrei Bobrov]

fixed

[Stephen23]

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]

Hi Stephen! Thank you for your comment. Another my mistake. Here should be used asin instead acos.

Fixed 2.

Answer 3

[Bruno Luong]

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.