IDL之由方向和距离求坐标(二分逼近法)

15195775117

问题由来：从已知坐标往X方向L米处的坐标是多少？

目前，我只知道距离怎么求，求坐标没有现成函数，于是我利用二分法来逼近精确坐标。
思路是，我先选坐标原点和大约1万米外的点(近似为经纬度0.1)为直线两端，
然后选择端点的中间点为新端点，
然后新端点和旧端点继续折中出新端点，......，
直到其中一个端点与原点的距离在L±1m的精度范围内。

代码与解释：

;vec是单位矢量，代表方向,笛卡尔坐标，从正东逆时针，[0,1]表示正北
function neighbor,core,vec,distance

  ;不确保为双精度，返回值很可能是6位有效数字
  x1=double(core[0]);经度
  y1=double(core[1]);纬度
  x2=x1+0.1*vec[0];以大约1万米外的位置为初始“远端点”
  y2=y1+0.1*vec[1]
  mx=(x1+x2)/2;折中点
  my=(y1+y2)/2
  
  constraint=1;精确度，单位米，不达到这个精度循环不停止

  repeat begin
    d1=map_2points(x1,y1,core[0],core[1],/METERS);近端点与原点的距离
    d2=map_2points(x2,y2,core[0],core[1],/METERS);远端点与原点的距离
    ;如果折中点较远，则折中点为新的远端点：
    if(map_2points(mx,my,core[0],core[1],/METERS) gt distance+constraint)then begin
      x2=mx
      y2=my
      mx=(x1+x2)/2
      my=(y1+y2)/2
    endif
    ;如果折中点较近，则折中点为新的近端点：
    if(map_2points(mx,my,core[0],core[1],/METERS) lt distance+constraint)then begin
      x1=mx
      y1=my
      mx=(x1+x2)/2
      my=(y1+y2)/2
    endif
    ;如果远近2个端点其中之一与原点的距离在误差范围内就结束循环
  endrep until (abs(d1-distance) le constraint or abs(d2-distance) le constraint)


  d1=map_2points(x1,y1,core[0],core[1],/METERS)
  d2=map_2points(x2,y2,core[0],core[1],/METERS)
  x=d1 gt d2?x2:x1;二端点近者为返回值
  y=d1 gt d2?y2:y1
  return,[x,y]
end

测试：
IDL> neighbor([120.0,30.0],[0,1],500)
% Compiled module: NEIGHBOR.
% Compiled module: MAP_2POINTS.
       120.00000000000000       30.004492187566939

该算法的python版见《python之由方向和距离求坐标(二分逼近法)》：
http://bbs.06climate.com/forum.php?mod=viewthread&tid=96382&extra=