此文转载
例题:
已知巴黎的地理位置为 东经2°20'14'',北纬 48°50'11',华盛顿的地理位置为 西经 22°03'56'',北纬 38°55'17''
求两城市之间的距离。
解答:
用低精度和高精度两种方法解答,两种方法的区别在于,低精度将地球看作球形,高精度将地球看作椭球。
首先将经纬度转换成十进制的度。
定义函数deg2dec,用来完成这一转换
deg2dec <-
function(h,m,s){
if(h < 0)
{m = - m
s = -s
}
res = h + m/60 + s/3600
return(res)
}
##巴黎
L1 = deg2dec(-2,20,14)
phi1 = deg2dec(48, 50, 11)
##华盛顿
L2 = deg2dec(77,03,56)
phi2 = deg2dec(38,55,17)
## 第一种低精度方案,
## 由于地球可以的扁率较小,可以近似看作一个球体,利用球面三角学公式可以快速的求出两点之间的距离。写成计算机程序,如下所示:
#低精度公式
geodist <-
function(L1, phi1, L2, phi2)
{
Ri = 6371
d = ( Ri * acos(sin(phi1*(pi/180))*sin(phi2*(pi/180)) + cos(phi1*(pi/180))*cos(phi2*(pi/180))*cos((L1 - L2)*(pi/180))))
res <- round(d, 3)
return(res)
}
#低精度
geodist(L1, phi1, L2, phi2)
答案:
华盛顿和巴黎之间的距离是 6165.597km
##第二种 高精度方案,高精度的方案的误差低于+_50m.
## 采用较为精密的公式,考虑到地球是一个椭球,a是地球的长轴,f是扁率,用如下公式计算
hageodist <-
function(L1, phi1, L2, phi2){
a = 6378.14
f = 1/298.257
F = (phi1+phi2)/2
G = (phi1 - phi2)/2
ramda <- (L1 - L2)/2
S = (sin(G*pi/180)^2)*(cos(ramda*pi/180)^2) + (cos(F*pi/180)^2)*(sin(ramda*pi/180)^2)
C= (cos(G*pi/180)^2)*(cos(ramda*pi/180)^2) + (sin(F*pi/180)^2)*(sin(ramda*pi/180)^2)
omega = atan(sqrt(S/C))
R = sqrt(S*C)/omega
D = 2*omega*a
H1 = (3*R-1)/(2*C)
H2 = (3*R+1)/(2*S)
res = D*(1 + f*H1*(sin(F*pi/180)^2)*(cos(G*pi/180)^2) - f*H2*(cos(F*pi/180)^2)*(sin(G*pi/180)^2))
return(round(res,3))
}
hageodist(L1, phi1, L2, phi2)
答案:
华盛顿和巴黎之间的距离是6181.628km.
在计算这一距离时,高精度和低精度算法给出的结果之间差距达16km。可见在天文计算中,低精度是远远不能满足要求的。
例题和公式参考 Jean Meeus 的 Astronomical Algorithms