Latitude and longitude and height projection conversion code (VB)

By: Program It seems that there is no easy-to-use library for coordinate conversion. In fact, the conversion between the two is commonly used. It is too difficult to write the definition of the coordinate system by yourself. I found the formula, I also wrote it by myself according to the formula. It is said that there is an error in this conversion. After comparison, it will be enough for dozens of meters.
'*************************************** *********************************
'Gauss kerlenge and Longitude and latitude Coordinate value conversion Code
'Writen by Rodger yuan 9 5 2006
' Reference Document
'V0 0.1
'Is used in Longitude and latitude The conversion between the coordinate and the Gaussian curlenge coordinate.
'The Gauss kerlenge is a projection based on Elliptical body And Datum Different Also, Beijing 54 and Xi'an 80 are commonly used.
'In this projection method, the coordinate after projection is the plane coordinate system, in meters
'The current parameter coordinate system uses the ing coordinate system, X is the vertical coordinate, and Y is the horizontal coordinate.
'The return parameter is of the custom type and has double precision.
'Before calling the conversion function, you must call the initialization process for initialization.
'-------------------------------------------------------------------------------
'Public sub Init (byval tuoqiucanshu as canshu, byval daihao as integer)
'Description: used to initialize conversion Parameters
'Tuoqiucanshu Enumeration type. It provides three elliptical parameters: Beijing 54, Xi'an 80, and WGS84.
'Daihao integer is a Gaussian lüge projection with a number of six degrees. The value ranges from 1 ~ 60
'-------------------------------------------------------------------------------
'Public sub initex (byval de as double, byval DN as double, byval K0 as double)
'Note: this parameter is used for further initialization of conversion parameters (not provided currently)
'De East offset
'De North offset
'K0 Proportional Factor
'-------------------------------------------------------------------------------
'Public function compute getgk (byval w as double, byval J as double) as pointd

'*************************************** *********************************
'Define basic variables
Dim A as double'Elliptical bodyLong Half Axis
Dim B as double'Elliptical bodyHalf a week
Dim F as double 'flat Rate
Dim e as double 'first eccentric Heart Rate
Dim eq as double 'second eccentric Heart Rate

Dim DH as integer 'with no.
Dim Fe as double 'East offset
Dim FN as double 'North offset
Dim l0 as double 'Central longitude
Dim K0 as double 'proportional factor

Const PI as double = 3.14159265358979
Public Enum canshu
Beijing54 = 0
Xian80 = 1
WGS84 = 2
End Enum
Public type pointd
X as double
Y as double
End type

Public sub Init (byval tuoqiucanshu as canshu, byval daihao as integer)
Select case tuoqiucanshu
'La sovsky (adopted in beijing 54) 6378245 6356863.0188
'Iag 75 (used in Xi'an 80) 6378140 6356755.2882
'Wgs 84 6378137 6356752.3142

Case 0: 'Beijing May 4
A = 6378245
B = 6356863.0188
Case 1: 'xi'an
A = 6378140
B = 6356755.2882
Case 2: 'wgs84
A = 6378137
B = 6356752.3142
End select
F = (a-B)/
E = sqr (1-(B/A) ^ 2)
Eq = sqr (A/B) ^ 2-1)

If daihao <1 or daihao> 60 then exit sub
DH = daihao

L0 = (6 * DH-3) x pi/180
K0 = 1
Fe = 500000 + DH * 1000000
Fn = 0

End sub
Public sub initex (byval de as double, byval DN as double, byval dk0 as double)

End sub
Public Function compute getgk (byval w as double, byval J as double) as pointd
'GivenLongitude and latitudeCoordinates, which are converted to high-grams projection coordinates.
Dim by as double
Dim lx as double
Dim TC as double
Dim cc as double
Dim AC as double
Dim MC as double
Dim NC as double

dim RX as double
dim ry as double
dim resultp as pointd
by = W * PI/180
Lx = J * PI/ 180
Tc = math. tan (by) ^ 2
cc = EQ ^ 2 * Cos (by) ^ 2
AC = (Lx-l0) * Cos ()
MC = A * (1-e ^ 2/4-3 * E ^ 4/64-5 * E ^ 6/256) * By-(3 * E ^ 2/8 + 3 * E ^ 4/32 + 45 * E ^ 6/1024) * sin (2 *) + (15 * E ^ 4/256 + 45 * E ^ 6/1024) * sin (4 * by)-(35 * E ^ 6/3072) * sin (6 * ))
NC = A/sqr (1-e ^ 2 * (sin (by) ^ 2)
RX = K0 * (MC + NC * Tan () * (AC ^ 2/2 + (5-TC + 9 * CC + 4 * Cc ^ 2) * AC ^ 4/24) + (61-58 * TC + t ^ 2 + 270 * Cc-330 * TC * CC) * AC ^ 6/720)
ry = Fe + K0 * NC * (AC + (1-TC + CC) * AC ^ 3/6 + (5-18 * TC + Tc ^ 2 + 14 * Cc-58 * TC * CC) * AC ^ 5/120)
resultp. X = RX
resultp. y = ry
required getgk = resultp
end function

Public Function gkgetjw (byval X as double, byval y as double) as pointd
'The height projection coordinate is given and convertedLongitude and latitudeCoordinates
Dim by as double
Dim lx as double
Dim E1 as double
Dim Fi as double
Dim MF as double
Dim BF as double
Dim TF as double
Dim cf as double
Dim NF as double
Dim RF as double
Dim D as double

Dim RW as double 'latitude
Dim RJ as double 'longitude
Dim resultp as pointd
Ye = y
Xn = x

E1 = (1-B/A)/(1 + B/)
Mf = (Xn-FN)/K0
FI = mf/(A * (1-e ^ 2/4-3 * E ^ 4/64-5 * E ^ 6/256 ))
BF = Fi + (3 * E1/2-27 * E1 ^ 3/32) * sin (2 * fi) + (21 * E1 ^ 2/16-55 * E1 ^ 4/32) * sin (4 * fi) + (151 * E1 ^ 3/96) * sin (6 * fi)
TF = tan (BF) ^ 2
Cf = EQ ^ 2 * Cos (BF) ^ 2
NF = A/sqr (1-e ^ 2 * sin (BF) ^ 2)
Rf = A * (1-e ^ 2)/sqr (1-e ^ 2 * sin (BF) ^ 2) ^ 3)
D = (Ye-Fe)/(k0 * NF)

RW = BF-(NF * Tan (BF)/Rf) * (d ^ 2/2-(5 + 3 * TF + CF-9 * tF * Cf) * d ^ 4/24 + (61 + 90 * TF + 45 * tF ^ 2) * d ^ 6/720)
RJ = l0 + 1/cos (BF) * (D-(1 + 2 * TF + CF) * d ^ 3/6 + (5 + 28 * TF + 6 * CF + 8 * tF * CF + 24 * tF ^ 2) * d ^ 5/120)
resultp. X = RW * 180/pI
resultp. y = RJ * 180/pI
gkgetjw = resultp
end function