球面三角法による距離計算 の変更点

• 追加された行はこの色です。
• 削除された行はこの色です。
• 球面三角法による距離計算 へ行く。
• 球面三角法による距離計算 の差分を削除

[[距離の計算]]
#author("2023-02-19T15:15:28+09:00","default:irrp","irrp")
→測地系と投影

 Private Const mATN1P45 = 0.0174533 'atn(1)/45
 Private Const ER = 6378.14
 Private Const PI = 3.14159265358979
 'LHT座標値例
 '北緯36度30分 → 365000
 '南緯40度45分 → - 407500
 '北緯136度30分 → 1365000
 '南緯140度45分 → -1407500
 Option Explicit
 
 '位置と距離と方位を指定して目的位置を求める
 Public Sub GetPos_Kyorihoui(ByVal Mypx&, ByVal Mypy&, ByVal Kyori#, ByVal Houi#, Topx&, Topy&)
 'Mypx　基準位置経度 Mypy 基準位置緯度　LHT座標
 'Kyori 距離　km　Houi 方位　度
 'Topx　目的位置経度 Topy 目的位置緯度　LHT座標
 '方位は0-359.99.. 北が0、東が90、南が180、西が270となる。
 
     Dim cb#, sb#, clc#, c#
     Dim xx#, ido#, kdo#
     ido = (900000 - Mypy) * (Atn(1) * 4) / 1800000
     c = Kyori / ER
     If c > (Atn(1) * 4) * 2 Then c = c - (Atn(1) * 4) * 2
     cb = Cos(c) * Cos(ido) + Sin(c) * Sin(ido) * Cos(Houi * mATN1P45)
     If cb = 1 Then
         sb = 0
     Else
         sb = Sqr(1 - cb * cb)
     End If
     If Sin(ido) * sb = 0 Then
         clc = 0
     Else
         clc = (Cos(c) - Cos(ido) * cb) / (Sin(ido) * sb)
         If clc > 1 Then clc = 1
         If clc < -1 Then clc = -1
     End If
     xx = cb
     If Sqr(-xx * xx + 1) = 0 Then
         ido = Atn(1) * 2
     Else
         ido = Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1)
     End If
     Topy = 900000 - ido * 1800000 / (Atn(1) * 4)
     xx = clc
     If xx = 1 Then
         kdo = 0
     ElseIf xx = -1 Then
         kdo = Atn(1) * 4
     Else
         If Sin(Houi * mATN1P45) >= 0 Then
             kdo = Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1)
         Else
             kdo = -(Atn(-xx / Sqr(-xx * xx + 1)) + 2 * Atn(1))
         End If
     End If
     If c > Atn(1) * 4 Then kdo = -kdo
     kdo = kdo * 1800000 / (Atn(1) * 4)
     Topx = Mypx + kdo
     If Topx > 1800000 Then Topx = Topx - 3600000
     If Topx < -1800000 Then Topx = Topx + 3600000
 End Sub
 
 '2点の位置を指定して距離と方位を求める。
 Public Sub GetKyoriHoui(ByVal Mypx&, ByVal Mypy&, ByVal Topx&, ByVal Topy&, Kyori#, Houi#)
 'Mypx　位置１経度 Mypy 位置２緯度　LHT座標
 'Topx　位置２経度 Topy 位置２緯度　LHT座標
 'Kyori 計算結果距離　km　Houi 計算結果方位　度
     Dim xx#, ram#, sig#, ssig#
     Dim mx#, my#, tx#, ty#
     On Error GoTo getkher
 
     '２点の位置を設定
     mx = Mypx
     tx = Topx
     my = Mypy
     ty = Topy
     If mx = tx And my = ty Then '010711
         Kyori = 0
         Houi = 0
         Exit Sub
     End If
 
     'ラジアン変換
     mx = mx * PI / 1800000
     my = my * PI / 1800000
     tx = tx * PI / 1800000
     ty = ty * PI / 1800000
 
     '距離を求める
     ram = tx - mx
     xx = Sin(my) * Sin(ty) + Cos(my) * Cos(ty) * Cos(ram)
     If xx = -1# Then
          sig = PI
     ElseIf xx = 1# Then
          sig = 0
     Else
          sig = Atn(-xx / Sqr(-xx * xx + 1)) + PI / 2
     End If
     Kyori = ER * sig
 
     '方位を求める
     ssig = Sin(sig)
     If ssig = 0 Then
          Houi = 0
     Else
          xx = Cos(ty) * Sin(ram) / ssig
          If -xx * xx + 1 < 0 Then
               If Sin(ram) > 0 Then
                    Houi = 90
               Else
                     Houi = 270
               End If
          ElseIf Sqr(-xx * xx + 1) = 0 Then
              If Sin(ram) > 0 Then
                    Houi = 90
              Else
                   Houi = 270
              End If
         Else
              Houi = Atn(xx / Sqr(-xx * xx + 1))
             If 0 > Cos(my) * Sin(ty) - Sin(my) * Cos(ty) * Cos(ram) Then
                  Houi = PI - Houi
             End If
              Houi = Houi * 180 / PI
             If Houi < 0 Then Houi = Houi + 360
          End If
     End If
 Exit Sub
 '参考
 'Arcsin(x) = Atn(x / Sqr(-x * x + 1))
 'Arccos(X) = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)
 getkher:
 'Errmes Err, "s_Getkyorihoui"
 Resume Next
 End Sub