高斯坐标正反算代码

高斯投影正反算1. 什么是高斯投影高斯投影是一种常用的地图投影方法，它将地球表面的经纬度坐标转换为平面坐标，常用于地理信息系统（GIS）和测绘工程中。

高斯投影分为正算和反算两个过程。

•正算：将经纬度坐标转换为平面坐标。

•反算：将平面坐标转换为经纬度坐标。

2. 高斯投影正算2.1 原理高斯投影正算的原理是根据椭球体上某一点处的曲率半径、子午线弧长和东西方向上的距离，计算该点在平面上的x、y坐标。

2.2 具体步骤高斯投影正算的具体步骤如下：1.根据给定的椭球体参数（长半轴a、短半轴b），计算椭球体第一偏心率e。

2.根据给定的中央子午线经度λ0，计算λ - λ0 的差值Δλ。

3.计算曲率半径N和子午线弧长A0。

4.根据给定的纬度φ和经度λ，计算Δφ和Δλ。

5.计算子午线弧长A1、A2、A3和A4。

6.计算平面坐标x和y。

2.3 Python代码实现下面是使用Python实现高斯投影正算的示例代码：import math# 输入参数a = 6378137.0 # 长半轴b = 6356752.314245 # 短半轴e = math.sqrt(1 - (b/a)**2) # 第一偏心率λ0 = math.radians(120) # 中央子午线经度，单位为弧度# 输入经纬度坐标φ = math.radians(30) # 纬度，单位为弧度λ = math.radians(121) # 经度，单位为弧度# 计算ΔλΔλ = λ - λ0# 计算曲率半径N和子午线弧长A0N = a / math.sqrt(1 - e**2 * math.sin(φ)**2)A0 = a * (1 - e**2) / (1 - e**2 * math.sin(φ)**2)**1.5# 计算Δφ和ΔλΔφ = φ - φ0# 计算子午线弧长A1、A2、A3和A4A1 = A0 + N * math.tan(φ) / 2 * Δλ**2 * math.cos(φ)A2 = A0 + N * math.tan(φ) / 24 * (5 - math.tan(φ)**2 + 9 * e2 * math.cos(φ)**2 + 4 * e2**2 * math.cos(φ)**4) * Δλ**4 * math.cos(φ)A3 = A0 + N * math.tan(φ) / 720 * (61 - 58 * math.tan(φ)**2 + math.tan(φ)** 4) * Δλ**6 * math.cos(φ)A4 = A0 + N * math.tan(φ) / 40320 * (1385 - 3111*math.tan(φ)**2 + 543*math.tan(φ)**4 - math.tan(φ)**6) \* Δλ**8 \* math.cos(phi)# 计算平面坐标x和yx = A1 + A2 + A3 + A4y = N / math.cos(phi) \(Δλ - Δλ**3/6*(1+math.tan(phi))**2/N/A0^2 \+ Δλ^5/120*(5+28*math.tan(phi)^2+24*math.tan(phi)^4)*N/A0^4/N/A0^3)# 输出结果print("平面坐标(x, y):", x, y)3. 高斯投影反算3.1 原理高斯投影反算的原理是根据平面坐标和中央子午线经度，计算对应的经纬度坐标。

C#⾼斯正算⾼斯反算⾼斯换带等⾸先，你要确定椭球参数：C#代码1. a = 6378140; //西安80椭球 IGA752. e2 = 0.006694384999588;3. m0 = a * (1 - e2);4. m2 = 3.0 / 2 * e2 * m0;5. m4 = 5.0 / 4 * e2 * m2;6. m6 =7.0 / 6 * e2 * m4;7. m8 = 9.0 / 8 * e2 * m6;8. a0 = m0 + m2 / 2 + (3.0 / 8.0) * m4 + (5.0 / 16.0) * m6 + (35.0 / 128.0) * m8;9. a2 = m2 / 2 + m4 / 2 + 15.0 / 32 * m6 + 7.0 / 16 * m8;10. a4 = m4 / 8 + 3.0 / 16 * m6 + 7.0 / 32 * m8;11. a6 = m6 / 32 + m8 / 16;12. a8 = m8 / 128;13. xx = 0;14. yy = 0;15. _x = 0;16. _y = 0;17. BB = 0;18. LL = 0;下⾯才开始正题：⾼斯正算：把经纬度坐标转换为平⾯坐标C#代码1. void GaussPositive(double B, double L, double L0)2. {3. double X, t, N, h2, l, m, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;5. Bdu = (int)B;6. Bfen = (int)(B * 100) % 100;7. Bmiao = (B - Bdu - Bfen * 0.01) * 10000.0;8. B = Bdu * PI / 180.0 + (Bfen / 60.0) * PI / 180.0 + Bmiao / 3600.0 * PI / 180.0;9. Ldu = (int)L;10. Lfen = (int)(L * 100) % 100;11. Lmiao = (L - Ldu - Lfen * 0.01) * 10000.0;12. L = Ldu * PI / 180.0 + (Lfen / 60.0) * PI / 180 + Lmiao / 3600.0 * PI / 180.0;13. l = L - L0 * PI / 180;14. X = a0 * B - Math.Sin(B) * Math.Cos(B) * ((a2 - a4 + a6) + (2 * a4 - 16.0 / 3.0 * a6) * Math.Sin(B) * Math.Sin(B) + 16.0 / 3.0 * a6 * Math.Pow(Math.Sin(B), 4)) + a8 / 8.0 * Math.Sin(8 * B);15. t = Math.Tan(B);16. h2 = e2 / (1 - e2) * Math.Cos(B) * Math.Cos(B);17. N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));18. m = Math.Cos(B) * l;19. xx = X + N * t * ((0.5 + (1.0 / 24.0 * (5 - t * t + 9 * h2 + 4 * h2 * h2) + 1.0 / 720.0 * (61 - 58 * t * t + Math.Pow(t, 4)) * m * m) * m * m) * m * m);20. yy = N * ((1 + (1.0 / 6.0 * (1 - t * t + h2) + 1.0 / 120.0 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * h2 - 58 * h2 * t * t) * m * m) * m * m) * m);21. yy = yy + 500000;22. }⾼斯反算：把平⾯坐标转换为经纬度坐标：C#代码1. void GaussNegative(double x, double y, double L0)2.3. double Bf, Vf, l, tf, hf2, Nf, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;4. int Bdu, Bfen, Ldu, Lfen;5. y = y - 500000;6. Bf = hcfansuan(x);7. Vf = Math.Sqrt(1 + e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf));8. tf = Math.Tan(Bf);9. hf2 = e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf);10. Nf = a / Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));11. BB = (Bf - 0.5 * Vf * Vf * tf * (Math.Pow(y / Nf, 2) - 1.0 / 12 * (5 + 3 * tf * tf + hf2 - 9 * hf2 * tf * tf) * Math.Pow(y / Nf, 4) + 1.0 / 360 * (61 + 90 * tf * tf + 45 * tf * tf) * Math.Pow(y / Nf, 6))) * 180 / PI;12. Bdu = (int)BB;13. Bfen = (int)((BB - Bdu) * 60);14. Bmiao = ((BB - Bdu) * 60 - Bfen) * 60;15. BB = Bdu + 0.01 * Bfen + 0.0001 * Bmiao;16. l = 1.0 / Math.Cos(Bf) * (y / Nf - 1.0 / 6.0 * (1 + 2 * tf * tf + hf2) * Math.Pow(y / Nf, 3) + 1.0 / 120.0 * (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * hf2 + 8 * hf2 * tf * tf) * Math.Pow(y / Nf, 5)) * 180.0 / PI;17. LL = L0 + l;18. Ldu = (int)LL;19. Lfen = (int)((LL - Ldu) * 60);20. Lmiao = ((LL - Ldu) * 60 - Lfen) * 60;21. LL = Ldu + 0.01 * Lfen + 0.0001 * Lmiao;⾥⾯涉及到的弧长反算：C#代码1. double hcfansuan(double pX)2. {3. double Bf0 = pX / a0;4. double Bf1, Bf2;5. Bf1 = Bf0;6. Bf2 = (pX - F(Bf1)) / a0;7. while ((Bf2 - Bf1) > 1.0E-11)8. {9. Bf1 = Bf2;10. Bf2 = (pX - F(Bf1)) / a0;11. }12. return Bf1;13. }⾼斯换带就⽐较简单了：C#代码1. void GaussZone(double x, double y, double L0, double L0new)2. {3. GaussNegative(x, y, L0);4. GaussPositive(BB, LL, L0new);5. }。

高斯五点法坐标正反算程序（For5800）一、主程序(正反算程序)程序名:GSZFS1.7→Dim Z: “ZS=1,FS=2”? T:T=1=>Goto 1:T=2=>Goto 2 (T=1正算，T=2反算)2.“1.ZX,2.YX”?E （左线输1，右线输2，如有其它线路可再添加）3.Lbl 1:?Z:?D:Prog“GS-SJK”:Prog"GS-ZS":Rec(Abs(D),F+90D/Abs(D+0.00001)):X+I→X:Y+J→Y :Cls:“X=”:Locate 3,1,X:“Y=”:Locate 3,2,Y:“FWJ=”:F ►DMS◢4.Goto 15.Lbl 2:”X”?C:”Y”?T:?Z: Prog“GS-FS”:Cls:“Z=”:Locate 3,1,Z:“D=”:Locate 3,2,J◢6.Goto 2二、正算子程序，程序名: GS-ZS1.0.5(R-P) ÷H→S:180÷π→Z[7]:Abs(Z-O)→W2.0.1184634425→A3.0.2393143352→B4.0.2844444444→N5.0.046910077→K6.0.2307653449→L7.0.5→M8.G+Z[7]QKW(P+KWS)→Z[1]9.G+Z[7]QLW(P+LWS)→Z[2]10.G+Z[7]QMW(P+MWS)→Z[3]11.G+Z[7]Q(1-L)W(P+(1-L)WS)→Z[4]12.G+Z[7]Q(1-K)W(P+(1-K)WS)→Z[5]13.U+W(Acos(Z[1])+Bcos(Z[2])+Ncos(Z[3])+Bcos(Z[4])+Acos(Z[5]))→X14.V+W(Asin(Z[1])+Bsin(Z[2])+Nsin(Z[3])+Bsin(Z[4])+Asin(Z[5]))→Y15.G+Z[7]QW(P+WS)→F三、反算子程序，程序名: GS-FS1.Lbl 1:Prog”GS-SJK”: Prog"GS-ZS"2.Pol(C-X+0.00001,T-Y+0.00001):Rec(I,J-F):Z+I→Z3.Abs(I) >0.001=>Goto 1四、数据库子程序，程序名: GS-SJK1.E=1=>Goto 1:E=2=>Goto 22.Lbl 1:I f Z≥25900AndZ≤26615.555:Then25900→O:11587.421→U:1847.983→V:101°09’23.1”→G:715.555→H:1÷10^(45)→P:1÷10^(45)→R:0→Q:IfEnd3.……4.Lbl 2:……（同Lbl 1……）子程序“SJK”字母说明：O-线元起点桩号U-起点X坐标V-起点Y坐标G－线元起点桩号切线方位角H－线元长度P－线元起点曲率R－线元终点曲率Q－线元偏转方向，以线路的前进方向区分左右，左偏Q=-1，右偏Q=1，当线元为直线时Q=0。

高斯投影坐标正反算编程报告1.编程思想高斯投影坐标正反算的编程需要涉及大量公式。

为了使程序更清晰，各模块的数据重用性更强，本文采用结构化编程思想。

程序由四大块组成。

大地测量家庭作业。

Cpp文件用于存储main（）函数，是整个程序的入口。

尝试通过结构化编程简化main（）函数。

myfunction.h和myfunction.cpp用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。

正蒜。

H和zhengsuan CPP用于存储zhengsuan类。

在正算类中，声明了用于高斯投影坐标正演计算的所有变量。

成员变量的初始化和正向计算在类的构造函数中执行。

通过get函数获得相应的阳性结果。

fansuan.h和fansuan.cpp用于存放fansuan类，类似于zhengsuan类，fansuan类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。

通过get函数获得相应的反算结果。

2.计算模型高斯投影正算公式十、十、nn232244？？？？sinbcosb？Lsimbcosb（5？t？9？4？）l2？？？224??? 4n5246？？sinbcosb（61？58t？t）l720？？？6y？？nn3223？？cosb？Lcosb（1？t？？）L6.3n5242225？？cosb（5？18吨？14？58吨）l120？？？五tf2mfnf5f高斯投影反算公式B男朋友？？tfy？2tf24mfn3f？5.3t2f224？？2f？9? ftfy？720mfn46y61？90t2？45泰夫？？yy32l？？1.2t2f？？f3nfcosbf6nfcosbf？？？y524222?5?28t?24t?6??8?fffftf5120nfcosbf?3.程序框图一开始输入b，l求定带号n，中央纬度l0，纬度差l正算按照实用公式计算x,y换算为国家统一坐标x，y输出x,y输入国家统一坐标x,y由y取定带号n，并换算出x，y求出中央经线l0反算按照实用公式计算b，ll=l0+l求出大地经度l输出b，l结束4.计算结果25.附录：程序代码/////主函数入口大地测量家庭作业。

#include <stdio.h>#include <math.h>int main (void){ int o;printf ("\n============================欢迎使用高斯投影计算程序============================\n");printf(" 说明\n");printf("==此程序可实现克拉索夫斯基椭球和1975国际椭球上进行高斯三度带和六度带正算和反算==\n");printf ("\n");printf ("=========================请选择正反算。

1为正算，2为反算。

=======================\n");scanf ("%d",&o);if (o==1){int L0,i,e,a,b,c,f,g,L1,u,w;float d,h;double pi = 3.141592653;double P = 206264.806247;double x,y,l,N,B,a0,a4,a6,a3,a5,B2,C,L;printf ("=======================请输入分带方法。

1为3度带，2为6度带=======================\n");scanf ("%d",&u);printf ("=====================请输入经纬度。

例：1，00，00；31，00，00====================\n");scanf ("%d,%d,%d;%d,%d,%d",&a,&b,&c,&e,&f,&g);if (u==1){L1 = (e/3)*3;};if (u==2){L1 = (((e+f/60+g/3600)/6)+1)*6-3;};d = a+(b/60)+((c/60)/60);B = (d * pi)/180;C = cos(B)*cos(B);L = e*3600+f*60+g;B2 = a*3600+b*60+c;L0 = L1*3600;l = (L-L0)/P;printf ("================请选择椭球。

高斯投影正、反算代码//高斯投影正、反算//////6度带宽 54年北京坐标系//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y){int ProjNo=0; int 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.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=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/30 72)*sin(6*latitude1));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+X0; yval = yval+Y0;*X = xval;*Y = yval;}//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude){int ProjNo; int 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.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=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-X0; yval = 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)*si n(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));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;}。

高斯投影坐标反算c语言代码#include<stdio.h>#include<math.h>#include<conio.h>main(){printf("#####################################################\n");printf("# 角度输入说明：如26°12′45.2″输入为26,12,45.2 #\n");printf("#####################################################\n");double x,y;int j,L0;printf("请输入高斯投影坐标（自然坐标），中间用逗号隔开：\n");scanf("%lf,%lf",&x,&y);//自然坐标输入printf("请输入中央子午线L0:\n");scanf("%d,%d,%lf",&L0); //中央子午线输入printf("请选择参考椭球：1.北京1954参考椭球。

\n 2.西安1980参考椭球。

\n");printf("选择的参考椭球为：");scanf("%d",&j);//选择椭球参数if(j==1){long double Bf0=0.157046064172*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005051773759*sin(Bf0)-0.000029837302*pow(sin(Bf0),3)+0.00000023818 9*pow(sin(Bf0),5));long double t=tan(Bf);long double m=0.00673852541468*pow(cos(Bf),2);long double V=1+m;long double N=6378245.000/sqrt(1-0.00669342162297*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d3)*60);long double m3=((L-d3)*60-f3)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //北京1954参考椭球}if(j==2){long double Bf0=0.157048687473*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005052505593*sin(Bf0)-0.000029847335*pow(sin(Bf0),3)+0.0000002416* pow(sin(Bf0),5)+0.0000000022*pow(sin(Bf0),7));long double t=tan(Bf);long double m=0.006739501819473*pow(cos(Bf),2);long double V=1+m;long double N=6378140.000/sqrt(1-0.006694384999588*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d2)*60);long double m3=((L-d2)*60-f2)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //西安1980参考椭球}getch();}。

CCoordinate CEarthPlaneDialog::EarthToPlane(const CEarth & earth){double B = earth.m_dB.m_realRadians;double L = earth.m_dL.m_realRadians;L = SetRad_PIPI(L);double L0 = SetL0(L);double X, N, t, t2, m, m2, ng2;double sinB, cosB;X = A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 * B);sinB = sin(B);cosB = cos(B);t = tan(B);t2 = t * t;N = a / sqrt(1 - e2 * sinB * sinB);m = cosB * (L - L0);m2 = m * m;ng2 = cosB * cosB * e2 / (1 - e2);double x = X + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2);double y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0));y += 500000;CCoordinate coor(x,y);return coor;}//将大地坐标系转换为平面直角坐标系CEarth CEarthPlaneDialog::PlaneToEarth(const CCoordinate & coor, UINT Sign){double x = coor.m_dX;double y = coor.m_dY;double L0 = (6.0 * Sign - 3.0) / 180 * PI;double sinB, cosB, t, t2, N ,ng2, V, yN;double preB0, B0;double eta;y -= 500000;B0 = x / A1;do{preB0 = B0;B0 = B0 * PI / 180.0;B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1;eta = fabs(B0 - preB0);}while(eta > 0.0000000000000001);B0 = B0 * PI / 180.0;sinB = sin(B0);cosB = cos(B0);t = tan(B0);t2 = t * t;N = a / sqrt(1 - e2 * sinB * sinB);ng2 = cosB * cosB * e2 / (1 - e2);V = sqrt(1 + ng2);yN = y / N;double B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2;double L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB;CEarth earth(B,L,0,'r');return earth;}//将平面直角坐标系转换为大地坐标系void CEarthPlaneDialog::OnRADIOPlane1954(){// TODO: Add your control notification handler code herea = 6378245;f = 298.3;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111134.8611;A2 = -16036.4803;A3 = 16.8281;A4 = -0.0220;m_iReferenceFrame = 0;}void CEarthPlaneDialog::OnRADIOPlane1980(){// TODO: Add your control notification handler code herea = 6378140;f = 298.257;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111133.0047;A2 = -16038.5282;A3 = 16.8326;A4 = -0.0220;m_iReferenceFrame = 1;}void CEarthPlaneDialog::OnRADIOPlaneWGS84(){// TODO: Add your control notification handler code herea = 6378137;f = 298.257223563;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111133.0047;A2 = -16038.5282;A3 = 16.8326;A4 = -0.0220;m_iReferenceFrame = 2;}double CEarthPlaneDialog::SetL0(double L){L = L / PI * 180;if( L > 0 ){for(m_iCincture = 1 ; !(( L >= (6 * m_iCincture - 6 )) & ( L < ( 6 * m_iCincture ))); m_iCincture++){}L0 = 6 * m_iCincture - 3;}else if( L < 0 ){L = -L;for(m_iCincture = 1 ; !(( L > (6 * m_iCincture - 6 )) & ( L <= ( 6 * m_iCincture ))); m_iCincture++){}m_iCincture = 61 - m_iCincture;L0 = m_iCincture * 6 - 363;}else{L0 = 3;m_iCincture = 1;}double temp = ((double)L0) / 180 * PI;return temp;}double CEarthPlaneDialog::SetRad_PIPI(double Rad) {if(Rad > 0){double i = floor(Rad / PI / 2);Rad -= 2 * PI * i;}if(Rad < 0){Rad = -Rad;double i = floor(Rad / PI / 2);Rad -= 2 * PI * i;Rad = -Rad + PI * 2;}if(Rad > PI)Rad -= 2 * PI;return Rad;}。