高斯投影坐标正反算VB程序

⽤（VB）实现测量坐标转换系统坐标转换系统（VB）东华理⼯⼤学Theory 北京54⾼斯坐标转换西安80⾼斯坐标转换系统1.0版多点处理结果核⼼源代码：'⾼斯坐标转换成⼤地坐标过程Public Sub GausReverse(a As Double, f As Double, x() As Double, y() As Double, RB() As Double, RL() As Double, k As Integer)Dim i As Integer, fxb As Double, fxbl As Double, fybl As Doublee = Sqr(2 *f - f ^ 2)C = a / Sqr(1 - e ^ 2)e2 = e / Sqr(1 - e ^ 2)beita0 = 1 - (3 / 4) * e2 ^ 2 + (45 / 64) * e2 ^ 4 - (175 / 256) * e2 ^ 6 + (11025 / 16384) * e2 ^ 8 beita2 = beita0 - 1beita4 = (15 / 32) * e2 ^ 4 - (175 / 384) * e2 ^ 6 + (3675 / 8192) * e2 ^ 8beita6 = (-35 / 96) * e2 ^ 6 + (735 / 2048) * e2 ^ 8beita8 = (315 / 1024) * e2 ^ 8For i = 1 To kB0 = x(i) / (C * beita0)Dofxb = 0fxbl = 0fybl = 0t = Tan(B0)yita = e2 * Cos(B0)n = a / Sqr(1 - e ^ 2 * (Sin(B0)) ^ 2)a2 = (1 / 2) * n * Sin(B0) * Cos(B0)a3 = (1 / 6) * n * (Cos(B0)) ^ 3 * (1 - t ^ 2 + yita ^ 2)a4 = (1 / 24) * n * Sin(B0) * (Cos(B0)) ^ 3 * (5 - t ^ 2 + 9 * yita ^ 2 + 4 * yita ^ 4)a5 = (1 / 120) * n * (Cos(B0)) ^ 5 * (5 - 18 * t ^ 2 + t ^ 4 + 14 * yita ^ 2 - 58 * yita ^ 2 * t ^ 2) a6 = (1 / 720) * n * Sin(B0) * (Cos(B0)) ^ 5 * (61 - 58 * t ^ 2 + t ^ 4)fxb = fxb + (C * beita6 + C * beita8 * (Cos(B0)) ^ 2) * (Cos(B0)) ^ 2fxb = (fxb + C * beita4) * (Cos(B0)) ^ 2fxb = (fxb + C * beita2) * Sin(B0) * Cos(B0)fxbl = a2 * l0 ^ 2 + a4 * l0 ^ 4 + a6 * l0 ^ 6fybl = a3 * l0 ^ 3 + a5 * l0 ^ 5RB(i) = (x(i) - fxb - fxbl) / (C * beita0)a1 = (a * Cos(RB(i))) / Sqr(1 - e ^ 2 * (Sin(RB(i)) ^ 2))RL(i) = (y(i) - fybl) / a1If Abs(RB(i) - B0) <= 0.0000000001 And Abs(RL(i) - l0) <= 0.0000000001 ThenRL(i) = zrl + l0Exit DoElseB0 = RB(i)l0 = RL(i)End IfLoopNext iEnd Sub'⼤地坐标B，L转换为⾼斯坐标x,y的过程Public Sub BLHGaus(RB() As Double, RL() As Double, GX() As Double, GY() As Double, a As Double, f As Double, k As Integer)Dim l0 As Double, fxb As Double, gxbl As Double, fybl As Doublebeita0 = 1 - (3 / 4) * e2 ^ 2 + (45 / 64) * e2 ^ 4 - (175 / 256) * e2 ^ 6 + (11025 / 16384) * e2 ^ 8 beita2 = beita0 - 1beita4 = (15 / 32) * e2 ^ 4 - (175 / 384) * e2 ^ 6 + (3675 / 8192) * e2 ^ 8beita6 = (-35 / 96) * e2 ^ 6 + (735 / 2048) * e2 ^ 8beita8 = (315 / 1024) * e2 ^ 8For i = 1 To kl0 = RL(i) - zrln = a / Sqr(1 - e ^ 2 * (Sin(RB(i))) ^ 2)t = Tan(RB(i))yita = e2 * Cos(RB(i))a1 = (a * Cos(RB(i))) / Sqr(1 - e ^ 2 * (Sin(RB(i)) ^ 2))a2 = (1 / 2) * n * Sin(RB(i)) * Cos(RB(i))a3 = (1 / 6) * n * (Cos(RB(i))) ^ 3 * (1 - t ^ 2 + yita ^ 2)a4 = (1 / 24) * n * Sin(RB(i)) * (Cos(RB(i))) ^ 3 * (5 - t ^ 2 + 9 * yita ^ 2 + 4 * yita ^ 4) a5 = (1 / 120) * n * (Cos(RB(i))) ^ 5 * (5 - 18 * t ^ 2 + t ^ 4 + 14 * yita ^ 2 - 58 * yita ^ 2 * t ^ 2)a6 = (1 / 720) * n * Sin(RB(i)) * (Cos(RB(i))) ^ 5 * (61 - 58 * t ^ 2 + t ^ 4)fxb = 0fxbl = 0fybl = 0fxb = fxb + (C * beita6 + C * beita8 * (Cos(RB(i))) ^ 2) * (Cos(RB(i))) ^ 2fxb = (fxb + C * beita4) * (Cos(RB(i))) ^ 2fxb = (fxb + C * beita2) * Sin(RB(i)) * Cos(RB(i))fxbl = a2 * l0 ^ 2 + a4 * l0 ^ 4 + a6 * l0 ^ 6fybl = a3 * l0 ^ 3 + a5 * l0 ^ 5GX(i) = C * beita0 * RB(i) + fxb + fxblGY(i) = a1 * l0 + fyblNext iEnd Sub'三维直⾓坐标XYZ转换成⼤地坐标BLH过程Public Sub BLHXYZ1(SX() As Double, SY() As Double, SZ() As Double, RB() As Double, RL() As Double, RH() As Double, k As Integer, a As Double, f As Double)Dim i As IntegerDim N0 As Double, H0 As Double, B0 As Double, sb As Double, Ni As Doublesb = a * Sqr(1 - e ^ 2)pi = 4 * Atn(1)For i = 1 To kRL(i) = Atn(Abs(SY(i) / SX(i)))If SY(i) > 0 And SX(i) < 0 ThenRL(i) = pi - RL(i)End IfN0 = aH0 = Sqr(SX(i) ^ 2 + SY(i) ^ 2 + SZ(i) ^ 2) - Sqr(a * sb)B0 = Atn(SZ(i) / ((Sqr(SX(i) ^ 2 + SY(i) ^ 2)) * (1 - e ^ 2 * N0 / (N0 + H0))))DoNi = a / Sqr(1 - e ^ 2 * (Sin(B0)) ^ 2)RH(i) = Sqr(SX(i) ^ 2 + SY(i) ^ 2) / Cos(B0) - NiRB(i) = Atn(SZ(i) / ((Sqr(SX(i) ^ 2 + SY(i) ^ 2)) * (1 - e ^ 2 * Ni / (Ni + RH(i)))))If Abs(RB(i) - B0) < 0.0000000001 And Abs(RH(i) - H0) < 0.0000000001 ThenExit DoElseB0 = RB(i)H0 = RH(i)End IfLoopNext iEnd Sub'最⼩⼆乘法求解七参数布尔萨模型Public Sub SloveBuersa(XD() As Double, XG() As Double, R() As Double, k As Integer) Dim a0 As Double, a1 As Double, a2 As DoubleFor i = 1 To ka0 = XG(i, 1)a1 = XG(i, 2)a2 = XG(i, 3)Call Cjuzhen(a0, a1, a2, G())For j = 1 To 3Next sNext jNext iFor i = 1 To kFor j = 1 To 3LC(3 * (i - 1) + j) = XD(i, j) - XG(i, j) Next jNext iFor i = 1 To 3 * kFor j = 1 To 7ET(j, i) = EC(i, j)Next jNext iFor i = 1 To 7For j = 1 To 7CTC(i, j) = 0For s = 1 To 3 * kCTC(i, j) = CTC(i, j) + ET(i, s) * EC(s, j) Next sNext jNext iFor i = 1 To 7For j = 1 To 7CTC1(i, j) = CTC(i, j)Next jNext iCall Comm.Reverse(CTC(), 7)For i = 1 To 7CTL(i) = 0For j = 1 To 3 * kCTL(i) = CTL(i) + ET(i, j) * LC(j)Next jNext iFor j = 1 To 7R(i) = R(i) + CTC(i, j) * CTL(j)Next jNext iEnd Sub'布尔萨模型系数矩阵⼀部分Public Sub Cjuzhen(x As Double, y As Double, z As Double, D() As Double) For i = 1 To 3D(i, i) = 1Next iD(1, 4) = xD(1, 5) = 0D(1, 6) = -zD(1, 7) = yD(2, 4) = yD(2, 5) = zD(2, 6) = 0D(2, 7) = -xD(3, 4) = zD(3, 5) = -yD(3, 6) = xD(3, 7) = 0End Sub'使⽤七参数求解新坐标系下的坐标Public Sub UseBuersa(XD() As Double, XG() As Double, EC() As Double, R() As Double, k As Integer) Dim a1 As Double, a2 As Double, a3 As Double, CR(1000) As DoubleFor i = 1 To ka1 = XG(i, 1)a2 = XG(i, 2)a3 = XG(i, 3)Call Cjuzhen(a1, a2, a3, G())For j = 1 To 3For s = 1 To 7EC(3 * (i - 1) + j, s) = G(j, s)Next iFor i = 1 To 3 * kCR(i) = 0For j = 1 To 7CR(i) = CR(i) + EC(i, j) * R(j)Next jNext iFor i = 1 To kFor j = 1 To 3XD(i, j) = XG(i, j) + CR(3 * (i - 1) + j)Next jNext iEnd Sub'弧度化为度分秒'弧度化成⾓度Public Sub RuJiao(ByVal rudu As Double, jiaodu As Double) Dim ja As Integer, jb As Integer, jc As Double pi = 4 * Atn(1)jiaodu = rudu * 180 / pija = Fix(jiaodu)jb = Fix((jiaodu - ja) * 60)jc = ((jiaodu - ja) * 60 - jb) * 60jiaodu = ja + jb / 100 + jc / 10000End Sub'矩阵求逆Public Sub Reverse(Ba, n%)Dim k%, K1%, j%, i%Dim C As Double, Aa(100, 200)For i = 1 To nFor j = 1 To nAa(i, j) = Ba(i, j)Next jNext iFor i = 1 To nAa(i, j + n) = 0End IfNext jNext iFor k = 1 To nFor j = k To nIf Aa(j, k) <> 0 Then GoTo 200Next jMsgBox "逆矩阵不存在": Exit Sub 200: For i = 1 To 2 * n C = Aa(k, i)Aa(k, i) = Aa(j, i)Aa(j, i) = CNext iC = 1 / Aa(k, k)For j = 1 To 2 * nAa(k, j) = C * Aa(k, j)Next jFor K1 = 1 To nIf K1 <> k ThenC = -Aa(K1, k)For j = 1 To 2 * nAa(K1, j) = Aa(K1, j) + C * Aa(k, j)Next jEnd IfNext K1Next kFor i = 1 To nFor j = n + 1 To 2 * nAa(i, j - n) = Aa(i, j)Next jNext iNext iEnd Sub'⾓度化成弧度Public Sub JiaoHu(ByVal jiaodu As Double, hudu As Double) Dim ja As Double, jb As Double, jc As Double pi = 4 * Atn(1)ja = Fix(jiaodu)jb = Fix((jiaodu - ja) * 100)jc = ((jiaodu - ja) * 100 - jb) * 100jiaodu = ja + jb / 60 + jc / 3600hudu = jiaodu * pi / 180End Sub。

高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。

在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系)，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。

但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。

此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。

//6度带宽54北京坐标系//高斯投影由大地坐标(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)*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));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;}//高斯投影由经纬度(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)-(3 5*e2*e2*e2/3072)*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;}NN卯酉圈曲率半径，测量学里面用N表示M为子午线弧长，测量学里用大X表示fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示R为底点所对的曲率半径，测量学里用Nf表示。

坐标正算、反算计算方法及在Excel 中的VBA 编程测量中经常需要将某点相对坐标系坐标转换成线路的里程、偏距，或根据线路某一里程偏距计算出对应的相对坐标系坐标，为寻求一种快速简单高效的计算方法，本文对线路正算反算的原理进行了阐述，并结合Excel VBA 编程，将编程和Excel 的拖拽的功能相结合，编制出实用计算表，特别适用于需要大量计算边桩、围护桩的情况。

关键词：坐标方位角坐标正算坐标反算 V AB 编程循环迭代直接算法一、坐标方位角的反算1.坐标方位角反算如图1所示，已知点A 、B 的坐标，求直线AB坐标方位角α。

图1坐标方位角反算直线AB 之间的坐标增量：AB B AAB B Ax x x y y y ∆=−∆=−当0,0AB AB x y ∆>∆>时，角α位于第一象限角：arctan ABABy x α∆=∆当0,0AB AB x y ∆<∆>时，角α位于第二象限角：arctan 180AB ABy x α∆=+°∆当0,0AB AB x y ∆<∆<时，角α位于第三象限角：arctan 180AB ABy x α∆=+°∆当0,0AB AB x y ∆>∆<时，角α位于第二象限角：arctan360AB AB y x α∆=+°∆2.坐标方位角反算的VBA 编程可用VBA 将上述过程定义为一个名为angel()的函数，代码如下：Function angel(x0As Double, y0 As Double, x1 As Double, y1 As Double) As Double dx = x1- x0dy = y1- y0If dx > 0 And dy > 0 Thenangel = Atn(dy / dx)End IfIf dx < 0 And dy > 0 Thenangel = Atn(dy / dx) + 3.14159265358979End IfIf dx < 0 And dy < 0 Thenangel = Atn(dy / dx) + 3.14159265358979End IfIf dx > 0 And dy < 0 Thenangel = Atn(dy / dx) + 3.14159265358979 * 2End IfEnd Function二、直线段坐标正算与反算1.直线段正算图2直线段计算已知HZ 点坐标（x1,y1）、里程N HZ ，ZH 点坐标（x2,y2），正算时已知P 点对应的中桩里程Np 和偏距e （规定沿着线路前进方向，左边偏距为负，右边偏距为正），Np>N HZ ，求P 点对应的坐标。

#include<stdio.h>#include<math.h>#include<stdlib.h>double trans1() //由度分秒到弧度{double B1,B11,B12,B13,B111;scanf("%lf°%lf′%lf″",&B11,&B12,&B13);B111=fabs(B11); //B11可能为负值B1=B111+B12/60.0+B13/3600.0;B1=B1*atan(1)/45.0;if(B11<0)B1=-B1;return B1;}void trans2(double B) //由弧度到度分秒并输出角度值{int a,b;double B0;B0=fabs(B); //B可能为负值double c;B0=B0*45.0/atan(1);a=int(B0);b=int((B0-a)*60);c=(B0-a)*3600-b*60;if((int)(c)==60) //为了避免出现59′60″这种情况，不过好像不起作用，不知道为什么,原来是int没有加括号{b=b+1;c=0.0;}if(b==60){b=0;a=a+1;}if(B<0)a=-a;printf("%d°%d′%.4f″\n",a,b,c);}void Gauss1() //高斯投影坐标正算，克拉索夫斯基椭球{double B,L; //输入经纬度并转换成弧度printf("请输入大地纬度B=\n");B=trans1();printf("请输入大地经度L=\n");L=trans1();double L0,l;printf("请以度分秒的形式（如10°10′10″）输入中央子午线L0=\n");L0=trans1();l=L-L0; //经度与中央子午线之差……坑爹的书上居然又错了!!114°20′0″-117°0′0″=3°20′0″FUCK! 原来是我误会了……人家的L0=111°double R; //平面子午线收敛角γ（仅与B和l有关）R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l*l)*l*l*cos(B)*cos(B))* l*sin(B);double N,a0,a4,a6,a3,a5; //求解一些常数以及系数N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);double x,y; //求解Gauss投影平面坐标想，x,yx=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);printf("平面坐标系中的坐标\nx=%.4fM\n",x);printf("平面坐标系中的坐标\ny=%.4fM\n",y);printf("平面子午线收敛角R=\n");trans2(R);}void Gauss2() //Gauss投影坐标反算，克拉索夫斯基椭球{double x,y,L0; //输入平面坐标x,y及中央子午线L0=printf("请输入平面坐标系中的坐标x=\n");scanf("%lf",&x);printf("请输入平面坐标系中的坐标y=\n");scanf("%lf",&y);printf("请以度分秒的形式输入（如10°10′10″）输入中央子午线L0=\n");L0=trans1();double beta,Bf,Nf,Z;beta=x/6367558.4969;Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(beta)*co s(beta))*1e-10*cos(beta)*sin(beta);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*cos(B f);Z=y/(Nf*cos(Bf));double b2,b3,b4,b5;b2=(0.5+0.003369*cos(Bf)*cos(Bf))*cos(Bf)*sin(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);double B,L,l,R;B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;l=(1-(b3-b5*Z*Z)*Z*Z)*Z;L=L0+l;R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l*l)*l*l*cos(B)*cos(B))* l*sin(B);printf("大地纬度B=\n");trans2(B);printf("大地经度L=\n");trans2(L);printf("平面子午收敛角R=\n");trans2(R);}void Gauss3() //换带计算，克拉索夫斯基椭球{double x1,y1,L10;printf("请输入平面坐标系中的原坐标x1=\n");scanf("%lf",&x1);printf("请输入平面坐标系中的原坐标y1=\n");scanf("%lf",&y1);printf("请以度分秒的形式输入（如10°10′10″）原中央子午线L10=\n");L10=trans1();double beta,Bf,Nf,Z;beta=x1/6367558.4969;Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(beta)*co s(beta))*1e-10*cos(beta)*sin(beta);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*cos(B f);Z=y1/(Nf*cos(Bf));double b2,b3,b4,b5;b2=(0.5+0.003369*cos(Bf)*cos(Bf))*cos(Bf)*sin(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);double B,L,l1;B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;l1=(1-(b3-b5*Z*Z)*Z*Z)*Z;L=L10+l1;double L20,l2,R;printf("请以度分秒的形式（如10°10′10″）输入新中央子午线L20=\n");L20=trans1();l2=L-L20; //R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l2*l2)*l2*l2*cos(B)*cos (B))*l2*sin(B);//子午线收敛角γ（仅与B和l有关）double N,a0,a4,a6,a3,a5; //求解一些常数以及系数N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);double x2,y2;x2=6367558.4969*B-(a0-(0.5+(a4+a6*l2*l2)*l2*l2)*l2*l2*N)*sin(B)*cos(B);y2=(1+(a3+a5*l2*l2)*l2*l2)*l2*N*cos(B);printf("平面坐标系中的坐标x2=%.4fM\n",x2);printf("平面坐标系中的坐标y2=%.4fM\n",y2);printf("平面子午线收敛角R=\n");trans2(R);printf("大地纬度B=\n");trans2(B);printf("大地经度L=\n");trans2(L);}void main(){int k;printf("请输入k=\n");printf("注：k=1表示正算（BL→xy）k=2表示反算（xy→BL）k=3表示换带（x1y1→x2y2）\n");scanf("%d",&k);if(k==1)Gauss1();if(k==2)Gauss2();if(k==3)Gauss3();}。

高斯正反算实验报告一、实验目的：通过对高斯正反算的软件编程，增强学生对《大地测量学基础》课程的理解，使学生牢固掌握高斯正反算的基本原理和公式，熟悉高斯正反算的基本技能和计算方法。

培养学生利用计算机编程的能力、综合运用知识的能力及理论联系实际的能力二、实验要求：要求学生综合运用测绘知识、大地测量知识、数学知识和计算机知识，设计数学模型和程序算法，编制程序实现数据的自动化处理。

三、实验环境Visual studio 2008界面如下：输入数据计算:代码如下：Public Class高斯投影正反算Private Sub高斯投影正反算_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load'自动化控件显示初始值RadioButton2.Checked = TrueRadioButton3.Checked = TrueTextBox1.Focus()TextBox1.Text = "0"TextBox2.Text = "0"TextBox3.Text = "0"TextBox4.Text = "0"TextBox5.Text = "0"TextBox6.Text = "0"RichTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标"End SubPrivate Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click'获取输入数据Dim bb1 As Double, bb2 As Double, bb3 As Double, ll1 As Double, ll2 As Double, ll3 As Doublebb1 = TextBox1.Textbb2 = TextBox2.Textbb3 = TextBox3.Textll1 = TextBox4.Textll2 = TextBox5.Textll3 = TextBox6.TextDim B As Double, L As Double, x As Double, y As Double, X1 As Double, Y1 As Double'X1 Y1Dim N As IntegerDim L0 As Double, l1 As DoubleDim p As Doublep = 206264.80625'检查输入格式的正确性If (bb1 >= 0 And bb1 < 90 And bb2 >= 0 And bb2 < 60 And bb3 >= 0 And bb3 < 60) Then B = bb1 * 3600 + bb2 * 60 + bb3ElseRichTextBox1.Text = "纬度输入格式不正确！"ReturnEnd IfIf (ll1 >= 0 And ll1 < 360 And ll2 >= 0 And ll2 < 60 And ll3 >= 0 And ll3 < 60) Then L = ll1 * 3600 + ll2 * 60 + ll3ElseRichTextBox1.Text = "经度输入格式不正确！"ReturnEnd IfDim b1 As Double'b1b1 = B / p'获取用户选项值Dim tyd As IntegerIf (RadioButton2.Checked = True) Thentyd = 6Elsetyd = 3End IfDim ty As IntegerIf (RadioButton3.Checked = True) Thenty = 1Elsety = 2End If'按带求带号以及中央子午线经度()If (tyd = 6) ThenN = Int(L / (6 * 3600.0)) + 1L0 = (6 * N - 3) * 3600.0l1 = (L - L0) / pElseIf (tyd = 3) ThenN = Int(L / (3 * 3600.0) + 0.5)L0 = 3 * N * 3600.0l1 = (L - L0) / pEnd If'采用克氏椭球参数的计算公式If (ty = 1) ThenDim c As DoubleDim c1 As DoubleDim c2 As DoubleDim l2 As DoubleDim n1 As DoubleDim a0 As DoubleDim a3 As DoubleDim a4 As DoubleDim a5 As DoubleDim a6 As DoubleDim i As Integeri = 1c = Math.Pow(Math.Cos(b1), 2)c1 = Math.Sin(b1) * Math.Cos(b1)c2 = Math.Cos(b1)l2 = Math.Pow(l1, 2)n1 = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c 'n1a0 = 32140.404 - (135.3302 - (0.7092 - 0.004 * c) * c) * ca4 = (0.25 + 0.00252 * c) * c - 0.04166a6 = (0.166 * c - 0.084) * ca3 = (0.3333333 + 0.001123 * c) * c - 0.1666667a5 = 0.0083 - (0.1667 - (0.1968 + 0.004 * c) * c) * cx = 6367558.4969 * b1 - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n1) * c1y = (1 + (a3 + a5 * l2) * l2) * l1 * n1 * c2X1 = xDim y11 As Doubley11 = y + 500000.0If y11 / i > 1 ThenY1 = N * i * 10 + y11 'Y坐标加代号i = i * 10End IfDim tuoqiu As Stringtuoqiu = "采用克氏椭球参数,"RichTextBox1.Text = "按经差" + tyd.ToString("0") + "°进行投影分带，" + tuoqiu + "其计算结果为：" + vbCrLf + "x=" + x.ToString("0.00000000") + vbCrLf + "y=" +y.ToString("0.00000000") + vbCrLf + "国家统一坐标为："+ vbCrLf + "X="+ X1.ToString("0.00000000") + vbCrLf + "Y=" + N.ToString + Y1.ToString("0.00000000")End If'采用1975国际椭球参数的计算公式()If (ty = 2) ThenDim c As DoubleDim c1 As DoubleDim c2 As DoubleDim l2 As DoubleDim n1 As DoubleDim a0 As DoubleDim a3 As DoubleDim a4 As DoubleDim a5 As DoubleDim a6 As DoubleDim i As Integeri = 1c = Math.Cos(b1) * Math.Cos(b1)c1 = Math.Sin(b1) * Math.Cos(b1)c2 = Math.Cos(b1)l2 = Math.Pow(l1, 2)n1 = 6399596.652 - (21565.045 - (108.996 - 0.603 * c) * c) * ca0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * c) * c) * ca4 = (0.25 + 0.00253 * c) * c - 0.04167a6 = (0.167 * c - 0.083) * ca3 = (0.3333333 + 0.001123 * c) * c - 0.1666667a5 = 0.00878 - (0.1702 - 0.20382 * c) * cx = 6367452.1328 * b1 - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n1) * c1y = (1 + (a3 + a5 * l2) * l2) * l1 * n1 * c2X1 = xDim y11 As Doubley11 = y + 500000.0If y11 / i > 1 ThenY1 = N * i * 10 + y11 'Y坐标加代号i = i * 10End IfDim tuoqiu As Stringtuoqiu = "采用1975国际椭球参数,"RichTextBox1.Text = "按经差" + tyd.ToString("0") + "°进行投影分带，" + tuoqiu + "其计算结果为：" + vbCrLf + "x=" + x.ToString("0.00000000") + vbCrLf + "y=" +y.ToString("0.00000000") + vbCrLf + "国家统一坐标为："+ vbCrLf + "X="+ X1.ToString("0.00000000") + vbCrLf + "Y=" + N.ToString + Y1.ToString("0.00000000")End IfEnd SubPrivate Sub Button2_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button2.ClickDim y As Double, X1 As Double, Y1 As Double, B As Double, L As DoubleX1 = TextBox7.Text '获取输入数据Y1 = TextBox8.TextDim L0 As DoubleDim i As Integeri = 1Do While Y1 / i >= 10y = Y1 - Int(Y1 / i) * i - 500000L0 = 6 * Int(Y1 / i) - 3i = i * 10Loop'对Y坐标处理并求出中央子午线经度'按6°带克氏椭球反算Dim bt As Double, BT1 As Double, c3 As Double, c4 As Double, Bf As Double, c5 As Double, c6 As Double, Nf As Double, Z As DoubleDim b2 As Double, b3 As Double, b4 As Double, b5 As Double, z2 As DoubleDim p As Doublep = 206264.80625bt = X1 / 6367558.4969 * pBT1 = X1 / 6367558.4969c3 = Math.Cos(BT1) * Math.Cos(BT1)c4 = Math.Sin(BT1) * Math.Cos(BT1)Bf = (bt + (50221746 + (293622 + (2350 + 22 * c3) * c3) * c3) * c4 * Math.Pow(10, -10) * p) / pc5 = Math.Pow(Math.Cos(Bf), 2)c6 = Math.Sin(Bf) * Math.Cos(Bf)Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * c5) * c5) * c5Z = y / (Nf * Math.Cos(Bf))b2 = (0.5 + 0.003369 * c5) * c6b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5b5 = 0.2 - (0.1667 - 0.0088 * c5) * c5z2 = Math.Pow(Z, 2)B = (Bf * p - (1 - (b4 - 0.12 * z2) * z2) * z2 * b2 * p) / 3600.0L = L0 + ((1 - (b3 - b5 * z2) * z2) * Z * p) / 3600.0RichTextBox2.Text = "按6°带克氏椭球反算后，结果为：" + vbCrLf + "B=" +B.ToString("0.0000000000000") + vbCrLf + "L=" + L.ToString("0.0000000000000")End SubPrivate Sub Button3_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button3.ClickTextBox1.Focus()TextBox1.Text = "0"TextBox2.Text = "0"TextBox3.Text = "0"TextBox4.Text = "0"TextBox5.Text = "0"TextBox6.Text = "0"RichTextBox1.Text = " "End SubPrivate Sub Button4_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button4.ClickTextBox7.Text = " "TextBox8.Text = " "RichTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标"End SubPrivate Sub Button5_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button5.ClickEndEnd SubEnd Class。

vb坐标正算程序
VB坐标正算程序是一种用于计算点位坐标的工具，通常用于工程测量、地理信息系统等领域。

通过输入已知点的坐标和距离、方位角等数据，程序可以准确计算出目标点的坐标。

下面将详细介绍VB 坐标正算程序的使用方法和注意事项。

打开VB坐标正算程序，通常程序界面会包括输入框和计算按钮。

在输入框中，需要输入已知点的坐标、距离和方位角等数据。

确保输入的数据准确无误，否则会影响计算结果的准确性。

接下来，点击计算按钮，程序会根据输入的数据进行计算，最终显示目标点的坐标。

在计算过程中，程序会考虑各种因素，如坐标系、单位制等，确保计算结果符合实际需求。

在使用VB坐标正算程序时，需要注意以下几点：
1. 确保输入的数据准确无误，包括坐标、距离、方位角等信息。

2. 注意选择合适的坐标系和单位制，以保证计算结果的准确性。

3. 在计算过程中，及时保存已知点和目标点的数据，方便后续查阅和使用。

4. 如遇到计算结果不符合预期的情况，及时检查输入数据和计算方法，找出错误并进行修正。

总的来说，VB坐标正算程序是一种方便快捷的工具，可以帮助工程师、测量员等专业人士准确计算点位坐标。

通过正确使用该程序，
可以提高工作效率，减少人为误差，确保测量数据的准确性和可靠性。

希望以上介绍能帮助大家更好地理解和应用VB坐标正算程序。

vb坐标正算程序
VB坐标正算程序是一种非常实用的工具，它可以帮助我们快速计算出某个点的坐标。

在实际工作中，我们经常需要用到这种工具，比如在地图制作、测量、建筑设计等领域。

VB坐标正算程序的实现原理是利用数学公式来计算出目标点的坐标。

具体来说，我们需要知道已知点的坐标、距离和方位角，然后根据三角函数公式来计算出目标点的坐标。

这个过程需要用到VB 语言的数学函数和逻辑运算符，因此需要一定的编程基础。

在编写VB坐标正算程序时，我们需要注意以下几点：
1. 确定计算公式：根据已知点的坐标、距离和方位角，确定计算目标点坐标的公式。

这个公式需要考虑到坐标系的不同，比如笛卡尔坐标系和极坐标系的计算公式是不同的。

2. 输入数据的格式：在编写程序时，需要考虑输入数据的格式，比如坐标的单位、距离的单位、方位角的单位等。

这些单位需要在程序中进行转换，以保证计算的准确性。

3. 界面设计：为了方便用户使用，我们需要设计一个简洁明了的界面，让用户能够方便地输入数据和查看计算结果。

界面设计需要考虑到用户的使用习惯和操作流程，以提高用户的体验。

4. 错误处理：在编写程序时，需要考虑到可能出现的错误情况，比
如输入数据错误、计算公式错误等。

我们需要在程序中加入相应的错误处理机制，以避免程序崩溃或计算结果错误。

VB坐标正算程序是一种非常实用的工具，它可以帮助我们快速计算出目标点的坐标。

在编写程序时，我们需要考虑到计算公式、输入数据的格式、界面设计和错误处理等方面，以保证程序的准确性和易用性。

一、输入界面1.单击VB中的“运行”快捷键，弹出图1所示的运行界面图12.在图1中选择按钮，进入坐标正算模式，如图2图2在图2中输入已知点坐标、已知点至未知点的边长和坐标方位角、保留的小数位数，点名可不输入，输入完后，选择“计算”按钮。

如要计算多个点坐标，则计算完一个点后用鼠标单击“刷新”按钮，重复以上操作即可。

若想退出坐标正算功能，用鼠标左键单击“退出”按钮，在弹出的提示框中选择“是”，如图3图33. 在图1中选择按钮，进入坐标反算模式，如图4图4在图4中输入两个已知点坐标、保留的小数位数，点名可不输入，输入完后，选择“计算”按钮。

如要计算多个点坐标，则计算完一个点后用鼠标单击“刷新”按钮，重复以上操作即可。

若想退出坐标反算功能，用鼠标左键单击“退出”按钮，在弹出的提示框中选择“是”，如图3二、输入，输出数据表格中的方位角以小数的形式输入，例如“321°18′56″“的输入格式为321.1856，若度数为321°18′56.5″的输入格式为321.18565，坐标的保留位数保留三位，别的以此类推方位角以小数的形式输入，例如“164°02′02″”的输入格式为164.0202，若度数为164°02′02.5″的输入格式为164.02025，表格中的方位角的小数保留位数为0位三、源代码Private Sub Command1_Click()If Option1.Value = True And Option2.Value = False Thenzhengsuan '当选择坐标正算按钮时调用坐标正算程序End IfIf Option1.Value = False And Option2.Value = True Thenfansuan '当选择坐标反算按钮时调用坐标反算程序End IfEnd Sub'坐标正算程序Private Sub zhengsuan()Dim s As DoubleDim a As DoubleDim x As DoubleDim y As DoubleDim m As DoubleDim n As DoubleDim degree As DoubleDim minute As DoubleDim second As DoubleDim rad As DoubleDim bbt As IntegerDim result As DoubleDim p As Doubles = Val(Text5.Text)a = Val(Text6.Text)Rad_do a, degree, minute, second, rad '调用将度分秒转化为弧度的程序x = Val(Text2.Text) + s * Cos(rad)y = Val(Text3.Text) + s * Sin(rad)p = Fix(x)x = x - pbbt = Val(Text4.Text)Sheru m, n, bbt, x, result '调用奇进偶舍程序Label12.Caption = result + pp = Fix(y)y = y - pSheru m, n, bbt, y, result '调用奇进偶舍程序Label13.Caption = result + pEnd Sub'坐标反算程序Private Sub fansuan()Dim s As DoubleDim x As DoubleDim y As DoubleDim sing As DoubleDim m As DoubleDim n As DoubleDim f As DoubleDim result As DoubleDim bbt As IntegerDim degree As DoubleDim minute As DoubleDim second As DoubleDim rad As DoubleDim p As Doublex = Val(Text5.Text)y = Val(Text6.Text)m = y - Val(Text3.Text)n = x - Val(Text2.Text)If n <> 0 Thenm = m / nrad = Atn(m)sing = Sgn(rad)End Ifdegree_m_s degree, minute, second, rad '调用将弧度转化为度分秒的程序bbt = Val(Text8.Text)m = second * 10 ^ bbtp = Fix(second)second = second - pSheru m, n, bbt, second, result '调用奇进偶舍程序second = result + pfangweijiao degree, minute, second, (x - Val(Text2.Text)), (y - Val(Text3.Text)), sing '调用计算坐标方位角的程序zhuanhua minute, second, degreesecond = second * 10 ^ bbtLabel13.Caption = degree & "." & minute & secondIf minute < 10 ThenLabel13.Caption = degree & "." & "0" & minute & secondEnd IfIf second < 10 ^ (bbt + 1) ThenLabel13.Caption = degree & "." & minute & "0" & secondEnd IfIf minute < 10 And second < 10 ^ (bbt + 1) ThenLabel13.Caption = degree & "." & "0" & minute & "0" & secondEnd If'计算边长Ss = (y - Val(Text3.Text)) ^ 2 + (x - Val(Text2.Text)) ^ 2s = Sqr(s)bbt = Val(Text4.Text)p = Fix(s)s = s - pSheru m, n, bbt, s, result '调用奇进偶舍程序Label12.Caption = result + pEnd Sub'将度分秒化为弧度Private Sub Rad_do(ByVal a As Double, ByVal degree As Double, ByVal minute As Double, ByVal second As Double, ByRef rad As Double)degree = a \ 1a = a - Fix(a)minute = Fix(a * 100)second = (a * 100 - minute) * 100rad = (3600 * degree + 60 * minute + second) / 206264.8063End Sub'将弧度转化为度分秒的程序Private Sub degree_m_s(ByRef degree As Double, ByRef minute As Double, ByRef second As Double, ByRef rad As Double)rad = rad * 180 / 3.1415926535degree = rad \ 1rad = (rad - degree) * 60minute = rad \ 1second = (rad - minute) * 60End Sub'奇进偶舍程序Private Sub Sheru(ByVal m As Double, ByVal n As Double, ByVal bbt As Double, ByVal x As Double, ByRef result As Double)m = x * 10 ^ bbtn = m - m \ 1If n < 0.5 Thenresult = (m \ 1) / 10 ^ bbtElseIf n > 0.5 Thenresult = (m \ 1 + 1) / 10 ^ bbtElseIf (m \ 1) Mod 2 Thenresult = (m \ 1 + 1) / 10 ^ bbtElseresult = (m \ 1) / 10 ^ bbtEnd IfEnd IfEnd Sub'计算坐标方位角的程序Private Sub fangweijiao(ByRef degree As Double, ByRef minute As Double, ByRef second As Double, ByVal m, ByVal n As Double, ByVal sing As Integer)Dim i As IntegerIf sing = 1 ThenIf n < 0 And m < 0 Thendegree = degree + 180End IfElseIf sing = -1 ThenIf n > 0 And m < 0 Thendegree = degree + 179minute = minute + 59second = second + 60ElseIf n < 0 And m > 0 Thendegree = degree + 359minute = minute + 59second = second + 60End IfEnd Ifi = Sgn(n)If m <> 0 ThenIf m < 0 And n = 0 Thendegree = "180"ElseIf m > 0 And n = 0 Thendegree = "0"End IfElseIf i = 1 Thendegree = "90"ElseIf i = -1 Thendegree = "270"ElseMsgBox "您输入了两个相同的点，请重新输入！"End IfEnd IfEnd Sub'当分秒超过60时须向上一级进位及方位角度数超过360°须减360°的程序Private Sub zhuanhua(ByRef minute As Double, ByRef second As Double, ByRef degree As Double)If second >= 60 Thenminute = minute + 1second = second - 60ElseIf second < 0 Thenminute = minute - 1second = second + 60End IfIf minute >= 60 Thendegree = degree + 1minute = minute - 60ElseIf minute < 0 Thendegree = degree - 1minute = minute + 60End IfIf degree >= 360 Thendegree = degree - 360End IfEnd Sub'退出程序Private Sub Command3_Click()If MsgBox("是否退出？", vbYesNo, "提示") = vbYes ThenUnload MeEnd IfEnd Sub'刷新程序Private Sub Command4_Click()Text2.Text = ""Text3.Text = ""Text4.Text = 3Text5.Text = ""Text6.Text = ""Text8.Text = 1Label12.Caption = ""Label13.Caption = ""End Sub'设置窗体的大小，使窗体充满整个屏幕并对label6和label7赋值Private Sub Form_Load()Me.Height = Screen.HeightMe.Width = Screen.WidthMe.Left = 0Me.Top = 0Label6.Caption = "边长（S）"Label7.Caption = "方位角（a）"Label17.Caption = "陈亮编程"Label18.Caption = " 2011.09.05"End Sub'设置坐标正算时Label6.Caption和Label7.Caption的值及设计小数点位数Private Sub Option1_Click()If Option1.Value = True And Option2.Value = False ThenLabel14.Caption = Text1.Text & "—>未知点" & Text7.TextLabel6.Caption = "边长（S）"Label7.Caption = "方位角（a）"Label10.Caption = "X"Label11.Caption = "Y"Text4.Text = 3Text8.Text = 1End IfEnd Sub''设置坐标反算时Label0.Caption和Label1.Caption的值及设计小数点位数Private Sub Option2_Click()If Option1.Value = False And Option2.Value = True ThenLabel6.Caption = "X"Label7.Caption = "Y"Label14.Caption = ""Label10.Caption = "边长（S）"Label11.Caption = "方位角（a）"Text4.Text = 3Text8.Text = 1End IfEnd Sub。

高斯消去法Private Sub Command1_Click()Dim a(1 To 10, 1 To 11) As DoubleDim x(1 To 10) As DoubleDim Sum As DoubleDim n As IntegerDim i As IntegerDim j As IntegerDim k As Integern = Val(InputBox("输入未知量个数(最多10个)")) For i = 1 To nFor j = 1 To n + 1a(i, j) = Val(InputBox("输入增广矩阵")) Next jNext iFor i = 1 To nFor j = 1 To n + 1Print a(i, j);Next jPrintNext iRem 消元过程For k = 1 To n - 1For i = k + 1 To nFor j = k + 1 To n + 1a(i, j) = a(i, j) - a(i, k) * a(k, j) / a(k, k) Next jNext iNext kRem 回代过程x(n) = a(n, n + 1) / a(n, n)Print "x"; n; "="; x(k);PrintFor k = n - 1 To 1 Step -1Sum = 0For j = k + 1 To nSum = a(k, j) * X(j) + SumNext jx(k) = (a(k, n + 1) - Sum) / a(k, k)Print "x"; k; "="; x(k);PrintNext kEnd Sub'将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度Public Function DoToHu(ByVal DoFenMiao As Double) As Single Dim du%, fen%, miao%, angle#du = Fix(DoFenMiao)DoFenMiao = (DoFenMiao - du) * 100fen = Fix(DoFenMiao)miao = (DoFenMiao - fen) * 100angle = du + fen / 60 + miao / 3600DoToHu = angle * PI / 180End Function'弧度化为度.分秒的形式：输入弧度值，输出度.分秒（各占两位）Public Function HuToDo(ByVal Hu As Double) As SingleDim du%, fen%, miao%Hu = Hu * 180 / PIdu = Fix(Hu)Hu = (Hu - du) * 60fen = Fix(Hu)Hu = (Hu - fen) * 60miao = Fix(Hu + 0.5)If miao = 60 Thenfen = fen + 1miao = 0End IfIf fen = 60 Thendu = du + 1fen = 0End IfHuToDo = du + fen / 100 + miao / 10000End Function矩阵运算Dim a() As LongDim b() As LongDim c() As LongDim m As LongDim n As LongDim p As LongDim q As LongPrivate Sub Command1_Click()On Error GoTo lab:m = InputBox("请输入A矩阵行数", "提示")n = InputBox("请输入A矩阵列数", "提示")Picture1.ClsPicture3.ClsReDim a(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To na(i, j) = InputBox("请输入矩阵a(" & i & "," & j & ")数值", "提示") Picture1.Print a(i, j);Next jPicture1.PrintNext ilab:End SubPrivate Sub Command10_Click()EndEnd SubPrivate Sub Command2_Click()On Error GoTo lab:p = InputBox("请输入B矩阵行数", "提示")q = InputBox("请输入B矩阵列数", "提示")Picture2.ClsPicture3.ClsReDim b(1 To p, 1 To q) As LongFor i = 1 To pFor j = 1 To qb(i, j) = InputBox("请输入矩阵b(" & i & "," & j & ")数值", "提示") Picture2.Print b(i, j);Next jPicture2.PrintNext ilab:End SubPrivate Sub Command3_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf m <> p Or n <> q ThenMsgBox "请输入行数和列数相同的矩阵才可相加", vbOKOnly, "提示" End IfIf m = p And n = q ThenLabel1.Caption = "+"ReDim c(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To nc(i, j) = a(i, j) + b(i, j)Picture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End SubPrivate Sub Command4_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf m <> p Or n <> q ThenMsgBox "请输入行数和列数相同的矩阵才可相减", vbOKOnly, "提示" End IfIf m = p And n = q ThenLabel1.Caption = "-"ReDim c(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To nc(i, j) = a(i, j) - b(i, j)Picture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End SubPrivate Sub Command5_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf n <> p ThenMsgBox "请输入A矩阵列数和B矩阵行数相等的矩阵再做乘积", vbOKOnly, "提示" End IfIf n = p ThenLabel1.Caption = "x"ReDim c(1 To m, 1 To q) As LongFor i = 1 To mFor j = 1 To qFor k = 1 To nc(i, j) = c(i, j) + a(i, k) * b(k, j)Next kPicture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End Sub逆矩阵ReDim a(1 To n, 1 To n * 2)ReDim x(1 To n)For i = 1 To nFor j = 1 To n * 2a(i, j) = Val(InputBox("输入增广矩阵第" & i & "行第" & j & "列元素"))Picture1.Print a(i, j);Next jPicture1.PrintNext iFor k = 1 To n - 1 '消元For i = k + 1 To nFor j = k + 1 To n * 2a(i, j) = a(i, j) - a(i, k) * a(k, j) / a(k, k) Next jNext iNext kFor j = n + 1 To 2 * na(n, j) = a(n, j) / a(n, n)Next jFor k = n To 1 Step -1 '反消元For i = k - 1 To 1 Step -1For j = k + 1 To n * 2a(i, j) = (a(i, j) - a(i, k) * a(k, j)) / a(i, i) Next jNext iNext kFor i = 1 To nFor j = n + 1 To 2 * nPicture2.Print a(i, j); Space(3);Next jPicture2.PrintNext ilab:End Sub坐标正反算Private Sub Command1_Click()XX = Val(Text3.Text) - Val(Text1.Text)YY = Val(Text4.Text) - Val(Text2.Text)If XX = 0 And YY = 0 ThenMsgBox ("A与B点不能为同一点，请重新输入") Text1.Text = ""Text2.Text = ""Text3.Text = ""Text4.Text = ""ElseIf YY = 0 And XX > 0 ThenText5.Text = "00"Text6.Text = "00"Text7.Text = "00"ElseIf XX = 0 And YY > 0 ThenText5.Text = 90Text6.Text = "00"Text7.Text = "00"ElseIf XX < 0 And YY = 0 ThenText5.Text = 180Text6.Text = "00"Text7.Text = "00"ElseIf XX = 0 And YY < 0 ThenText5.Text = 270Text6.Text = "00"Text7.Text = "00"Elsek = Abs(YY / XX)R = Atn(k) * 180 / 3.1415926If XX > 0 And YY > 0 Thena = RElseIf XX < 0 And YY > 0 Thena = 180 - RElseIf XX < 0 And YY < 0 Thena = 180 + RElseIf XX > 0 And YY < 0 Then ‘AB在第四象限a = 360 – REnd Ifb = Format(a, "0.00000") ‘取小数后5位Text5.Text = Fix(b) ‘求度数c = Fix((b - Text5.Text) * 60)Text6.Text = c ‘求分d = Round(((b - Text5.Text) * 60 - c) * 60)Text7.Text = d ‘求秒End IfEnd SubMatlabX=[107563.8100,107620.9521,109989.2770,111411.7664,109584.1244,109506.4700];Y=[571684.5200,568999.8520,568164.3993,569885.3537,572397.4802,570099.6000];a=[0,1,0,0,1,1;1,0,1,0,0,1;0,1,0,1,0,1;0,0,1,0,1,1;1,0,0,1,0,1;1,1,1,1,1,0][m,n]=size(a)plot(Y(1),X(1),'r^',Y(2),X(2),'bo',Y(3),X(3),'bo',Y(4),X(4),'bo',Y(5),X(6),'bo',Y(6),X(6),'r^') hold onfor i=1:mfor j =1:nif a(i,j)==1line([Y(i),Y(j)],[X(i),X(j)])hold onendendendtext(Y(1),X(1),'AA')text(Y(2),X(2),'11')text(Y(3),X(3),'22')text(Y(4),X(4),'33')text(Y(5),X(6),'44')text(Y(6),X(5),'BB')grid。

高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到.这不进行详细的推导.只是列出基本公式指导编程的进行。

二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句.本文中是利用C++语言进行基本的设计。

高斯正算中对椭球参数和带宽的选择主要运用了选择语句。

而高斯反算中除了选择语句的应用.在利用迭代算法求底点纬度还应用了循环语句。

编程中还应特别注意相关的度分秒和弧度之间的相互转换.这是极其重要的。

（2）相关流程图1）正算2）反算三．编程的相关代码（1）正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]) {FILE *r=fopen("a.txt","r"); assert(r!=NULL);FILE *w=fopen("b.txt","w"); assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L 0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系：北京54=1.西安80=2.WGS84=3：");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一.第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*m6/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B) )*100.0)*60.0+(g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L) )*100.0)*60.0+(g[i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i] .L0))*100.0)*60.0+(g[i].L0*100-(int)(g[i].L0*100))*100 .0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6 +a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽：3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n *n*n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/7 20;printf("\n所选坐标系的转换结果：x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}（2）反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量：x,y坐标和中央子午线经度//main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L 0)!=EOF)//文件为空.无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径.子午圈曲率半径// double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球：1=克拉索夫斯基椭球.2=1975国际椭球.3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球.求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B *C)-g[i].y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);//利用公式求解经纬度// int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果：%d度%d分%lf秒 %d 度%d分%lf秒 \n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*m6/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*s in(8*B1)/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2;}。