# Distance between two points on the sphere.

42 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

### 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

### Categories

### Products

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!