Can't reproduce distance value between sources obtained with astropy

Question

I have two sources with equatorial coordinates (ra, dec) and (ra_0, dec_0) located at distances r and r_0, and I need to calculate the 3D distance between them.

I use two approaches that should give the same result as far as I understand, but do not.

The first approach is to apply astropy's separation_3d function. The second approach is to use the expression that gives the distance between two sources with spherical coordinates:

as shown here.

In the MCVE below, the values returned are:

91.3427173002 pc
93.8470493776 pc

Shouldn't these two values be equal?

MCVE:

from astropy.coordinates import SkyCoord
from astropy import units as u
import numpy as np

# Define some coordinates and distances for the sources.
c1 = SkyCoord(ra=9.7*u.degree, dec=-50.6*u.degree, distance=1500.3*u.pc)
c2 = SkyCoord(ra=7.5*u.degree, dec=-47.6*u.degree, distance=1470.2*u.pc)

# Obtain astropy's distance between c1 & c2 coords.
print c1.separation_3d(c2)

# Obtain distance between c1 & c2 coords using explicit expression.
ra_0, dec_0, r_0 = c1.ra.radian, c1.dec.radian, c1.distance
ra, dec, r = c2.ra.radian, c2.dec.radian, c2.distance
alpha_delta_par = np.sin(dec) * np.sin(dec_0) * np.cos(ra - ra_0) +\
    np.cos(dec) * np.cos(dec_0)
d_pc = np.sqrt(r**2 + r_0**2 - 2*r*r_0*alpha_delta_par)
print d_pc

Comyns

I don't have the module installed, but throwing a suggestion. Are you certain SkyCoord() knowns exactly in what units you're feeding it the angle values? – Nineteen

It doesn't, which is why I explicit the units used using *u.degree and *u.pc. – Comyns

Answer 1

This is a problem with coordinate systems, and the difference between declination (astral coordinates) and polar angle θ (spherical coordinates) :-)

Astral coordinates define declination as north of the celestial equator, while spherical coordinates define polar angle θ as downward from from vertical.

If you change your alpha_delta_par to account for this 90° difference by adding np.pi/2 to all your declination terms, you get

alpha_delta_par = np.sin(np.pi/2 + dec)*np.sin(np.pi/2 + dec0)*np.cos(ra - ra0) +\
np.cos(np.pi/2 + dec)*np.cos(np.pi/2 + dec0)

Which gives the correct result: 91.3427173002 pc.

_{^{Turns out physicists usually use the symbol θ as the polar angle and mathemeticians usually use φ; I went with θ because I followed my heart. I'm not making this up I swear.}}

Blaineblainey