已知球面经纬度求方位角和反方位角（awk一行代码实现）

已知球面经纬度求方位角和距离

---

一个常见的错误

假如你在广州，先朝东北走2000km，然后朝西南走2000km，你不会回到起点，而是到达深圳或者东莞。
这是因为地球是一个球面，方位角和反方位角加起来并不是360度。所以要朝着南偏西约49度的方向走2000km才能回到原点。

---

原理

球面上两点和北极点组成球面三角形，黎曼几何告诉我们，这个三角形内角和大于180度。而且起点到终点的方位角和终点到起点的方位角只和也不是180度。
已知球面的起点（lon1，lat1）和终点（lon2，lat2），求他们的方位角azi的公式是：
cos(azi)=sin(lon1-lon2)cos(lat2)/(cos(lat1)sin(lat2)-sin(lat1)cos(lat2)cos(lon1-lon2))

---

awk一行代码实现求起点到终点的方位角

#以下代码输入经纬度单位是度，输出方位角的单位也是度
echo $lon1 $lat1 $lon2 $lat2 | awk '{
  radlon1=$1*0.017453;
  radlat1=$2*0.017453; 
  radlon2=$3*0.017453;
  radlat2=$4*0.017453; 
  deltalon=($3-$1)*0.017453; 
  y=sin(deltalon)*cos(radlat2); 
  x=cos(radlat1)*sin(radlat2)-sin(radlat1)*cos(radlat2)*cos(deltalon); 
  az=(atan2(y,x)/0.017453+360)%360; 
  print "Azimuth from (",$1,$2,") to (",$3,$4,") is",az," deg."}' | cat

求反方位角时，把起点重点的位置互换即可。

尤其注意
以上代码代表点2（lon2，lat2）相对于点1（lon1，lat1）的方位角，因此，da必须是lon2-lon1,不能反过来。