高斯投影正反算的直接计算结果是什么

高斯投影坐标正反算

一、基本思想：

高斯投影正算公式就是由大地坐标（L ，B ）求解高斯平面坐标（x ，y ），而高斯投影反算公式则是由高斯平面坐标（x ，y)求解大地坐标（L ，B).

二、计算模型:

基本椭球参数：

椭球长半轴a

椭球扁率f

椭球短半轴:(1)b a f =-

椭球第一偏心率

第 1 页

高斯计_高精度_高效率_高品质

全面覆盖电离辐射， 非电离辐射及直流磁场的检测。低频电磁场还可在线监测， 多种功能供用户选配

点击立即咨询，了解更多详情

咨询

深圳市柯雷科技开发 广告

：e = 椭球第二偏心率

：e '=高斯投影正算公式：此公式换算的精度为0.001m

64256

442234

22)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+''++-''+''⋅''+

=ρηηρρ 52224255

32233

)5814185(cos 120)1(cos 6cos l t t t B N l t B

第 2 页

N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ

其中：角度都为弧度

B 为点的纬度,0l L L ''=-，L 为点的经度,0L 为中央子午线经度；

N 为子午圈曲率半径,1222

(1sin )N a e B -=-；

tan t B =； 222cos e B η'=

180

3600ρπ''=*

其中X 为子午线弧长：

2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ⎡⎤=--++-+⎢⎥⎣⎦

第 3 页

02468,,,,a a a a a 为基本常量，按如下公式计算：

2004682426844686868

83535281612815722321637816323216128m a m m m m m m a m m m a m m m m a m a ⎧=++++⎪⎪⎪=+++⎪⎪⎪=++⎨⎪⎪=+⎪⎪⎪=⎪⎩

02468,,,,m m m m m 为基本常量，按如下公式计算：

22222020426486379(1)5268

m a e m e m m e m m e m m e m =-====

高斯投影反算公式：此公式换算的精度为0.0001’'。

()()()()22222432465

第 4 页

3

2235

2422250

53922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f f

f f f f f f f

t t B B y t t y

M N M N t y t t y

M N y y l t N B N B y t t t N B L l L ηηηηη=-

+++--++=-+++++++=+

其中: 0L 为中央子午线经度。

第 5 页

f B 为底点纬度，也就是当x X =时的子午线弧长所对应的纬度。按照子午线弧长公式：68240sin 2sin 4sin 6sin82468

a a a a X a B B B B B =-+-+，迭代进行计算； 初始开始时设：10f B X a =

以后每次迭代按下式计算:

10

6824(())()sin 2sin 4sin 6sin82468i

i f f i

i i i i f

f f f f B X X F B a a a a a F B B B B B +=-=-+-+

重复迭代至1i

第 6 页

i

f f B B ε+-<；为止。

1222

(1sin )f f N a e B -=-

3

2222(1)(1sin )f f M a e e B -=-- tan f f t B =；

222cos f

f e B η'=

海福特椭球（1910)我国52年以前基准椭球 a=6378388m b=6356911。9461279m α=0.33670033670

克拉索夫斯基椭球（1940 Krassovsky)

第 7 页

北京54坐标系基准椭球 a=6378245m b=6356863。018773m α=0.33523298692

1975年I 。U 。G.G 推荐椭球（国际大地测量协会1975) 西安80坐标系基准椭球

a=6378140m b=6356755.2881575m α=0.0033528131778

WGS-84椭球（GPS 全球定位系统椭球、17届国际大地测量协会） WGS —84 GPS 基准椭球

a=6378137m b=6356752.3142451m α=0。00335281006247

三、程序代码函数:

/*＊＊**＊＊＊**＊＊高斯投影正算函数***＊＊**＊＊**＊*＊*

第 8 页

输入 : double a ，f 椭球参数，B,L 为大地坐标,L0为中央子午线的经度，单位为弧度，x,y 为高斯平面坐标,y 加上了500000常量

返回:none

＊＊＊**＊*＊*＊*＊＊＊＊*＊＊＊＊*＊*＊*＊*＊**＊＊*＊＊*＊＊*＊＊＊/

void gaosiforward (double a ,double f ，double B ，double L ，double L0，double ＆x ，double ＆y ） ｛

double b , c ，e1， e2； //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

double l ， W ,N ， M ， daihao ；//W 为常用辅助函数，N 为子午圈曲率半径，M 为卯酉圈曲率半径

第 9 页

double X //子午线弧长，高斯投影的坐标

double ruo ， ita , sb , cb ，t

double m [5］，n ［5]；

//计算一些基本常量

{

b =a ＊(1-f )；

e1=sqrt (a ＊a -b ＊b )/a

e2=sqrt （a *a -b ＊b ）/b

c =a ＊a /b ；

m ［0］=a *(1-e1＊e1） m ［1］=3＊（e1＊e1*m [0])/2。0

m[2］=5*(e1＊e1*m［1]）/4。0

第 10 页

m[3]=7*（e1*e1*m［2]）/6.0

m[4］=9*(e1*e1*m[3]）/8。0；

n[0]=m[0］+m［1]/2+3*m［2］/8+5＊m［3]/16+35*m[4]/128

n[1］=m［1］/2+m[2］/2+15*m[3]/32+7＊m[4]/16；

n[2］=m［2]/8+3＊m[3]/16+7*m［4]/32

n［3]=m［3]/32+m[4]/16；

n［4］=m［4]/128 /////by kjh 2014。5。22 把改成了

｝

//由纬度计算子午线弧长

第 11 页

{

X=n［0]*B—sin（B)*cos(B）＊（（n［1]-n[2］+n[3]）+（2*n［2]-（16*n[3］/3.0))＊sin(B)＊sin(B）+16*n[3]*pow(sin(B)，4)/3。0)

}

l=L—L0；//弧度

ita=e2＊cos（B)

sb=sin(B）；

cb=cos(B）；

W=sqrt（1-e1*e1*sb*sb）；

N=a/W；

第 12 页

t=tan(B)

ruo=（180/Pi)＊3600

x=(X+N＊sb*cb*l＊l/2+N＊sb＊cb*cb＊cb*（5—t*t+9＊ita*ita+4＊ita*ita＊ita＊ita）＊l*l*l ＊l/24+N＊sb*cb*cb*cb*cb*cb＊（61-58*t*t+t＊t*t*t)＊l*l＊l*l*l*l/720）

y=（N＊cb*l+N*cb*cb*cb＊（1-t*t+ita＊ita）＊l*l*l/6+N*cb*cb＊cb*cb*cb＊（5—18＊t＊t+t＊t*t＊t+14＊ita＊ita-58＊ita＊ita*t＊t)＊l*l*l＊l*l/120）；

y=y+500000

}

/＊＊＊＊＊＊*＊*＊＊*＊*高斯反算函数***＊****＊＊*****

第 13 页

输入： double a ，f 椭球参数， x，y为高斯平面坐标，L0为中央子午线的经度； B,L为大地坐标，单位为弧度

＊返回:none

＊*＊＊＊＊＊＊*＊＊＊＊**＊*＊*＊**＊＊＊*＊＊＊/

void gaosibackward(double a,double f,double x，double y,double L0,double＆B,double &L)

｛

double b, c,e1, e2//短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

double Bf,itaf,tf,Nf，Mf,Wf；

double l

第 14 页

double m[5］,n[5］

y=y-500000

//计算一些基本常量

｛

b=a＊（1—f）；

e1=sqrt(a*a—b*b）/a

e2=sqrt(a*a-b＊b)/b

c=a*a/b

m［0］=a＊（1-e1＊e1)

m[1]=3*（e1*e1*m［0］)/2.0；

m[2]=5*（e1＊e1*m［1］)/4。0；

第 15 页

m[3］=7*(e1＊e1*m[2］）/6。0；

m[4]=9*(e1＊e1＊m［3］)/8.0；

n［0]=m［0］+m［1]/2+3＊m［2]/8+5*m［3]/16+35*m［4］/128；

n[1]=m［1］/2+m[2]/2+15＊m［3]/32+7*m[4］/16

n[2］=m［2］/8+3*m［3］/16+7＊m[4]/32；

n[3]=m［3]/32+m[4］/16

n［4]=m［4]/128；

｝

//计算Bf

第 16 页

｛

double Bf1,Bfi0，Bfi1，FBfi；

Bf1=x/n[0]

Bfi0=Bf1

Bfi1=0

FBfi=0；

int num=0；

do

｛

num=0

FBfi=0.0-n[1］＊sin（2＊Bfi0）/2.0+n[2］*s

第 17 页

in(4*Bfi0）/4。0-n［3]*sin(6*Bfi0）/6。0；

Bfi1=（x—FBfi)/n［0]

if (fabs（Bfi1—Bfi0）>(Pi*pow(10.0,-8）/（36*18）)）

{

num=1；

Bfi0=Bfi1；

｝

｝while（num==1）

Bf=Bfi1；

}

tf=tan（Bf）；

第 18 页

Wf=sqrt（1-e1*e1＊sin（Bf）＊sin（Bf)）

Nf=a/Wf；

Mf=a＊（1—e1*e1)/(Wf＊Wf＊Wf)；

itaf=e2*cos（Bf)

B=Bf-tf＊y＊y/（2*Mf＊Nf）+tf*(5+3*tf＊tf+itaf＊itaf—9＊itaf*itaf*tf＊tf）*pow（y，4)/(24＊Mf＊pow(Nf，3））—tf＊（61+90*tf*tf+45*pow（tf,4）)*pow（y,6）/（720*Mf*pow（Nf，5)）

l=y/（Nf*cos(Bf））-(1+2*tf＊tf+itaf＊itaf)＊pow（y,3)/（6＊pow(Nf,3)＊cos（Bf))+(5+28＊

tf*tf+24＊pow(tf,4）+6*itaf*itaf+8*itaf＊itaf*tf*tf)*pow(y，5）/（120＊pow(Nf，5)*

第 19 页

cos(Bf）)

L=l+L0

｝

2014-5—22

’输入: double a ，f 椭球参数，B，L为大地坐标，L0为中央子午线的经度，单位为弧度，x，y为高斯平面坐标,y加上了常量

Private Function gaosiforward（ByVal a As Double, ByVal f As Double, ByVal B As Double，ByVal L As Double，ByVal L0 As Double）As Double（)

Dim x， y， xy(2) As Double

Dim bb, c, e1， e2 As Double’短半轴，极点

第 20 页

处的子午线曲率半径，第一偏心率，第二偏心率

Dim ll, W, N， M， daihao As Double’W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径

Dim xx As Double'子午线弧长，高斯投影的坐标

Dim ruo， ita， sb, cb, t As Double

Dim mm（5), nn（5) As Double

bb = a * （1 — f）

e1 = Math。Sqrt(a ＊ a - bb ＊ bb） / a

e2 = Math.Sqrt（a ＊ a — bb ＊ bb） / bb

c = a ＊ a / bb

第 21 页

mm（0) = a ＊（1 - e1 ＊ e1）

mm(1） = 3 ＊（e1 * e1 ＊ mm(0)） / 2.0

mm(2) = 5 * （e1 ＊ e1 ＊ mm（1)) / 4.0

mm（3) = 7 ＊（e1 ＊ e1 * mm(2)） / 6。0

mm（4） = 9 * （e1 ＊ e1 * mm（3）) / 8.0

nn(0) = mm(0） + mm（1) / 2 + 3 ＊ mm（2) / 8 + 5 * mm（3) / 16 + 35 ＊ mm（4） / 128nn(1) = mm(1) / 2 + mm(2） / 2 + 15 ＊ mm（3) / 32 + 7 ＊ mm（4） / 16

nn（2) = mm(2） / 8 + 3 ＊ mm（3) / 16 + 7 ＊ mm（4) / 32

nn(3) = mm（3) / 32 + mm(4） / 16

nn（4） = mm（4) / 128

第 22 页

xx = nn（0) ＊ B — Sin（B) * Cos（B）＊（（nn(1） - nn(2） + nn(3)) + (2 ＊ nn（2) - (16 ＊ nn(3） / 3。0）） * Sin（B）＊ Sin（B） + 16 * nn(3） * Pow(Sin（B), 4） / 3.0）

ll = L — L0 ’弧度

ita = e2 ＊ Cos（B)

sb = Sin（B)

cb = Cos(B)

W = Sqrt（1 - e1 ＊ e1 * sb * sb）

N = a / W

t = Tan(B）

ruo = (180 / PI） * 3600

第 23 页

x = (xx + N ＊ sb ＊ cb * ll ＊ ll / 2 + N ＊ sb ＊ cb ＊ cb ＊ cb * （5 - t ＊ t + 9 ＊ ita ＊ ita + 4 * ita ＊ ita * ita * ita) * ll ＊ ll * ll * ll / 24 + N * sb ＊ cb ＊ cb * cb ＊ cb ＊ cb * (61 — 58 ＊ t * t + t * t * t * t) * ll * ll ＊ ll * ll ＊ ll * ll / 720)

第 24 页

百度文库

搜索

百度文库10亿海量资料，查找管理一应俱全

打开APP

//高斯投影正、反算

//////6度带宽 54年北京坐标系

//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)

void GaussProjCal(double longitude, double latitude, double *X, double *Y)

{

int ProjNo=0int ZoneWide////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval

double a,f, e2,ee, NN, T,C,A, M, iPI

iPI = 0.0174532925199433////3.1415926535898/180.0

ZoneWide = 6////6度带宽

a=6378245.0f=1.0/298.3//54年北京坐标系参数

////a=6378140.0f=1/298.257//80年西安坐标系参数

ProjNo = (int)(longitude / ZoneWide)

longitude0 = ProjNo * ZoneWide + ZoneWide / 2

longitude0 = longitude0 * iPI

latitude0=0

longitude1 = longitude * iPI //经度转换为弧度

latitude1 = latitude * iPI //纬度转换为弧度

e2=2*f-f*f

ee=e2*(1.0-e2)

NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1))

T=tan(latitude1)*tan(latitude1)

C=ee*cos(latitude1)*cos(latitude1)

A=(longitude1-longitude0)*cos(latitude1)

M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2

*e2/1024)*sin(2*latitude1)

+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*sin(6*l

atitude1))

xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120)

yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24

+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720)

X0 = 1000000L*(ProjNo+1)+500000L

Y0 = 0

xval = xval+X0yval = yval+Y0

*X = xval

*Y = yval

}

//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)

void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) 字串9

{

int ProjNoint ZoneWide////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval

double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI

iPI = 0.0174532925199433////3.1415926535898/180.0

a = 6378245.0f = 1.0/298.3//54年北京坐标系参数

////a=6378140.0f=1/298.257//80年西安坐标系参数

ZoneWide = 6////6度带宽

ProjNo = (int)(X/1000000L) //查找带号

longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2

longitude0 = longitude0 * iPI //中央经线

X0 = ProjNo*1000000L+500000L

Y0 = 0

xval = X-X0yval = Y-Y0//带内大地坐标

e2 = 2*f-f*f

e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2))

ee = e2/(1-e2)

M = yval

u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256))

fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(

4*u)

+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u)

C = ee*cos(fai)*cos(fai)

T = tan(fai)*tan(fai)

NN = a/sqrt(1.0-e2*sin(fai)*sin(fai))字串1

R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin

(fai)*sin(fai)))

D = xval/NN

//计算经度(Longitude) 纬度(Latitude)

longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D

*D*D*D*D/120)/cos(fai)

latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24

+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720)

//转换为度 DD

*longitude = longitude1 / iPI

*latitude = latitude1 / iPI

}

NN卯酉圈曲率半径，测量学里面用N表示

M为子午线弧长，测量学里用大X表示 字串2

fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示 字串4

R为底点所对的曲率半径，测量学里用Nf表示

我有VB的，自己很多年前写的，一直用，但是正算－>反算－>正算后，Y坐标与原来的差了0.5－0.7mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！代码如下：

'高斯坐标正算

Private Sub DadiZs()

Dim t As Double, Itp As Double, X0 As Double, N As Double, L0 As Double

Dim V As Double, ll As Double, W As Double, M As Double

Lat = Radian(Lat)

Lon = Radian(Lon)

L0 = Radian(Lo)

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 6356863.01877305

ep = 0.006693421622966

ep1 = 0.006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7799 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6973 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

'X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7798 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6972 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

Else

a = 6378140 '75椭球参数

b = 6356755.28815753

ep = 0.006694384999588

ep1 = 0.006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 111133.0047 * (Lat * 180 / Pi) - (32009.8575 * Sin(Lat) + 133.9602 * (Sin(Lat)) ^ 3 + 0.6976 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

End If

ll = Lon - L0

t = Tan(Lat)

Itp = ep1 * Cos(Lat) ^ 2

W = Sqr(1 - ep * Sin(Lat) ^ 2)

V = Sqr(1 + ep1 * Cos(Lat) ^ 2)

M = c / V ^ 3

N = a / W

'x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp + 4 * Itp * Itp) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp - 330 * t ^ 2 * Itp) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320

x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp ^ 2 + 4 * Itp ^ 4) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp ^ 2 - 330 * t ^ 2 * Itp ^ 2) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320

y = N * Cos(Lat) * ll + N * (1 - t * t + Itp) * (Cos(Lat)) ^ 3 * ll ^ 3 / 6 + N * (5 - 18 * t * t + t ^ 4 + 14 * Itp - 58 * Itp * t * t) * (Cos(Lat)) ^ 5 * ll ^ 5 / 120 + N * (61 - 479 * t ^ 2 + 179 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 7 * ll ^ 7 / 5040

r = Sin(Lat) * ll + Sin(Lat) * (Cos(Lat)) ^ 2 * ll ^ 3 * (1 + 3 * Itp + 2 * Itp ^ 2) / 3 + Sin(Lat) * (Cos(Lat)) ^ 4 * ll ^ 5 * (2 - t * t) / 15

r = Degree(r)

y = y + 500000#

End Sub

'高斯反算

Private Sub DadiFs()

Dim t As Double, Itp As Double, X0 As Double, Bf As Double, N As Double

Dim v As Double, ll As Double, W As Double, M As Double, L0 As Double

L0 = Radian(Lo)

X0 = x * 0.000001

y = y - 500000#

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 6356863.01877305

ep = 0.006693421622966

ep1 = 0.006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 <3 Then

Bf = 9.04353301294 * X0 - 0.00000049604 * X0 ^ 2 - 0.00075310733 * X0 ^ 3 - 0.00000084307 * X0 ^ 4 - 0.00000426055 * X0 ^ 5 - 0.00000010148 * X0 ^ 6

ElseIf X0 <6 Then

Bf = 27.11115372595 + 9.02468257083 * (X0 - 3) - 0.00579740442 * (X0 - 3) ^ 2 - 0.00043532572 * (X0 - 3) ^ 3 + 0.00004857285 * (X0 - 3) ^ 4 + 0.00000215727 * (X0 - 3) ^ 5 - 0.00000019399 * (X0 - 3) ^ 6

End If

Else

a = 6378140 '75椭球参数

b = 6356755.28815753

ep = 0.006694384999588

ep1 = 0.006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 <3 Then

Bf = 9.04369066313 * X0 - 0.00000049618 * X0 ^ 2 - 0.00075325505 * X0 ^ 3 - 0.0000008433 * X0 ^ 4 - 0.00000426157 * X0 ^ 5 - 0.0000001015 * X0 ^ 6

ElseIf X0 <6 Then

Bf = 27.11162289465 + 9.02483657729 * (X0 - 3) - 0.00579850656 * (X0 - 3) ^ 2 - 0.00043540029 * (X0 - 3) ^ 3 + 0.00004858357 * (X0 - 3) ^ 4 + 0.00000215769 * (X0 - 3) ^ 5 - 0.00000019404 * (X0 - 3) ^ 6

End If

End If

Bf = Bf * Pi / 180#

t = Tan(Bf)

Itp = ep1 * Cos(Bf) ^ 2

W = Sqr(1 - ep * Sin(Bf) ^ 2)

v = Sqr(1 + ep1 * Cos(Bf) ^ 2)

M = c / v ^ 3

N = a / W

Lat = Bf - 0.5 * v ^ 2 * t * ((y / N) ^ 2 - (5 + 3 * t * t + Itp - 9 * Itp * t * t) * (y / N) ^ 4 / 12 + (61 + 90 * t * t + 45 * t ^ 4) * (y / N) ^ 6 / 360)

ll = ((y / N) - (1 + 2 * t * t + Itp) * (y / N) ^ 3 / 6 + (5 + 28 * t * t + 24 * t ^ 4 + 6 * Itp + 8 * Itp * t * t) * (y / N) ^ 5 / 120) / Cos(Bf)

r = y * t / N - y ^ 3 * t * (1 + t * t - Itp) / (3 * N ^ 3) + y ^ 5 * t * (2 + 5 * t * t + 3 * t ^ 4) / (15 * N ^ 5)

Lat = Degree(Lat)

Lon = Degree(L0 + ll)

r = Degree(r)

End Sub

有了正反算，换带也就完成了！

用到的子程序：

Public Const Pi = 3.14159265358979, p = 206264.806

Public Cktq As String

'角度化弧度

Public Function Radian(a As Double) As Double

Dim Ro As Double

Dim c As Double

Dim Fs As Double

Dim Ib As Integer

Dim Ic As Integer

If a <0 Then a = -a: t = 1

Ro = Pi / 180#

Ib = Int(a)

c = (a - Ib) * 100#

Ic = Int(c + 0.000000000001)

Fs = (c - Ic) * 100#

If t = 1 Then Radian = -(Ib + Ic / 60# + Fs / 3600#) * Ro Else Radian = (Ib + Ic / 60# + Fs / 3600#) * Ro

End Function

'弧度化角度

Public Function Degree(a As Double) As Double

Dim Bo As Double

Dim Fs As Double

Dim Im As Integer

Dim Id As Integer

If a <0 Then a = -a: t = 1

Bo = a

Call DMS(Bo, Id, Im, Fs)

If t = 1 Then Degree = -(Id + Im / 100# + Fs / 10000#) Else Degree = Id + Im / 100# + Fs / 10000#

End Function

Public Sub DMS(a As Double, Id As Integer, Im As Integer, Fs As Double)

Dim Bo As Double

Dim c As Double

c = a

c = 180# / Pi * c

Id = Int(c)

Bo = (c - Id) * 60

Im = Int(Bo)

Fs = (Bo - Im) * 60

End Sub

'取位计算

Public Function Qw(a As Double, Ws As Integer) As Double

Qw = Int(a * 10 ^ Ws + 0.5) / 10 ^ Ws

End Function

---