前言：
雖然正統的 FORTH語言不太使用浮點運算，但 gFORTH 真的是個不錯的科學計算平台。它的浮點運算在我的 MacOS 上來看，是雙準確度，非常精確的，真的越用越喜歡。
最近心血來潮，用它來測試一些地理計算的數學運算，例如地球上任兩點距離的計算，真的覺得非常方便。站在推廣 FORTH 語言的立場，來灌溉一下提供一下 FORTH 程式碼囉。
球面上兩點的距離：
這裏真的蠻有趣的，我們大家所熟悉的是歐幾里德幾何學。但因為我們生活在一個接近圓球的地球的表面上。所以我們的地理平面居然不是平的，而是彎曲有弧度的一個平面。
自然，我們要在這個彎曲的平面上標定座標，跟位置，及距離或是角度等。都無法使用歐幾里德幾何學的方式了。
而平面都彎曲掉了，當然所有的幾何學公式都完全不一樣囉。現在整個平面彎曲成一個圓球了，所以叫做球面幾何學囉。
所以，根據我們的生活直覺，我們一直以為我們所住的平面是平的，所以我們一直所使用來計算兩地的距離的歐幾里德距離公式是對的。但事實的真相卻全然不是那麼一回事，這真的是讓人蠻震驚且十分有趣的事情呢！
來比較一下，
在平面幾何，歐幾里德的平面空間下，兩點的距離 P(x1,y1) Q(x2,y2) 大家都再熟悉不過了. 就每個分量相減的平方和相加開根號 sqrt((x1-x2)^2+(y1-y2)^2)
它們最短的距離就是「直線距離」
在球面幾何，呵呵呵，這空間是彎曲的，所以兩點的距離 PQ，如下圖所示，因為在球的表面上，所以會是通過圓心的一的圓，被稱呼為大圓，在圓球表面的一個弧長囉。
(這裏所有的圖，懶得畫了，都借用維基百科啦)
緯度，經度：
在平面上，要標定平面上任一點的位置，我們可以用直角座標 x, y
而在彎曲的平面，一個球面上，如何標定上面一個點的位置呢？因為現在平面是個彎曲的一個球表面。所以球表面上面的點，無法直接用直角座標來描述了。(因為所有座標的垂直格線的間隔距離現在都不是像直角座標這樣永遠等距的囉，而是往這個球平面的南北極，間隔距離越來越靠近，然後最後相交於球平面上的南北極點。)
因爲球上不同位置下，垂直格線的間隔距離都不一樣，如果能改用角度來描述會比較容易些。所以現在是用以圓心出來，對應一個通過圓心參考軸的一個夾角(圓心角)來描述物體在這個球平面上的位置。
所以改用：
緯度：球平面的半圓(赤道)位置參考線，往上為正(北半球，或以N記號標示)，往下為負(南半球，或為正值但以S記號標示)，物體所在位置跟此半圓的夾角。所以範圍為 +/- 180度，或 0-180度 N 或 S。
經度：以某一條通過南北極的線當0度參考線(地球上此線通過英國倫敦天文台)，往右邊為正(或以E記號標示)，往左邊為負(或為正值但以W記號標示)，物體所在位置跟0度參考線的夾角。所以範圍為 +/- 360度，或 0-360度 E或 W。
就這樣，利用經度跟緯度這兩個參數，可以準確的標定一個彎曲球狀平面上的任何位置。就像是一般沒有彎曲的平面座標x 跟y
球面三角學：
在一般平面下，大家有歐幾里德幾何學。裡面有大家很熟悉畢氏定理，或兩點直線距離公式，三角形三內角和為180度，三角函數的正弦定理，餘弦定理...
但現在這個平面彎曲成一個球了，自然，這已經不是歐幾里德幾何學囉，所有的幾何定理會以不同的形貌出現。
這種在球面形狀的彎曲平面上的三角學，被稱之為球面三角學。
餘弦定理：
平面下的餘弦定理
球面上的餘弦定理
這裡要注意的，平面幾何下， a, b, c 一般是長度
在這裏球面幾何下， a, b, c 是球上對圓心的角度夾角(圓心角)，單位可是 rad 喲！
這裡球面上兩點的距離計算，答案就呼之欲出囉。
例如，假如我們想要計算 AB 這兩點球面上的距離，只要把 c 這個 AB 兩點對圓心的夾角(圓心角)先求出來，AB距離 = 球半徑 * 兩點夾角 = R * c
透過上面的餘弦定理,
cos c = cos a cos b + sin a sin b cos C
所以 c = acos( cos a cos b + sin a sin b cos C )
今天, 定義 C 的位置為球的北極，A位置的緯度為 LAT.a, 經度為 LONG.a ，B位置的緯度為 LAT.b, 經度為 LONG.b ，
所以
a = pi/2 - Lat.b ， b = pi/2 - Lat.a
所以
c = acos( cos(pi/2 - LAT.b) cos(pi/2 - LAT.a) + sin(pi/2 - LAT.b) sin(pi/2 - LAT.a) cos(LONG.b - LONG.a)
所以最終公式
c = acos( sin(LAT.b) sin(LAT.a) + cos(LAT.b)cos(LAT.a)cos(LONG.b - LONG.a)
求出 c 後，兩點距離 AB = R * c (R 為地球半徑)
半正矢公式 haversin
haversin(theta) = sin(theta)^2
上述球平面上兩點距離的公式，當兩點間的距離很短的時候，因為圓心角很小，此時餘弦函數接近於1 ，這時候餘弦函數的反函數用計算機來計算的時候，如果只是使用單精確的浮點數來計算的話會存在極大的捨入誤差而導致計算結果不夠精確。
這時候會使用半正矢函數和兩角和的餘弦函數展開式對公式展開變成如下的形式
(這裡懶得打了，直接貼維基百科的圖 XD )
這個也是在航海上運用廣泛的半正矢公式，歷史上會將距離和半正矢函數值的關係直接製成表格，方便使用。
所以，最終，我們的程式採用這個半正矢公式來進行計算
Place A: (緯度LAT.a , 經度LONG.a)
Place B: (緯度LAT.b , 經度LONG.b)
c = 2 asin( sqrt( sin( (LAT.b - LAT.a)/2 )^2 + cos(LAT.a) cos(LAT.b) sin( (LONG.b - LONG.a)/2 )^2 ) )
距離AB = R(地球半徑) * c
地球的平均半徑大概是 6371.009 km，
所以 R = 6371.009
FORTH 程式解說：
純用堆疊來運算這樣的數學方程式，當然是可以的。但是可讀性會變得很低。畢盡還是希望有一定的數學可讀性啦。
所以這裡經緯度用兩個變數來定義
位置1 的經緯度為
LAT1: 位置一的緯度，單位為 度 degree
LONG1: 位置一的經度，單位為 度 degree
LAT2: 位置二的緯度，單位為 度 degree
LONG2: 位置二的經度，單位為 度 degree
三角函數都是用 rad 來計算的，所以底下度轉rad, 或 rad轉度的轉換式不可少
就 rad = degree * pi / 180
或 degree = rad * 180 / pi
: Deg>Rad ( f: deg -- rad) pi f* 180E0 f/ ;
: Rad>Deg ( f: rad -- deg) 180E0 f* pi f/ ;
將經緯度由度數轉換成 rad 後存入變數中
: LAT1! ( f: deg --) Deg>Rad LAT1 f! ;
: LAT2! ( f: deg --) Deg>Rad LAT2 f! ;
: LONG1! ( f: deg --) Deg>Rad LONG1 f! ;
: LONG2! ( f: deg --) Deg>Rad LONG2 f! ;
公式的主角，計算兩位置間弧角的數值
: v-sin ( f: -- d-sig)
D-LAT 2E0 f/ fsin 2E0 f**
LAT1 f@ fcos LAT2 f@ fcos f* D-LONG 2E0 f/ fsin 2E0 f** f*
f+
fsqrt fasin
2E0 f*
;
最後這個弧角乘上地球的平均半徑，就得到我們所要的答案，兩位置間的距離囉！
: distance ( f: -- d ) Earth-R v-sin f* ;
實際使用常常需要，度，分，秒在那邊轉來轉去的。所以定義下面的指令來幫助計算
: deg ;
: min ( f: deg min -- deg' ) 60E0 f/ f+ ;
: sec ( f: deg sec -- deg' ) 3600E0 f/ f+ ;
: N ;
: S ( f: deg -- -deg ) fnegate ;
: E ;
: W ( f: deg -- -deg ) fnegate ;
例如
120E0 deg 46E0 min 32E0 sec N
這樣的指令就可以很方便的把 120°46'32" 往北的方向適當地轉換成角度 120.775555555556 方便後面的計算
結果驗證：
(1)
新竹城隍廟的緯度跟經度： 24.804331, 120.966188
清泉張學良文化紀念館的緯度跟經度：24.573802, 121.105514
根據 google Map, 兩者距離為 29.29 km
我們這個程式所算出來的是 29.2442204395348 km. 兩者相當的一致
(2)
從網路上隨便找的
地點A的緯度跟經度：50°03'59"N, 05°42'53"W
地點B的緯度跟經度：58°38'38"N, 03°04'12"W
這兩地的距離是 968.9 km
我們這個程式所算出來的是 968.854915365142 km. 兩者相當的一致
得證，所以我們程式碼應該正確無誤。
最後，原始程式碼列表
\
\ distance between two places on earth
\ 2020.6.18 Frank Lin
\
\ set word search order
FORTH only FORTH also
\ define new vocabulary
vocabulary GeoTools
\ define new words on this vocabulary
GeoTools definitions
6371.009E0 fconstant Earth-R
fvariable LAT1
fvariable LONG1
fvariable LAT2
fvariable LONG2
: Deg>Rad ( f: deg -- rad) pi f* 180E0 f/ ;
: Rad>Deg ( f: rad -- deg) 180E0 f* pi f/ ;
: LAT1! ( f: deg --) Deg>Rad LAT1 f! ;
: LAT2! ( f: deg --) Deg>Rad LAT2 f! ;
: LONG1! ( f: deg --) Deg>Rad LONG1 f! ;
: LONG2! ( f: deg --) Deg>Rad LONG2 f! ;
: D-LAT ( f: -- d-LAT) LAT1 f@ LAT2 f@ f- ;
: D-LONG ( f: -- d-LONG ) LONG1 f@ LONG2 f@ f- ;
: v-sin ( f: -- d-sig)
D-LAT 2E0 f/ fsin 2E0 f**
LAT1 f@ fcos LAT2 f@ fcos f* D-LONG 2E0 f/ fsin 2E0 f** f*
f+
fsqrt fasin
2E0 f*
;
: distance ( f: -- d ) Earth-R v-sin f* ;
\ extra converting words
: deg ;
: min ( f: deg min -- deg' ) 60E0 f/ f+ ;
: sec ( f: deg sec -- deg' ) 3600E0 f/ f+ ;
: N ;
: S ( f: deg -- -deg ) fnegate ;
: E ;
: W ( f: deg -- -deg ) fnegate ;
\ add to search order
GeoTools also
.( verify by ...) cr
.( The distance between 24.804331, 120.966188 and 24.573802, 121.105514 is 29.29 km )
cr
cr
.( The check result is ... )
24.804331E0 LAT1! 120.966188E0 LONG1!
24.573802E0 LAT2! 121.105514E0 LONG2!
distance f.
cr
cr
cr
.( The distance between 50 3 59N, 5 42 53W and 58 38 38N, 3 4 12W is 968.9 km )
cr
cr
.( The check result is ... )
50E0 deg 3E0 min 59E0 sec N LAT1! 5E0 deg 42E0 min 53E0 sec W LONG1!
58E0 deg 38E0 min 38E0 sec N LAT2! 3E0 deg 4E0 min 12E0 sec W LONG2!
distance f.
cr
cr
cr
XXX
留言列表
