球面三角法による距離計算

→測地系と投影

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

最新の50件