坐标转换和单位转换的函数源代码

大家看看，请指正！

'函数：经纬度坐标转换成大地坐标
'输入参数：ZB,纬度,经度,坐标系统（如：80坐标）,中央经线
'输出参数：大地X,大地Y

Sub BLtoXY(Zb As ZbxPa, B As Double, L As Double, tyzbx As String, CenterLine As Single)
    Dim ZS As Double, Zc As Double
    Dim Ztan As Double, lc As Double, Nf As Double, n As Double
    Dim Q As Double, x As Double, y As Double
    Dim Jch
    Dim RstX As Double, RstY As Double
    Call SleParamer(tyzbx, Zb)
    ZS = Sin(B * W)
    Zc = Cos(B * W)
    Ztan = Tan(B * W)
    lc = (L - CenterLine) * 3600
    Nf = Zb.e1 * Zc * Zc
    n = Zb.ea / Sqr(1 - Zb.E * ZS * ZS)
    Q = Zb.tC0 * B * W + Zb.tC1 * Sin(2 * B * W) + Zb.tC2 * Sin(4 * B * W) + Zb.tC3 * Sin(6 * B * W) + Zb.tC4 * Sin(8 * B * W) + Zb.tC5 * Sin(10 * B * W)
    x = Q + n * Ztan * Zc ^ 2 * lc ^ 2 / (2 * PI * PI)
    x = x + n * Ztan * (5 - Ztan ^ 2 + 9 * Nf ^ 2 + 4 * Nf ^ 4) * Zc ^ 4 * lc ^ 4 / (24 * PI ^ 4)
    RstX = x + n * Ztan * (61 - 58 * Ztan ^ 2 + Ztan ^ 4 + 270 * Nf ^ 2 - 330 * Nf ^ 2 * Ztan ^ 2) * Zc ^ 6 * lc ^ 6 / (720 * PI ^ 6)
    RstX = RstX + (RstX < 1000000) * 0.000005
    y = n * Zc * (lc / PI)
    y = y + n * (1 - Ztan ^ 2 + Nf) * Zc ^ 3 * lc ^ 3 / (6 * PI ^ 3)
    RstY = y + n * (5 - 18 * Ztan ^ 2 + Ztan ^ 4 + 14 * Nf - 58 * Nf * Ztan ^ 2) * Zc ^ 5 * lc ^ 5 / (120 * PI ^ 5)
    RstY = RstY + Sgn(RstY) * 0.000005
    RstY = RstY '+ 500000  '大庆测试通过
    gRstX = RstY + CF + DaiHao: gRstY = RstX
    'Debug.Print RstX, RstY
End Sub

'2.==============================================================================================
'函数：大地坐标转换成经纬度坐标
'输入参数：坐标系(54|80)、自然值、度带、中央经线、带号Boolean、Zb、X、Y
'输出参数：经度、纬度

Sub XYtoBL(zbx As Integer, ByVal btt As String, ld As String, ByVal mCenterline As Single, ByVal DaiNO As Boolean, Zb As ZbxPa, x As Double, y As Double)
    '判断带号
    Dim Dh As Integer
    If DaiNO Then
        '有带号
        'YZbxtxt = "1980 北京坐标系"
        'WZbxtxt = "1954 北京坐标系"
        y = Val(Mid$(CStr(y), 3))
    Else
        y = y
    End If
    If btt = "50自然值" Then
        y = y - CF
    ElseIf btt = "D50自然值" Then
        y = Val(Mid$(Str(y), 4)) - CF
        Dh = Val(Mid(Str(y), 2, 2))
        If y < 0 And Abs(y) > 1000000 Then
            y = Val(Mid$(Trim(Str(y)), 3)) - CF
            '===================================
            Dh = Val(Mid(Trim(Str(y)), 1, 2))
            '===================================
        End If
        If ld = "6度带" Then
            mCenterline = OriDMStoDeg(Dh * 6 - 3)
        Else
            mCenterline = OriDMStoDeg(Dh * 3)
        End If
    Else
        y = y
    End If
    '----------------------------------------------------------------------------------
    Dim Zcos As Double, ZSin As Double
    Dim E As Double, z As Double, c As Double, t As Double, V As Double, U As Double
    Dim H As Double, n As Double, M As Double, c1 As Double, c2 As Double, c3 As Double
    Dim BF As Double, Q As Double
    Dim I As Integer
    Dim Tf, Nf, Mf, a1#, b1#, B, CF1#
    Dim ka As Double, kb As Double, kc As Double, ka1 As Double, kb1 As Double, kc1 As Double
    Dim kt1 As Double, kt2 As Double, kt3 As Double
    Call SleParamerGS(zbx, Zb)
    BF = x / (Zb.tA * W)
    Zcos = Cos(BF * W)
    ZSin = Sin(BF * W)
    ka = Zb.tB / (2 * Zb.tA)
    ka1 = ka
    kb = -Zb.tC / (4 * Zb.tA)
    kb1 = kb
    kc = Zb.tD / (6 * Zb.tA)
    kc1 = kc
    For I = 1 To 3
      kt1 = ka + ka * kb1 - 3 * ka * ka1 ^ 2 / 2 - 2 * kb * ka1
      kt2 = kb + ka * ka1
      kt3 = kc + ka * kb1 + ka * ka1 ^ 2 / 2 + 2 * kb * ka1
      ka1 = kt1
      kb1 = kt2
      kc1 = kt3
    Next
    Zb.tK1 = 2 * ka1 + 4 * kb1 + 6 * kc1
    Zb.tK2 = 8 * kb1 + 32 * kc1
    Zb.tK3 = 32 * kc1
    BF = BF * W + Zcos * (Zb.tK1 * ZSin - Zb.tK2 * ZSin ^ 3 + Zb.tK3 * ZSin ^ 5 + Zb.tK4 * ZSin ^ 7)
    Tf = Tan(BF)  ' 计算 B
    Nf = Zb.ea / Sqr(1 - Zb.E * Sin(BF) * Sin(BF))
    Mf = Zb.e1 * Cos(BF) ^ 2
    a1# = Tf * (1 + Mf) * (y / Nf) ^ 2 / 2
    b1# = Tf * (1 + Mf) * (5 + 3 * Tf ^ 2 + Mf - 9 * Mf * Tf ^ 2) * y ^ 4 / (24 * Nf ^ 4) 'Tf * (5 + 3 * Tf * Tf + 6 * Mf * Mf * (1 - Tf * Tf)) * y ^ 4 / (24 * Nf ^ 4)
    c1# = Tf * (1 + Mf) * (61 + 90 * Tf ^ 2 + 45 * Tf ^ 4) * y ^ 6 / (720 * Nf ^ 6) 'Tf * (61 + 45 * Tf * Tf * (2 + Tf * Tf)) * y ^ 6 / (720 * Nf ^ 6)
    B = BF - a1# + b1# - c1#
    gLat = B / W
    CF1# = Cos(BF) ' 计算 L
    a1# = y / (Nf * CF1#)
    b1# = (1 + 2 * Tf * Tf + Mf) * (y / Nf) ^ 3 / (6 * CF1#)
    c1# = (5 + 28 * Tf ^ 2 + 24 * Tf ^ 4 + Mf * (6 + 8 * Tf ^ 2)) * (y / Nf) ^ 5 / (120 * CF1#)
    gLgt = mCenterline + (a1# - b1# + c1#) / W
End Sub


'经纬度==〉高斯坐标(坐标转换的中间过程函数)
Public Sub SleParamer(tyzbx As String, Zb As ZbxPa)
    If tyzbx = "54" Then
        '54坐标
        Zb.ea = a01 '6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b01 '6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E01 '0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e01 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    Else
        '80坐标
        Zb.ea = a02 ' 6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b02 ' 6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E02 ' 0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e02 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    End If
    Zb.tA = Zb.ea * (1 - Zb.E) * (1 + 3 * Zb.E / 4 + 45 * Zb.E ^ 2 / 64 + 175 * Zb.E ^ 3 / 256 + 11025 * Zb.E ^ 4 / 16384)
    Zb.tB = Zb.ea * (1 - Zb.E) * (3 * Zb.E / 4 + 15 * Zb.E ^ 2 / 16 + 525 * Zb.E ^ 3 / 512 + 2205 * Zb.E ^ 4 / 2048)
    Zb.tC = Zb.ea * (1 - Zb.E) * (15 * Zb.E ^ 2 / 64 + 105 * Zb.E ^ 3 / 256 + 2205 * Zb.E ^ 4 / 4096 + 10395 * Zb.E ^ 5 / 16384)
    Zb.tD = Zb.ea * (1 - Zb.E) * (35 * Zb.E ^ 3 / 512 + 315 * Zb.E ^ 4 / 2048 + 31185 * Zb.E ^ 5 / 131072)
    Zb.tE = Zb.ea * (1 - Zb.E) * (315 * Zb.E ^ 4 / 16384 + 3465 * Zb.E ^ 5 / 65536)
    Zb.Tf = Zb.ea * (1 - Zb.E) * (693 * Zb.E ^ 5 / 131072)
    Zb.tC0 = Zb.tA
    Zb.tC1 = -Zb.tB / 2
    Zb.tC2 = Zb.tC / 4
    Zb.tC3 = -Zb.tD / 6
    Zb.tC4 = Zb.tE / 8
    Zb.tC5 = -Zb.Tf / 10
    Zb.cA = 1 + Zb.E / 2 + 3 * Zb.E ^ 2 / 4 + 5 * Zb.E ^ 3
    Zb.cB = Zb.E / 6 + 3 * Zb.E ^ 2 / 16 + 3 * Zb.E ^ 3 / 16
    Zb.cc = 3 * Zb.E ^ 2 / 80 + Zb.E ^ 3 / 16
    Zb.cD = (Zb.E ^ 3 / 112)
End Sub

'大地--〉经纬度(坐标转换的中间过程函数)
Public Sub SleParamerGS(cor_ID As Integer, Zb As ZbxPa)
    If cor_ID = 54 Then
        '54坐标
        Zb.ea = a01 '6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b01 '6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E01 '0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e01 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    Else
        '80坐标
        Zb.ea = a02 ' 6378245 ' Csb("a")   ' 6378140#
        Zb.eb = b02 ' 6356863.01877 'Csb("b")      ' 6356755.28815753
        Zb.E = str_E02 ' 0.00669342162297 'Csb("e2")    '0.006694384999588    '   e
        Zb.e1 = e02 ' 0.00673852541468 'Csb("e3") ' 0.006739501819473   '   e'
    End If
    Zb.tA = Zb.ea * (1 - Zb.E) * (1 + 3 * Zb.E / 4 + 45 * Zb.E ^ 2 / 64 + 175 * Zb.E ^ 3 / 256 + 11025 * Zb.E ^ 4 / 16384) ' + 43659 * Zb.E ^ 5 / 65536)
    Zb.tB = Zb.ea * (1 - Zb.E) * (3 * Zb.E / 4 + 15 * Zb.E ^ 2 / 16 + 525 * Zb.E ^ 3 / 512 + 2205 * Zb.E ^ 4 / 2048) '+ 72765 * Zb.E ^ 5 / 65536)
    Zb.tC = Zb.ea * (1 - Zb.E) * (15 * Zb.E ^ 2 / 64 + 105 * Zb.E ^ 3 / 256 + 2205 * Zb.E ^ 4 / 4096 + 10395 * Zb.E ^ 5 / 16384)
    Zb.tD = Zb.ea * (1 - Zb.E) * (35 * Zb.E ^ 3 / 512 + 315 * Zb.E ^ 4 / 2048 + 31185 * Zb.E ^ 5 / 131072)
    Zb.tE = Zb.ea * (1 - Zb.E) * (315 * Zb.E ^ 4 / 16384 + 3465 * Zb.E ^ 5 / 65536)
    Zb.Tf = Zb.ea * (1 - Zb.E) * (693 * Zb.E ^ 5 / 131072)
    Zb.tC0 = Zb.tA
    Zb.tC1 = -Zb.tB / 2
    Zb.tC2 = Zb.tC / 4
    Zb.tC3 = -Zb.tD / 6
    Zb.tC4 = Zb.tE / 8
    Zb.tC5 = -Zb.Tf / 10
    Zb.cA = 1 + Zb.E / 2 + 3 * Zb.E ^ 2 / 4 + 5 * Zb.E ^ 3
    Zb.cB = Zb.E / 6 + 3 * Zb.E ^ 2 / 16 + 3 * Zb.E ^ 3 / 16
    Zb.cc = 3 * Zb.E ^ 2 / 80 + Zb.E ^ 3 / 16
    Zb.cD = (Zb.E ^ 3 / 112)
End Sub

'==========================================
'函数功能：度分秒换算为度
'输入参数：度分秒(45.3030表示45度30分30秒)
'输出参数：度(45.5083333333333)

'==========================================
Public Function OriDMStoDeg(DMS As Double) As Double
    Dim Strdms$, d#, M#, Feng#
    Dim n As Integer
    Strdms$ = Str(DMS)
    n = InStr(Strdms$, ".")
    If n = 0 Then
        OriDMStoDeg = DMS
        Exit Function
    Else
        Strdms$ = Str(DMS) + "0000"
    End If
    d# = Val(Mid$(Strdms$, 1, n - 1))
    Feng# = (Val(Mid$(Strdms$, n + 1, 2))) / 60
    M# = (Val(Mid$(Strdms$, n + 3, 2) + "." + Mid$(Strdms$, n + 5))) / 3600
    OriDMStoDeg = (d# + Feng# + M#)
End Function

'==========================================
'函数功能：度换算为度分秒
'输入参数：度(45.5083333333333)
'输出参数：度分秒(45.3030表示45度30分30秒)

'==========================================
Public Function OriDegToDms(Deg As Double) As Double
    Dim TempDeg#, Zd#, AllF, Zf, strZf, p
     TempDeg# = Deg
     If InStr(Str(Deg), ".") = 0 Then
        OriDegToDms = Str(Deg) + ".0000"
        Exit Function
     End If
     Zd# = Int(TempDeg#)
     AllF = (TempDeg# - Zd#) * 60
     Zf = Int(AllF)
     If Zf < 10 Then
        strZf = "0" + Trim$(Str(Zf))
     Else
        strZf = Trim$(Str(Zf))
        If Val(strZf) = 60 Then Zd# = Zd# + 1
     End If
    Dim Allm, Strm
     Allm = (AllF - Zf) * 60
     If Allm < 10 Then
        Strm = "0" + Trim$(Str(Int(Allm)))
     Else
        If Allm >= 60 Then
            strZf = Str(Val(strZf) + 1)
            Strm = Str(Allm - 60)
            If Val(strZf) = 60 Then
                Zd# = Zd# + 1
                strZf = "00"
            End If
        Else
            Strm = Trim$(Str(Allm))
        End If
        p = InStr(Strm, ".")
        If p <> 0 Then
            Strm = Left(Strm, p - 1) + Mid(Strm, p + 1)
        End If
     End If
     OriDegToDms = CDbl(Str(Zd#) + "." + strZf + Strm)
End Function


<img src="images/post/smile/dvbbs/em09.gif" />

感谢楼上的兄弟的大公无私！

另外我想给斑竹提个建议：能不能专门开辟一个MO代码共享专栏，从实现MO最基本的功能开始，然后由大家把自己所想到的功能加上去，就像是Linux的发展一样！而且一定要声明代码是完全共享的!不知斑竹一下如何！
       也可以先在论坛中讨论一下嘛！

不过我觉得mo的功能并不是特别多,基本功能好多人都实现了,而且做得比较好,我很赞成open sourece,还希望大家多提点意见，做到更好的学习!多谢!<img src="images/post/smile/dvbbs/em05.gif" />

好，我看不懂，但是我支持

以下是引用cdhello在2003-8-26 16:36:15的发言：
谢谢GIS,有空一起做项目.

兄弟多联系，呵呵

先谢过，在慢慢看。

zbxpa的结构是什么?能不能说一下，不好意思，本人有些小笨，不能像使用联*电脑的人士一样，整天联想我的计算机到底怎么了。

斑竹，能不能告诉我在您的程序中bltoxy中的w代表什么含义？谢谢！

<P>总统，我爱你，就像猫咪爱老鼠。</P>
<P>虽然没看这些代码，我想肯定会给我不少帮助。</P>