高斯克吕格与经纬度坐标值转换代码

<FONT size=2>'Writen by Rodger Yuan 9 5 2006<BR>'参考文献<BR>'v0 0.1<BR>'用于在经纬度坐标和高斯克吕格坐标之间的转换。<BR>'高斯克吕格为一种投影，根据椭球体和基准面不同又有所区分，常用的北京54和西安80即<BR>'采用这种投影方式，投影后的坐标为平面坐标系，单位为米<BR>'现在参数的坐标系采用测绘坐标系，x为纵坐标，y为横坐标<BR>'返回参数为自定义类型，双精度点<BR>'调用转换函数前需要调用初始化过程进行初始化<BR>'-------------------------------------------------------------------------------<BR>'Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)<BR>'说明: 用于初始化转换参数<BR>'TuoqiuCanshu  枚举类型，提供北京54、西安80和WGS84三个椭球参数<BR>'Daihao  整型  为高斯克吕格投影六度分带带号，取值为1～60<BR>'-------------------------------------------------------------------------------<BR>'Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal k0 As Double)<BR>'说明: 用于进一步初始化转换参数（暂不提供）<BR>'dE  东偏移<BR>'dE  北偏移<BR>'k0  比例因子<BR>'-------------------------------------------------------------------------------<BR>'Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD<BR><BR></FONT>'************************************************************************<BR>'基本变量定义<BR>Dim a As Double '椭球体长半轴<BR>Dim b As Double '椭球体短半周<BR>Dim f As Double '扁率<BR>Dim e As Double '第一偏心率<BR>Dim eq As Double '第二偏心率<BR><BR>Dim dh As Integer '带号<BR>Dim FE As Double '东偏移<BR>Dim FN As Double '北偏移<BR>Dim L0 As Double '中央经度<BR>Dim k0 As Double '比例因子<BR><BR><BR>Const PI As Double = 3.14159265358979<BR>Public Enum Canshu<BR>    Beijing54 = 0<BR>    Xian80 = 1<BR>    WGS84 = 2<BR>End Enum<BR>Public Type PointD<BR>    X As Double<BR>    Y As Double<BR>End Type<BR><BR>Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)<BR>Select Case TuoqiuCanshu<BR>'Krassovsky （北京54采用） 6378245 6356863.0188<BR>'IAG 75（西安80采用） 6378140 6356755.2882<BR>'WGS 84 6378137 6356752.3142<BR><BR>Case 0: '北京五四<BR>    a = 6378245<BR>    b = 6356863.0188<BR>Case 1: '西安八零<BR>    a = 6378140<BR>    b = 6356755.2882<BR>Case 2: 'WGS84<BR>    a = 6378137<BR>    b = 6356752.3142<BR>End Select<BR>f = (a - b) / a<BR>e = Sqr(1 - (b / a) ^ 2)<BR>eq = Sqr((a / b) ^ 2 - 1)<BR><BR>If Daihao < 1 Or Daihao > 60 Then Exit Sub<BR>dh = Daihao<BR><BR>L0 = (6 * dh - 3) * PI / 180<BR>k0 = 1<BR>FE = 500000 + dh * 1000000<BR>FN = 0<BR><BR>End Sub<BR>Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal dk0 As Double)<BR><BR>End Sub<BR>Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD<BR>'给出经纬度坐标，转换为高克投影坐标<BR>Dim BY As Double<BR>Dim LX As Double<BR>Dim TC As Double<BR>Dim CC As Double<BR>Dim AC As Double<BR>Dim MC As Double<BR>Dim NC As Double<BR><BR>Dim rx As Double<BR>Dim ry As Double<BR>Dim resultP As PointD<BR>BY = W * PI / 180<BR>LX = J * PI / 180<BR>TC = Math.Tan(BY) ^ 2<BR>CC = eq ^ 2 * Cos(BY) ^ 2<BR>AC = (LX - L0) * Cos(BY)<BR>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 * BY) + (15 * e ^ 4 / 256 + 45 * e ^ 6 / 1024) * Sin(4 * BY) - (35 * e ^ 6 / 3072) * Sin(6 * BY))<BR>NC = a / Sqr(1 - e ^ 2 * (Sin(BY)) ^ 2)<BR>rx = k0 * (MC + NC * Tan(BY) * (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)<BR>ry = FE + k0 * NC * (AC + (1 - TC + CC) * AC ^ 3 / 6 + (5 - 18 * TC + TC ^ 2 + 14 * CC - 58 * TC * CC) * AC ^ 5 / 120)<BR>resultP.X = rx<BR>resultP.Y = ry<BR>JWgetGK = resultP<BR>End Function<BR><BR>Public Function GKgetJW(ByVal X As Double, ByVal Y As Double) As PointD<BR>'给出高克投影坐标，转换为经纬度坐标<BR>Dim BY As Double<BR>Dim LX As Double<BR>Dim e1 As Double<BR>Dim FI As Double<BR>Dim Mf As Double<BR>Dim Bf As Double<BR>Dim Tf As Double<BR>Dim Cf As Double<BR>Dim Nf As Double<BR>Dim Rf As Double<BR>Dim D As Double<BR><BR><BR>Dim RW As Double '纬度<BR>Dim RJ As Double '经度<BR>Dim resultP As PointD<BR>YE = Y<BR>XN = X<BR><BR><BR>e1 = (1 - b / a) / (1 + b / a)<BR>Mf = (XN - FN) / k0<BR>FI = Mf / (a * (1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256))<BR>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)<BR>Tf = Tan(Bf) ^ 2<BR>Cf = eq ^ 2 * Cos(Bf) ^ 2<BR>Nf = a / Sqr(1 - e ^ 2 * Sin(Bf) ^ 2)<BR>Rf = a * (1 - e ^ 2) / Sqr((1 - e ^ 2 * Sin(Bf) ^ 2) ^ 3)<BR>D = (YE - FE) / (k0 * Nf)<BR><BR><BR>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)<BR>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)<BR>resultP.X = RW * 180 / PI<BR>resultP.Y = RJ * 180 / PI<BR>GKgetJW = resultP<BR>End Function