El Rincón de Tolgalen

Hola,

os dejo unas funciones en Visual Basic 6 para cálculos geodésicos que espero os sean útiles.

Para las conversiones UTM-Decimal se emplean las ecuaciones de Coticchia-Surace. Teneis más información sobre su implementación en la web de Gabriel Ortiz

Option Explicit

'This software is free software; you can redistribute it and/or modify it under the terms of
'the GNU General Public License as published by the Free Software Foundation; either version 2
'of the License, or (at your option) any later version.
'
'This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
'without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
'See the GNU General Public License for more details.

'Copyright (C) 2006 Carlos Francisco Andión López candion@gmail.com

'**********************************************************************************
' Módulo de diversas funciones geodésicas
'**********************************************************************************

Const pi As Double = 3.14159265358979
Const e As Double = 2.71828182845905

'**********************************************************************************
' Funcion que convierte grados decimales en sexagesimales
'**********************************************************************************
Function DecToSex(ByVal Decimalgrad As Double) As String

    Dim grados As Double
    Dim minutos As Double
    Dim segundos As Double
    Dim tmp As Double

    tmp = Decimalgrad

    If Decimalgrad < 0 Then
        Decimalgrad = -Decimalgrad
    End If

    grados = Int(Decimalgrad)
    minutos = Int((Decimalgrad - grados) * 60)
    segundos = Left(((Decimalgrad - grados) * 60 - minutos) * 60, 10)

      If tmp < 0 Then
        DecToSex = "-" & grados & "º " & minutos & "´" & segundos
            Else
            DecToSex = grados & "º " & minutos & "´" & segundos
    End If

End Function

'**********************************************************************************
' Funcion que convierte grados sexagesimales en decimales
'**********************************************************************************
Function SexToDec(ByVal grados As Integer, ByVal minutos As Double, ByVal segundos As Double) As Double

    SexToDec = grados + minutos / 60 + segundos / 3600

End Function

'**********************************************************************************
' Funcion que convierte coordenadas UTM a Decimal
' Emplea las ecuaciones de Coticchia-Surace empleando el elipsoide de Hayford o internacional 1924
' Devuelve un array con las coordenadas (latitud,longitud)
'**********************************************************************************

Public Function Utm2Dec(ByVal X As Double, ByVal y As Double, ByVal zone As String, ByVal Hemisf As String) As Double()

    Dim SMayor As Double
    Dim SMenor As Double
    Dim Ex As Double
    Dim Ex2 As Double
    Dim Ex2c As Double
    Dim Rc As Double
    Dim Apl As Double
    Dim Lambda As Double
    Dim Phi As Double
    Dim Ny As Double
    Dim A As Double
    Dim A1 As Double
    Dim A2 As Double
    Dim J2 As Double
    Dim J4 As Double
    Dim J6 As Double
    Dim Alfa As Double
    Dim Beta As Double
    Dim Gamma As Double
    Dim b As Double
    Dim bphi As Double
    Dim Sigma As Double
    Dim Xi As Double
    Dim Eta As Double
    Dim Tau As Double
    Dim Sen_h_Xi As Double
    Dim DeltaLambda As Double
    Dim coord() As Double

    ReDim coord(2)

    ' Datos elipsoide hayford
    SMayor = 6378388#
    SMenor = 6356911.94613
    Ex = Sqr(SMayor ^ 2 - SMenor ^ 2) / SMayor
    Ex2 = Sqr(SMayor ^ 2 - SMenor ^ 2) / SMenor
    Ex2c = Ex2 ^ 2
    Rc = (SMayor ^ 2) / SMenor
    Apl = (SMayor - SMenor) / SMayor

    Lambda = zone * 6 - 183
    X = X - 500000

    Select Case Hemisf
        Case "N", "n"
            y = y

        Case "S", "s"
            y = y - 10000000
    End Select

     ' aplico las ecuaciones de Coticchia-Surace
    Phi = y / (6366197.724 * 0.9996)
    Ny = (Rc / (1 + Ex2c * (Cos(Phi)) ^ 2) ^ (1 / 2)) * 0.9996
    A = X / Ny
    A1 = Sin(2 * Phi)
    A2 = A1 * (Cos(Phi)) ^ 2
    J2 = Phi + (A1 / 2)
    J4 = (3 * J2 + A2) / 4
    J6 = (5 * J4 + A2 * (Cos(Phi)) ^ 2) / 3
    Alfa = (3 / 4) * Ex2c
    Beta = (5 / 3) * Alfa ^ 2
    Gamma = (35 / 27) * Alfa ^ 3
    bphi = 0.9996 * Rc * (Phi - Alfa * J2 + Beta * J4 - Gamma * J6)
    b = (y - bphi) / Ny
    Sigma = ((Ex2c * A ^ 2) / 2) * (Cos(Phi)) ^ 2
    Xi = A * (1 - Sigma / 3)
    Eta = b * (1 - Sigma) + Phi
    Sen_h_Xi = ((e ^ Xi) - (e ^ (-Xi))) / 2
    DeltaLambda = ArcTan(Sen_h_Xi / Cos(Eta))
    Tau = ArcTan(Cos(DeltaLambda) * Tan(Eta))

    ' calculo la longitud
    coord(1) = (DeltaLambda / pi * 180) + Lambda

    'calculo la latitud
    coord(0) = (Phi + (1 + Ex2c * (Cos(Phi)) ^ 2 - (3 / 2) * Ex2c * Sin(Phi) * Cos(Phi) * (Tau - Phi)) * (Tau - Phi)) / pi * 180

    Utm2Dec = coord

End Function

'**********************************************************************************
' Funcion que convierte coordenadas  Decimales en UTM
' Emplea las ecuaciones de Coticchia-Surace empleando el elipsoide de Hayford o internacional 1924
' Devuelve un array con las coordenadas (latitud,longitud)
'**********************************************************************************

Public Function Dec2Utm(ByVal X As Double, ByVal y As Double, ByVal Hemisf As String) As Double()

    Dim SMayor As Double
    Dim SMenor As Double
    Dim Ex As Double
    Dim Ex2 As Double
    Dim Ex2c As Double
    Dim Rc As Double
    Dim Apl As Double
    Dim Lambda As Double
    Dim Phi As Double
    Dim Ny As Double
    Dim A As Double
    Dim A1 As Double
    Dim A2 As Double
    Dim J2 As Double
    Dim J4 As Double
    Dim J6 As Double
    Dim Alfa As Double
    Dim Beta As Double
    Dim Gamma As Double
    Dim bphi As Double
    Dim Sigma As Double
    Dim Xi As Double
    Dim Eta As Double
    Dim DeltaLambda As Double
    Dim coord() As Double
    Dim huso As Integer

    ReDim coord(2)

    ' Datos elipsoide hayford
    SMayor = 6378388#
    SMenor = 6356911.94613
    Ex = Sqr(SMayor ^ 2 - SMenor ^ 2) / SMayor
    Ex2 = Sqr(SMayor ^ 2 - SMenor ^ 2) / SMenor
    Ex2c = Ex2 ^ 2
    Rc = (SMayor ^ 2) / SMenor
    Apl = (SMayor - SMenor) / SMayor

    ' calculo el meridiano central del huso y la distancia entre la posicion y el meridiano central
    huso = Int(y / 6 + 31)
    Lambda = huso * 6 - 183

    ' paso a radianes las coordenadas
    X = X * pi / 180
    y = y * pi / 180

    DeltaLambda = y - (Lambda * pi / 180)

    ' aplico las ecuaciones de Coticchia-Surace
    A = Cos(X) * Sin(DeltaLambda)
    Xi = (1 / 2) * log((1 + A) / (1 - A))
    Eta = ArcTan(Tan(X) / Cos(DeltaLambda)) - X
    Ny = (Rc / (1 + Ex2c * (Cos(X)) ^ 2) ^ (1 / 2)) * 0.9996
    Sigma = (Ex2c / 2) * Xi ^ 2 * Cos(X) ^ 2
    A1 = Sin(2 * X)
    A2 = A1 * Cos(X) ^ 2
    J2 = X + A1 / 2
    J4 = (3 * J2 + A2) / 4
    J6 = (5 * J4 + A2 * Cos(X) ^ 2) / 3
    Alfa = 3 / 4 * Ex2c
    Beta = 5 / 3 * Alfa ^ 2
    Gamma = 35 / 27 * Alfa ^ 3
    bphi = 0.9996 * Rc * (X - Alfa * J2 + Beta * J4 - Gamma * J6)

    ' calculo la longitud Y
        Select Case Hemisf
            Case "N", "n"
                coord(1) = Eta * Ny * (1 + Sigma) + bphi
            Case "S", "s"
                coord(1) = Eta * Ny * (1 + Sigma) + bphi + 10000000
        End Select

    'calculo la latitud X
    coord(0) = Xi * Ny * (1 + Sigma / 3) + 500000

    Dec2Utm = coord

End Function

'**********************************************************************************
' Calcula la Cotangente inversa
'**********************************************************************************
Public Function ArcTan(ByVal X As Double) As Double
    ArcTan = Atn(X)
End Function

'**********************************************************************************
'Calcula el seno inverso. El angulo debe estar en radianes
'**********************************************************************************
Public Function ArcSin(ByVal Angle As Double) As Double
  If Abs(Angle) <> 1# Then
    ArcSin = Atn(Angle / Sqr(-Angle * Angle + 1#))
  Else
    ArcSin = IIf(Angle = 1#, Atn(1#) * 2#, Atn(1#) * 6#)
  End If

End Function

'**********************************************************************************
'Calcula el coseno inverso. El angulo debe estar en radianes
'**********************************************************************************
Public Function ArcCos(ByVal Angle As Double) As Double
  If Abs(Angle) <> 1# Then
      ArcCos = Atn(-Angle / Sqr(-Angle * Angle + 1#)) + 2# * Atn(1#)
   Else
      ArcCos = IIf(Angle = 1#, 0#, Atn(1#) * 4#)
   End If

End Function