公共測量、土木設計、用地境界測量、土木建築用測量、土地調査、各種申請etc...株式会社 滝下測量設計事務所(京都府綾部市、福知山市)

'BLXY
Option Explicit

Const PI = 3.14159265358979

'楕円体の定数
Const a As Double = 6378137  '長半径
Const f = 1 / 298.257222101  '扁平率

'楕円体の諸公式
Const b  As Double = a * (1 - f)  '短半径
Const c  As Double = a * a / b    '極での曲率半径
Dim e As Double
Dim e2 As Double

'座標系の原点における縮尺係数
Const m0 = 0.9999

'子午線弧長計算用
Const A_ As Double = 1.00505250181309
Const B_ As Double = 0.005063108622224
Const C_ As Double = 0.000010627590263
Const D_ As Double = 0.000000020820379
Const E_ As Double = 0.000000000039324
Const F_ As Double = 0.000000000000071
Dim B1 As Double
Dim B2 As Double
Dim B3 As Double
Dim B4 As Double
Dim B5 As Double
Dim B6 As Double

'座標系原点
Dim LAT_0_list() As Variant '1系〜19系の原点北緯
Dim LNG_0_list() As Variant '1系〜19系の原点東経


'初期化
Sub BLXY_init()
    e = Sqr(2 * f - f ^ 2)            '第一離心率
    e2 = Sqr((a ^ 2 - b ^ 2) / b ^ 2) '第二離心率

'   W=  sqr(1.0-e*e * powf(sin(lat),2))
'   M=  c / powf(V,3)   '子午線曲率半径
'   R=  sqr(M*N)       '平均曲率半径

    B1 = a * (1 - e ^ 2) * A_
    B2 = a * (1 - e ^ 2) * (-B_ / 2)
    B3 = a * (1 - e ^ 2) * (C_ / 4)
    B4 = a * (1 - e ^ 2) * (-D_ / 6)
    B5 = a * (1 - e ^ 2) * (E_ / 8)
    B6 = a * (1 - e ^ 2) * (-F_ / 10)

'座標系の定義
    LAT_0_list = Array(0, 33, 33, 36, 33, 36, 36, 36, 36, 36, 40, 44, 44, _
                       44, 26, 26, 26, 26, 20, 26)
    LNG_0_list = Array(0, 129.3, 131, 132.1, 133.3, 134.2, 136, 137.1, _
                       138.3, 139.5, 140.5, 140.15, 142.15, 144.15, 142, _
                       127.3, 124, 131, 136, 154)

    'ラジアンにしておく
    Dim i As Integer
    For i = 1 To 19
        LAT_0_list(i) = DEG2RAD(DMS2DEG(LAT_0_list(i)))
        LNG_0_list(i) = DEG2RAD(DMS2DEG(LNG_0_list(i)))
    Next i

End Sub

'XYから経緯度を求める
'XY2BL(座標系、X、Y)
'戻り値(緯度,経度(DEG10進))
Function XY2BL(ByVal 座標系 As Integer, ByVal X As Double, ByVal Y As Double) As Double()

    If e = 0 Then Call BLXY_init
    If (座標系 < 1) Or (19 < 座標系) Then Exit Function

    Dim Φ1 As Double
    Φ1 = S2LAT(座標系, X) 'φ1:座標から基準子午線へ下らせる垂線の足の緯度

    Dim RetBuf(1) As Double
    Dim v As Double
    Dim N1 As Double
    Dim t1 As Double
    Dim eta1 As Double
    Dim ypm As Double
    Dim dLmd As Double
    Dim sign As Double
    v = Sqr(1 + e2 ^ 2 * (Cos(Φ1) ^ 2))
    N1 = c / v          '卯酉線曲率半径(引数φ1)
    t1 = Tan(Φ1)
    eta1 = Sqr(e2 ^ 2 * (Cos(Φ1) ^ 2))
    ypm = Abs(Y) / m0
    If Y < 0 Then
        sign = -1
    Else
        sign = 1
    End If

    Dim tmp(1 To 4) As Double

    '緯度
    tmp(1) = 1 / 2 * 1 / (N1 ^ 2) * t1 * (1 + eta1 ^ 2) * (ypm ^ 2)
    tmp(2) = 1 / 24 * 1 / (N1 ^ 4) * t1 * (5 + 3 * t1 ^ 2 + 6 * eta1 * _
             eta1 + 6 * t1 ^ 2 * eta1 ^ 2 - 3 * (eta1 ^ 4) - 9 * t1 * _
             t1 * (eta1 ^ 4)) * (ypm ^ 4)
    tmp(3) = 1 / 720 * 1 / (N1 ^ 6) * t1 * (61 + 90 * t1 ^ 2 + 45 * _
             (t1 ^ 4) + 107 * eta1 ^ 2 - 162 * t1 ^ 2 * eta1 ^ 2 - _
             45 * (t1 ^ 4) * eta1 ^ 2) * (ypm ^ 6)
    tmp(4) = 1 / 40320 * 1 / (N1 ^ 8) * t1 * (1385 + 3633 * t1 ^ 2 + 4095 * _
             (t1 ^ 4) + 1575 * (t1 ^ 6)) * (ypm ^ 8)

    RetBuf(0) = RAD2DEG(Φ1 - tmp(1) + tmp(2) - tmp(3) + tmp(4))

    '経度
    tmp(1) = 1 / (N1 * Cos(Φ1)) * ypm * sign
    tmp(2) = -1 / 6 * 1 / ((N1 ^ 3) * Cos(Φ1)) * (1 + 2 * t1 ^ 2 + eta1 * _
             eta1) * (ypm ^ 3) * sign
    tmp(3) = 1 / 120 * 1 / ((N1 ^ 5) * Cos(Φ1)) * (5 + 28 * t1 ^ 2 + 24 * _
            (t1 ^ 4) + 6 * eta1 ^ 2 + 8 * t1 ^ 2 * eta1 ^ 2) * (ypm ^ 5) * sign
    tmp(4) = -1 / 5040 * 1 / ((N1 ^ 7) * Cos(Φ1)) * (61 + 662 * t1 ^ 2 + _
             1320 * (t1 ^ 4) + 720 * (t1 ^ 6)) * (ypm ^ 7) * sign
    dLmd = tmp(1) + tmp(2) + tmp(3) + tmp(4)

    RetBuf(1) = RAD2DEG(LNG_0_list(座標系) + dLmd)

    XY2BL = RetBuf

End Function


'経緯度からXYを求める
'BL2XY(緯度,経度(DEG10進))
'戻り値(座標)
Function BL2XY(ByVal 座標系 As Integer, ByVal LAT As Double, ByVal LNG As Double) As Double()

    If e = 0 Then Call BLXY_init
    If (座標系 < 1) Or (19 < 座標系) Then Exit Function

    LAT = DEG2RAD(LAT)
    LNG = DEG2RAD(LNG)

    Dim dLmd As Double
    Dim sign As Double
    dLmd = LNG - LNG_0_list(座標系)  '経度の差(東方を正にとる
    If dLmd < 0 Then
        sign = -1
    Else
        sign = 1
    End If
    dLmd = Abs(dLmd)

    Dim RetBuf(1) As Double
    Dim eta As Double
    Dim t As Double
    Dim v As Double
    Dim N As Double
    Dim s0 As Double
    Dim S As Double
    eta = Sqr(e2 * e2 * (Cos(LAT) ^ 2))
    t = Tan(LAT)

    v = Sqr(1 + e2 * e2 * (Cos(LAT) ^ 2))
    N = c / v   '卯酉線曲率半径

    'X座標
    s0 = LAT2S(LAT_0_list(座標系)) '子午線弧長
    S = LAT2S(LAT)

    Dim tmp(1 To 4) As Double
    tmp(1) = (S - s0) + 0.5 * N * (Cos(LAT) ^ 2) * t * (dLmd ^ 2)
    tmp(2) = 1 / 24 * N * (Cos(LAT) ^ 4) * t * (5 - t ^ 2 + 9 * eta * _
             eta + 4 * (eta ^ 4)) * (dLmd ^ 4)
    tmp(3) = 1 / 720 * N * (Cos(LAT) ^ 6) * t * (-61 + 58 * t ^ 2 - _
             (t ^ 4) - 270 * eta * eta + 330 * t ^ 2 * eta * eta) * (dLmd ^ 6)
    tmp(4) = 1 / 40320 * N * (Cos(LAT) ^ 8) * t * (-1385 + 3111 * t * _
             t - 543 * (t ^ 4) + (t ^ 6)) * (dLmd ^ 8)
    RetBuf(0) = (tmp(1) + tmp(2) - tmp(3) - tmp(4)) * m0

    'Y座標
    tmp(1) = N * Cos(LAT) * dLmd * sign
    tmp(2) = 1 / 6 * N * (Cos(LAT) ^ 3) * (-1 + t ^ 2 - eta * eta) * (dLmd ^ 3) * sign
    tmp(3) = 1 / 120 * N * (Cos(LAT) ^ 5) * (-5 + 18 * t ^ 2 - _
             (t ^ 4) - 14 * eta * eta + 58 * t ^ 2 * eta * eta) * (dLmd ^ 5) * sign
    tmp(4) = 1 / 5040 * N * (Cos(LAT) ^ 7) * (-61 + 479 * t ^ 2 - _
             179 * (t ^ 4) + (t ^ 6)) * (dLmd ^ 7) * sign
    RetBuf(1) = (tmp(1) - tmp(2) - tmp(3) - tmp(4)) * m0

    BL2XY = RetBuf

End Function


'緯度を与えて赤道からの子午線弧長を求める
'引数=緯度(ラジアン)
Function LAT2S(ByVal Ido As Double) As Double
    LAT2S = B1 * Ido + B2 * Sin(2 * Ido) + B3 * Sin(4 * Ido) + _
            B4 * Sin(6 * Ido) + B5 * Sin(8 * Ido) + B6 * Sin(10 * Ido)
End Function

'赤道からの子午線弧長を与えて緯度を求める
Function S2LAT(ByVal 座標系 As Integer, ByVal X As Double) As Double

    If (座標系 < 1) Or (19 < 座標系) Then Exit Function

    Dim ansΦ_ As Double
    Dim ansΦ As Double
    Dim Φ As Double
    ansΦ_ = 0
    'φn:緯度(1回目は座標系の原点の緯度とし、2回目以降はφn+1=φnとする)
    Φ = LAT_0_list(座標系)

    Dim M As Double
    M = LAT2S(Φ) + (X / m0)
    Dim i As Integer
    Dim S As Double
    Dim tmp(1) As Double

    For i = 1 To 100
        S = LAT2S(Φ) 'φnの緯度に対する赤道からの子午線弧長

        tmp(0) = 2 * (S - M) * ((1 - e ^ 2 * (Sin(Φ) ^ 2) ^ 1.5))
        tmp(1) = 3 * e ^ 2 * (S - M) * Sin(Φ) * Cos(Φ) * (1 - e * _
                 e * (Sin(Φ) ^ 2) ^ 0.5) - 2 * a * (1 - e ^ 2)
        ansΦ = Φ + tmp(0) / tmp(1)

        '差がDEG 2"*10E-5以下になるまで反復計算
        If Abs(ansΦ - ansΦ_) < 0.0000000096962 Then Exit For
        Φ = ansΦ
        ansΦ_ = ansΦ
    Next i
    S2LAT = ansΦ

End Function

'DEGREE→RADIAN
Function DEG2RAD(ByVal DEG As Double) As Double
    DEG2RAD = (PI / 180) * DEG
End Function

'RADIAN→DEGREE
Function RAD2DEG(ByVal RAD As Double) As Double
    RAD2DEG = (180 / PI) * RAD
End Function

'Deg→DMMSS
Function DEG2DMS(ByVal Deg10 As Double) As Double

    Dim D As Double
    Dim M As Double
    Dim S As Double

    D = RDwn(Deg10)
    M = (RDwn(Frac(Deg10) * 60)) / 100
    S = Frac(Frac(Deg10) * 60) * 60 / 10000
    DEG2DMS = D + M + S

End Function
'DMS→DEG
Function DMS2DEG(ByVal DMMSS As Double) As Double

    Dim D As Double
    Dim M As Double
    Dim S As Double

    D = RDwn(DMMSS)
    M = RDwn(Frac(DMMSS) * 100) / 60
    S = Frac(Frac(DMMSS) * 100) * 100 / 3600
    DMS2DEG = D + M + S

End Function

'小数部取り出し
Function Frac(ByVal v As Double) As Double
    Frac = v - RDwn(v)
End Function

'切り捨て
Function RDwn(ByVal v As Double) As Double
    RDwn = CDbl(Application.RoundDown(v, 0))
End Function