第十七章 角距离 Angular Separation
1. 定义
角距离，也称为角分离、视距离、或视分离，在数学(特别是几何学和三角学)和自然科学(包括天文学、地质学等等)，从不同于两个点物体的位置（即第三点）观察这两个物体，由观测者指向这两个物体的直线之间所夹角度的大小。角距离(或分离)与角度本身是同义的，但意义却是对两个天体(对恒星，是当从地球观测)之间线距离的建议(通常是很大或未知的)。
2. 计算
2.1 余弦公式直接计算
$$\cos d=\sin δ_1 \sin δ_2 + \cos δ_1 \cos δ_2 \cos (α_1 -α_2)$$
式中$α_1,δ_1,α_2,δ_2$分别对应两个天体的赤经和赤纬。
当$d$接近于0或180度时，$\left|\cos d\right|$接近于1，并且其值随$d$变化很小，所以得到的$d$不精确。此时需用以下公式计算：$$d = \sqrt {(Δα\cos δ)^2+(Δδ)^2}$$
式中$Δα$是两个赤经的差，$Δδ$是两个赤纬的差，$δ$是两个赤纬的平均值。
// Sep returns the angular separation between two celestial bodies.
//
// The algorithm is numerically naïve, and while patched up a bit for
// small separations, remains unstable for separations near π.
// 计算两天体之间的角距，r为赤经，d为赤纬
func Sep(r1, d1, r2, d2 unit.Angle) unit.Angle {
sd1, cd1 := d1.Sincos()
sd2, cd2 := d2.Sincos()
cd := sd1*sd2 + cd1*cd2*(r1-r2).Cos() // (17.1) p. 109
if cd < base.CosSmallAngle {
return unit.Angle(math.Acos(cd))
}
// (17.2) p. 109
dm := (d1 + d2) / 2
return unit.Angle(math.Hypot((r2-r1).Rad()*dm.Cos(), (d2 - d1).Rad()))
}
2.2 根据观测数据插值计算最小角距离
注意，不能先计算出各组数据点的角距离，再对角距离插值计算。因为当两个天体很近时，两天体间是线性的靠近再线性的离开。此时需要先对数据点插值，再把得出的数据套用2.1的方法进行计算。或者引入 u,v 坐标，先将原始数据转为 u,v 坐标，然后对 u,v 坐标点进行插值计算，得出最小值
// MinSep returns the minimum separation between two moving objects.
//
// The motion is represented as an ephemeris of three rows, equally spaced
// in time. Jd1, jd3 are julian day times of the first and last rows.
// R1, d1, r2, d2 are coordinates at the three times. They must each be
// slices of length 3.
//
// Result is obtained by computing separation at each of the three times
// and interpolating a minimum. This may be invalid for sufficiently close
// approaches.
//
// 计算两个天体之间的最小角距
// 此方法是将数据点计算成角距，然后直接对角距3点插值，求取最小值，
// 当两个天体十分接近时，这个结果是不准确的，要使用方法MinSepRect
func MinSep(jd1, jd3 float64, r1, d1, r2, d2 []unit.Angle) (unit.Angle, error) {
if len(r1) != 3 || len(d1) != 3 || len(r2) != 3 || len(d2) != 3 {
return 0, interp.ErrorNot3
}
y := make([]float64, 3)
for x, r := range r1 {
y[x] = Sep(r, d1[x], r2[x], d2[x]).Rad()
}
d3, err := interp.NewLen3(jd1, jd3, y)
if err != nil {
return 0, err
}
_, dMin, err := d3.Extremum()
return unit.Angle(dMin), err
}
// MinSepRect returns the minimum separation between two moving objects.
//
// Like MinSep, but using a method of rectangular coordinates that gives
// accurate results even for close approaches.
// 计算两个天体之间的最小角距
// 此方法是引入直角坐标 u,v,先将数据点转换成u,v 表达，然后对u,v插值，求取最小值，
func MinSepRect(jd1, jd3 float64, r1, d1, r2, d2 []unit.Angle) (unit.Angle, error) {
if len(r1) != 3 || len(d1) != 3 || len(r2) != 3 || len(d2) != 3 {
return 0, interp.ErrorNot3
}
uv := func(r1, d1, r2, d2 unit.Angle) (u, v float64) {
sd1, cd1 := d1.Sincos()
Δr := r2 - r1
tΔr := Δr.Tan()
thΔr := (Δr / 2).Tan()
K := 1 / (1 + sd1*sd1*tΔr*thΔr)
sΔd := (d2 - d1).Sin()
u = -K * (1 - (sd1/cd1)*sΔd) * cd1 * tΔr
v = K * (sΔd + sd1*cd1*tΔr*thΔr)
return
}
us := make([]float64, 3, 6)
vs := us[3:6]
for x, r := range r1 {
us[x], vs[x] = uv(r, d1[x], r2[x], d2[x])
}
u3, err := interp.NewLen3(-1, 1, us)
if err != nil {
panic(err) // bug not caller's fault.
}
v3, err := interp.NewLen3(-1, 1, vs)
if err != nil {
panic(err) // bug not caller's fault.
}
up0 := (us[2] - us[0]) / 2
vp0 := (vs[2] - vs[0]) / 2
up1 := us[0] + us[2] - 2*us[1]
vp1 := vs[0] + vs[2] - 2*vs[1]
up := up0
vp := vp0
dn := -(us[1]*up + vs[1]*vp) / (up*up + vp*vp)
n := dn
var u, v float64
for limit := 0; limit < 10; limit++ {
u = u3.InterpolateN(n)
v = v3.InterpolateN(n)
if math.Abs(dn) < 1e-5 {
return unit.Angle(math.Hypot(u, v)), nil // success
}
up := up0 + n*up1
vp := vp0 + n*vp1
dn = -(u*up + v*vp) / (up*up + vp*vp)
n += dn
}
return 0, errors.New("MinSepRect: failure to converge")
}
2.3 利用半正矢的特点提高当角距很小时的精确程度
根据半正矢公式$hav(d) = hav(Δδ) + \cos δ_1 \cos δ_2 hav(Δα)$,式中 $Δα = α_1 - α_2，Δδ = δ_1 - δ_2$， 又由$hav(d)=\frac {1-\cos d}{2}=\sin^2(\frac d2)$,可以有效的在0，180°附近提高计算机的精度。
// SepHav returns the angular separation between two celestial bodies.
//
// The algorithm uses the haversine function and is superior to the naïve
// algorithm of the Sep function.
// 利用半正矢的特点提高当角距很小时的精确程度
func SepHav(r1, d1, r2, d2 unit.Angle) unit.Angle {
// using (17.5) p. 115
return unit.Angle(2 * math.Asin(math.Sqrt(base.Hav(d2-d1)+
d1.Cos()*d2.Cos()*base.Hav(r2-r1))))
}
2.4 Pauwels公式
令 \begin{cases} x &= \cos δ_1\sin δ_2 - \sin δ_1\cos δ_2\cos (α_2-α_1)\\[2ex] x &= \cos δ_2\sin (α_2-α_1)\\[2ex] z &= \sin δ_1\sin δ_2 + \cos δ_1\cos δ_2\cos (α_2-α_1)\\[2ex] \end{cases} 则有：$$d=\arctan (\frac {\sqrt {x^2+y^2}}{z})$$
数学上来说，这与余弦定理完全等价，只不过是巧妙的将余弦转化为正切，而对于计算机来说， 反正切比反正弦能获得更高的精确度
// SepPauwels returns the angular separation between two celestial bodies.
//
// The algorithm is a numerically stable form of that used in Sep.
// 当z小于0时，返回值应该在90-180度之间
// 该方法与直接余弦定理求角距在数学上是等价的，
// 但是对于计算机来说，arctan 比 arcsin能获得更高的精度
func SepPauwels(r1, d1, r2, d2 unit.Angle) unit.Angle {
sd1, cd1 := d1.Sincos()
sd2, cd2 := d2.Sincos()
cdr := (r2 - r1).Cos()
x := cd1*sd2 - sd1*cd2*cdr
y := cd2 * (r2 - r1).Sin()
z := sd1*sd2 + cd1*cd2*cdr
return unit.Angle(math.Atan2(math.Hypot(x, y), z))
}
2.5 Relative Position Angle
// RelativePosition returns the position angle of one body with respect to
// another.
//
// The position angle result is measured counter-clockwise from North.
// 1相对2的角距，从2的正北到1的角度
// https://en.wikipedia.org/wiki/Position_angle
func RelativePosition(r1, d1, r2, d2 unit.Angle) unit.Angle {
sΔr, cΔr := (r1 - r2).Sincos()
sd2, cd2 := d2.Sincos()
return unit.Angle(math.Atan2(sΔr, cd2*d1.Tan()-sd2*cΔr))
}