11
$\begingroup$

Suppose I have two points on a unit sphere whose spherical coordinates are $(\theta_1, \varphi_1)$ and $(\theta_2, \varphi_2)$. What is the great arc distance between these two points?

I found something from Wiki here but it is written in terms of latitude and longitude instead. I wonder latitude is $\pi/2 - \theta$ and longitude is $\varphi$. Am I right?

spherical coordinates

  • 3
    That is correct.2012-11-06
  • 0
    @RobertIsrael So I can just plug these into the equation on Wiki and compute the great arc distance, right?2012-11-06
  • 0
    latitude $\theta$ is measured from equator, not from the north pole as in your figure2013-05-08

3 Answers 3

12

Consider the two vectors to the points on the sphere, $${\bf v}_i=(r\sin\theta_i\cos\varphi_i,r\sin\theta_i\sin\varphi_i,r\cos\theta_i)$$ with $i=1,2$. Use the dot product to get the angle $\psi$ between them: $${\bf v}_1\cdot {\bf v}_2=r^2\left(\cos\theta_1\cos\theta_2+\sin\theta_1\sin\theta_2\cos\left(\varphi_1-\varphi_2\right)\right)=r^2\cos\psi.$$

Then the arclength is $$s=r\psi=r\cos^{-1}\left(\cos\theta_1\cos\theta_2+\sin\theta_1\sin\theta_2\cos\left(\varphi_1-\varphi_2\right)\right).$$

  • 0
    how exactly you got the dot product ?2014-02-28
  • 0
    What you get is not the formula from here: http://en.wikipedia.org/wiki/Great-circle_distance#Formulas . I can't understand how the differences in cos are derived since you compute the inner product2014-02-28
  • 0
    @curious the Wikipedia coordinates are defined differently than those in the question posted here.2014-02-28
  • 0
    Even if without the wikipedia definition i cant see how the formula of inner product between $v_1$ and $v_2$ gives the result you wrote2014-03-01
  • 2
    @curious $\sin^2\alpha+\cos^2\alpha=1$ and $\cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta)$.2014-03-01
  • 1
    I'm a bit confused, are \varphi_1 and \varphi_2 angles? If I use 0 for both \theta_1 and \theta_2, I get a 0 degree difference no matter what \varphi_1 and \varphi_2 are. That doesn't seem right to me.2016-03-08
  • 0
    @BT if $\theta_1=\theta_2=0$ then both points are the north pole and there is no distance between them regardless of $\varphi_1$, $\varphi_2$. I follow the angle definitions made by the diagram in the question.2016-03-10
  • 0
    @Jonathan Got it, someone else helped me figure that out. I was misunderstanding how spherical coordinates worked. Thanks!2016-03-10
  • 0
    if the points are diametrically opposite lying on the equator i.e. $\theta_1=\theta_2=0$ & $\phi_1-\phi_2=\pi$ then distance should be $\pi r$ but your formula doesn't hold good.2016-08-28
  • 0
    @Bhaskara-III with the coordinates as defined here, the two points you are describing are both at the north pole and the formula correctly shows a distance of zero.2016-09-01
1

I compare two formulas in Matlab, in terms of precision and running time.

Consider the sphere, in Cartesian coordinates:

$x^2 + y^2 + z^2 = \rho^2$

or in spherical coordinates:

$x = \rho \sin(\theta) \cos(\varphi)\\ y = \rho \sin(\theta) \sin(\varphi)\\ z = \rho \cos(\theta)$

Choose a bunch of points on the sphere, for example:

N = 1e6 + 1; rho_sph = rand; phi_sph = rand(N, 1) * pi; theta_sph = rand(N, 1) * pi;  xyz_sph = [rho_sph * sin(theta_sph) .* cos(phi_sph), ...            rho_sph * sin(theta_sph) .* sin(phi_sph), ...            rho_sph * cos(theta_sph)]; 

Now calculate the distance between each couple of points (1,2), (2,3), ..., (N-1, N).

1º formula

tic aa = cross(xyz_sph(1:end-1, :), xyz_sph(2:end, :)); bb = dot(xyz_sph(1:end-1, :)', xyz_sph(2:end, :)'); aa = arrayfun(@(n) norm(aa(n, :)), 1:size(aa, 1)); d1 = rho_sph * atan2(aa, bb)'; toc % Elapsed time is 12.376257 seconds. 

2º formula

tic d2 = rho_sph * acos(cos(theta_sph(1:end-1)) .* cos(theta_sph(2:end)) + ...      sin(theta_sph(1:end-1)) .* sin(theta_sph(2:end)) .* ...      cos(phi_sph(1:end-1) - phi_sph(2:end))); toc % Elapsed time is 0.264629 seconds. 

The first formula is considered more precise, but it's slower; the second is less accurate at the poles, but it's much faster. Otherwise the two solutions are practically equivalent:

plot(d1 - d2, '.') 
  • 0
    Does the expense of the first formula decrease if you store `xyz_sph(2:end, :)` as an intermediate value instead of computing it twice?2018-06-16
0

well, i see Wikipedia says there is serious rounding error in the formula of top answer. i found on google a new formula with the derivation. It seems slightly complicated but it is correct one. you can calculate values of distance & compare with percentage errors in both the formula.
i hope it will help you