[[距離の計算]] #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 位置1経度 Mypy 位置2緯度 LHT座標 'Topx 位置2経度 Topy 位置2緯度 LHT座標 'Kyori 計算結果距離 km Houi 計算結果方位 度 Dim xx#, ram#, sig#, ssig# Dim mx#, my#, tx#, ty# On Error GoTo getkher '2点の位置を設定 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