高斯投影正反算编程(可编辑修改word版)

大地测量高斯投影正反算程序代码课程：大地测量学基础姓名：林江伟学号：2008301610045班级： 0804 界面如下：输入数据计算：using System;using System.Collections.Generic; using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;namespace大地{public partial class Form1 : Form {double B, L;double x, y;double X, Y;int N;double L0, l1;double p=206264.80625;public Form1(){InitializeComponent();}private void Form1_Load(object sender, EventArgs e){//自动化控件显示初始值radioButton2.Checked=true;radioButton3.Checked = true;textBox1.Focus();textBox1.Text = textBox2.Text = textBox3.Text = textBox4.Text = textBox5.Text = textBox6.Text = "0";this.richTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标";}private void button1_Click(object sender, EventArgs e){//获取输入数据double bb1=Convert.ToDouble(this.textBox1.Text);double bb2=Convert.ToDouble(this.textBox2.Text);double bb3=Convert.ToDouble(this.textBox3.Text);double ll1=Convert.ToDouble(this.textBox4.Text);double ll2=Convert.ToDouble(this.textBox5.Text);double ll3=Convert.ToDouble(this.textBox6.Text);//检查输入格式的正确性if (bb1 >= 0 && bb1 <90 && bb2 >= 0 && bb2 < 60 && bb3 >= 0 && bb3 < 60){B = bb1 * 3600 + bb2 * 60 + bb3;}else{MessageBox.Show("纬度输入格式不正确！", "警告");return;}if (ll1 >= 0 && ll1<360 && ll2 >= 0 && ll2 < 60 && ll3 >= 0 && ll3 < 60){L = ll1 * 3600 + ll2 * 60 + ll3;}else{MessageBox.Show("经度输入格式不正确！", "警告");return;}double b = B/p;//获取用户选项值int tyd=radioButton1.Checked?3:6;int tq=radioButton3.Checked?1:2;//按带求带号以及中央子午线经度if (tyd == 6){N =(int)( L / (6*3600.0)) + 1;L0 = (6 * N - 3)*3600.0;l1 = (L-L0)/p;}if(tyd==3){N =(int)(L / (3*3600.0)+0.5);L0 = 3 * N * 3600.0;l1 =( L-L0)/p;}//采用克氏椭球参数的计算公式if (tq == 1){double c = Math.Pow( Math.Cos(b),2);double c1= Math.Sin(b) * Math.Cos(b);double c2=Math.Cos(b);double l2 = Math.Pow(l1,2);double n = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c;double a0 = 32140.404 - (135.3302 - (0.7092 - 0.0040 * c) * c) * c;double a4 = (0.25 + 0.00252 * c) * c - 0.04166;double a6 = (0.166 * c - 0.084) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.0083 - (0.1667 - (0.1968 + 0.0040 * c) * c) * c;x = 6367558.4969 * b - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) *c1; y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10){Y = N * i*10 + y1;}string tuoqiu = "采用克氏椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y + "\r\n" +"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n" + "Y=" + Y;}//采用1975国际椭球参数的计算公式if (tq == 2){double c = Math.Cos(b) * Math.Cos(b);double c1 = Math.Sin(b) * Math.Cos(b);double c2 = Math.Cos(b);double l2 = Math.Pow(l1, 2);double n = 6399596.652 - (21565.045 - (108.996 - 0.603 * c) * c) * c;double a0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * c) * c) * c;double a4 = (0.25 + 0.00253 * c) * c - 0.04167;double a6 = (0.167 * c - 0.083) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.00878 - (0.1702 - 0.20382 * c) * c;x = 6367452.1328 * b-(a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) * c1;y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10)//Y坐标加代号{Y = N * i * 10 + y1;}string tuoqiu = "采用1975国际椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y+"\r\n"+"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n"+"Y="+Y;}}private void button2_Click(object sender, EventArgs e){X = Convert.ToDouble(this.textBox7.Text);//获取输入数据Y = Convert.ToDouble(this.textBox8.Text);for (int i = 1; Y/i >=10; i = i * 10)//对Y坐标处理并求出中央子午线经度{y = Y - (int)(Y / i) * i-500000;L0 = 6 * (int)(Y / i) - 3;}//按6°带克氏椭球反算double bt = x / 6367558.4969*p;double BT = x / 6367558.4969;double c3=Math.Cos(BT)*Math.Cos(BT);double c4=Math.Sin(BT)*Math.Cos(BT);double Bf=(bt+(50221746+(293622+(2350+22*c3)*c3)*c3)*c4*Math.Pow(10,-10)*p)/p;double c5=Math.Pow(Math.Cos(Bf),2);double c6=Math.Sin(Bf)*Math.Cos(Bf);double Nf=6399698.902-(21562.267-(108.973-0.612*c5)*c5)*c5;double Z=y/(Nf*Math.Cos(Bf));double b2 = (0.5 + 0.003369 * c5) * c6;double b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5;double b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5;double b5=0.2-(0.1667-0.0088*c5)*c5;double z2=Math.Pow(Z,2);B = (Bf*p - (1 - (b4 - 0.12 *z2) * z2) * z2 * b2 * p)/3600.0;L = L0+((1 - (b3 - b5 * z2) * z2) * Z * p)/3600.0;this.richTextBox2.Text = "按6°带克氏椭球反算后，结果为："+ "\r\n"+"B="+ B + "\r\n" + "L=" + L;}}。

#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;}。

高斯投影正、反算及换带程序执行条件※数组投影选择T、换算点个数“Z=0 F≠0”、=0正算0、≠0反算※坐标系选择“54 ≠54”、=54换算为1954年北京坐标系输入54、≠54换算为1988年西安坐标系M、中央子午线经度（°′″）输入※大地坐标I、序列号B、L：大地纬度和经度（地理坐标）（°′″）※高斯平面坐标轴子午线I、序列号X、Y：高斯平面坐标（m） Z、轴子午线（°）输出※大地坐标子午收敛角N、序列号B、L：大地纬度和经度（地理坐标）（°′″） R、子午收敛角（°′″）※高斯平面坐标子午收敛角N、序列号X、Y：高斯平面坐标（m） R、子午收敛角（°′″）注：1、程序执行前必须进行数组定位。

如：Defm 10 T×2=5×2=102、Y坐标值要去掉带号及避免出现负值的500公里；4、本程序运算时，各已知数据、观测变量不会随之变化，可非常方便地进行各数据的核对;5、本程序在进行换带计算时采用的是间接换带计算法。

Prog GSXYDefm 10:TA“Z=0 F≠0”G“54 ≠54”Z:Fixm:I=0:「b」0:I=I+1◢J=2I-1:M=Z[J:L=Z[J+1:A=0=>Prog“3”:B=M:M=L+Z:Prog“3”:L=M:{BL}:M=B:Prog“2”: B=M:M=L:Prog“2”:L=M-Z:≠>X=M:Y=L:{XY}:B=X:L=Y⊿Z[J]=B:Z[J+1]=L:I<T=>Goto 0⊿G=54=>C=6399698.90178271:E=.006738525414684:≠>C=6399596.65198801:E=.006 739501819473⊿I=0:「b」0:I“N”=I+1◢J=2I-1:B=Z[J:L=Z[J+1:A≠0=>X=B:Y=L:Goto 2⊿S=sin B:G=54=>F=111134.8611B-(32 005.7799S+133.9238S∧3+.6973S∧5+.0039S∧7)cos B:≠>F=111133.0047B-(32009.857 S+133.9602S∧3+.6976S∧5+.0039S∧7)cos B⊿U=√Ecos B:V=√(1+U2:N=C÷V:W=tan B: M=cos B(Lπ÷180:X=F+NW(.5M2+1┛24(5-W2+9U2+4U∧4)M∧4+1┛720(61-58W2+W∧4)M∧6◢Y=N(M+1┛6(1-W 2+U 2)M ∧3+1┛120(5-18W 2+W ∧4+14U 2-58U 2W 2)M ∧5◢M=W ┛π(180M+60(1+3U 2+2U ∧4)M ∧3+12(2-W 2)M ∧5:Goto 3:「b 」2:W=E ﹣6X-3:G=54=>F=27.11115372595+9.024********W-.00579740442W 2-4.3532572E ﹣4W ∧3+4.857285E ﹣5W ∧4+2.15727E ﹣6W ∧5-1.9399E ﹣7W ∧6:≠>F=27.11162289465+9.024********W-.00579850656W2-4.3540029E ﹣4W ∧3+4.858357E ﹣5W ∧4+2.15769E ﹣6W ∧5-1.9404E ﹣7W ∧6⊿U=√Ecos F:V=√(1+U 2:Q=YV ÷C:W=tan F:M=F-(1+U 2)W ┛π(90Q 2-7.5(5+3W 2+U 2-9U 2W 2)Q ∧4+.25(61+90W 2+45W ∧4)Q ∧6:Prog “3”:B=M ◢M=Z+1┛(πcos F)(180Q-30(1+2W 2+U 2)Q ∧3+1.5(5+28W 2+24W ∧4)Q ∧5:Prog “3”:L=M ◢M=W ┛π(180Q-60(1+W 2-U 2)Q ∧3+12(2+5W 2+3W ∧4)Q ∧5:「b 」3:Prog “3”:R=M ◢ I<T=>Goto 1⊿“END ”概要说明：我国的经度范围西边自73°起，东边至135°，可分成6°带共11带或3°共22带。

高斯投影正反算程序设计一．程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。

（2），程序主要的算法来自于教材。

但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球（3{{//{p ublicdoubledLongitude;publicdoubledLatitude;publicdoubledHeight;}//笛卡尔坐标//CartesianCoordinatepublicstructCRDCARTESIAN{publicdoublex;publicdoubley;publicdoublez;}publicForm1(){InitializeComponent();}privatevoidbutton1_Click(objectsender,EventArgse) {try{tt=}{}{ee=}{ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0){a=6378137.0;}constdoublepai=3.1415926;doubleb=Math.Sqrt(a*a*(1-ee*ee));doublec=a*a/b;doubleepp=Math.Sqrt((a*a-b*b)/b/b);CRDGEODETICpcrdGeo;CRDCARTESIANpcrdCar;doublemidlong;//{}ai;//{}pai;//intnum;//带号midlong=0;//默认值，需要制定分带try{tt=}catch{MessageBox.Show("Choose3/6error!");return;if(pareTo("3度带")==0){num=Convert.ToInt32(deglon/6+1);midlong=(6*num-3)/180.0*pai;}if(pareTo("6度带")==0){num=Convert.ToInt32((deglon+1.5)/3);}a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;a6=m6/32.0+m8/16.0;a8=m8/128.0;doubleB=pcrdGeo.dLatitude;doublesb=Math.Sin(B);doublecb=Math.Cos(B);doubles2b=sb*cb*2;doubles4b=s2b*(1-2*sb*sb)*2;doubles6b=s2b*Math.Sqrt(1-s4b*s4b)+s4b*Math.Sqrt(1-s2b*s2b); X=a0*B-a2/2.0*s2b+a4*s4b/4.0-a6/6.0*s6b;pcrdCar.x=Nscnb*lp*lp/2.0+Nscnb*cosb*cosb*Math.Pow(lp,4)*(5-t*t+9*ita*ita+4*Math.Pow(ita,4))/24.0 +Nscnb*Math.Pow(cosb,4)*Math.Pow(lp,6)*(61-58*t*t+Math.Pow(t,4))/720.0+X;pcrdCar.y=Ncosb*Math.Pow(lp,1)+Ncosb*cosb*cosb*(1-t*t+ita*ita)/6.0*Math.Pow(lp,3)+Ncosb*Math.Pow(l p,5)*Math.Pow(cosb,4)*(5-18*t*t+Math.Pow(t,4)+14*ita*ita-58*ita*ita*t*t)/120.0;if(pcrdCar.y<0)pcrdCar.y+=500000;richTextBox1.Text="Results:\nX:\t"+Convert.ToString(pcrdCar.x)+"\nY:\t"+Convert.ToString(pcrdCar.y );}{try{tt=}{}{ee=}if(pareTo("WGS-84")==0){a=6378137.00;ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0) {a=6378137.0;ee}constdoublepai=doubleb=Math.Sqrt(a*a*(1-ee*ee));//求Xtry{tt=}{}if(pareTo("3度带")==0){midlong=num*3*pai/180;}if(pareTo("6度带")==0){midlong=(6*num-3)*pai/180;}b=Math.Sqrt(a*a*(1-ee*ee));c=a*a/b;epp=Math.Sqrt(a*a-b*b)/b;doublem0,m2,m4,m6,m8;doublea0,a2,a4,a6,a8;m0=a*(1-ee*ee);m2=3.0/2.0*m0*ee*ee;m4=5.0/4.0*ee*ee*m2;m6=7.0/6.0*ee*ee*m4;m8=9.0/8.0*ee*ee*m6;{}tf=Math.Tan(Bf);Vf=Math.Sqrt(1+epp*epp*Math.Cos(Bf)*Math.Cos(Bf));Nf=c/Vf;doubleynf=pcrdCar.y/Nf;pcrdGeo.dLatitude=Bf-1.0/2.0*Vf*Vf*tf*(ynf*ynf-1.0/12.0*Math.Pow(ynf,4)*(5+3*tf*tf+itaf*itaf-9*Math.Pow(itaf*tf,2))+1.0/360.0*(61+90*tf*tf+45*Math.Pow(tf,4))*Math.Pow(ynf,6));pcrdGeo.dLongitude=(ynf/Math.Cos(Bf)-(1+2*tf*tf+itaf*itaf)*Math.Pow(ynf,3)/6.0/Math.Cos(Bf)+ (5+28*tf*tf+24*Math.Pow(tf,4)+6*itaf*itaf+8*Math.Pow(itaf*tf,2))*Math.Pow(ynf,5)/120.0/Math.Cos(Bf ));pcrdGeo.dLongitude=pcrdGeo.dLongitude+midlong;//pcrdGeo.dLatitude=pcrdGeo.dLatitude;richTextBox2.Text="Results:\nLatitude:"+Convert.ToString(pcrdGeo.dLatitude)+"\nLongtitude:"+Conver t.ToString(pcrdGeo.dLongitude);}privatevoidlabel13_Click(objectsender,EventArgse){}}}米之间，。

高斯投影正反算编程一．高斯投影正反算基本公式（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].L0)!=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)/720;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].L0)!=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*sin(8*B1 )/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2; }。

§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。

8.3.1高斯投影坐标正算公式： B,l ⇒ x,y高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。

由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为l 的偶函数，y 为l 的奇函数；0330'≤l ，即20/1/≈''''ρl ，如展开为l 的级数，收敛。

ΛΛ+++=++++=553316644220l m l m l m y l m l m l m m x （8-33）式中Λ,,10m m 是待定系数，它们都是纬度B 的函数。

由第三个条件知：qyl x l y q x ∂∂-=∂∂∂∂=∂∂, (8-33)式分别对l 和q 求偏导数并代入上式ΛΛΛΛ----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m (8-34) 上两式两边相等，其必要充分条件是同次幂l 前的系数应相等，即ΛΛΛΛΛΛdq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)(8-35)是一种递推公式，只要确定了0m 就可依次确定其余各系数。

由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当0=l时有：0m X x == (8-36) 顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2==== 得：B V cB N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-37,38)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-= (8-39)依次求得6543,,,m m m m 并代入(8-33)式，得到高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ (8-42) 8.3.2高斯投影坐标反算公式x,y ⇒B,l投影方程：),(),(21y x l y x B ϕϕ== (8-43)满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。

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

在同一大地坐标系中(例如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)-(35 *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表示。

高斯投影坐标正反算编程报告The Standardization Office was revised on the afternoon of December 13, 2020高斯投影坐标正反算编程报告1. 编程思想进行高斯投影坐标正反算的编程需要牵涉到大量的公式，为了使程序条理更清楚，各块的数据复用性更强，这里采取了结构化的编程思想。

程序由四大块组成。

文件用于存放main()函数，是整个程序的入口。

通过结构化的编程尽力使main()函数变得简单。

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

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

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

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

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

2. 计算模型高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos lt t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ高斯投影反算公式()()()()22242552233642542222328624285cos 12021cos 6cos 459061720935242f f f f f ff ff ff f f ff ff ff f f f ff f ff f f t t t B N y t B N y B N y l y t t y NM t yt tNM t y N M t B B ηηηηη+++++++-=++--+++-=3. 程序框图4. 计算结果5.附录：程序代码517"<<endl;cin>>myX>>myY;Fansuan myFansuan1(myX,myY);();}////////////////////////////////////////////////////////////////////////#include ""//////////////////////////////////////////////////////////////////////// Construction/Destruction//////////////////////////////////////////////////////////////////////Fansuan::Fansuan(){}Fansuan::Fansuan(double X,double Y){this->X=X;this->Y=Y;//初始化x,yN=(int)(Y/1000000);//取出带号L0=6*N-3;//初始化该带号的中央经线经度x=X;y=Y-1000000*N-500000;beta=x/;//初始化beta，弧度单位betaSecond=beta*rouSecond;//初始化beta，秒单位betaDegree=betaSecond/3600;//初始化beta,整度数单位Bf=beta+(+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(b eta)*cos(beta))*(1e-10)*sin(beta)*cos(beta);//初始化Bf,弧度单位BfSecond=Bf*rouSecond;//初始化Bf，秒单位BfDegree=BfSecond/3600;//初始化Bf，整度数单位Nf= Z=y/(Nf*cos(Bf));b2=+*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b4=++*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b3= b5= Bsecond=BfSecond-(1-*Z*Z)*Z*Z)*Z*Z*b2*rouSecond;//计算大地经度B，单位为秒B=Bsecond/3600;//用整度数表示Blsecond=(1-(b3-b5*Z*Z)*Z*Z)*Z*rouSecond;//计算经度差l，单位为秒l=lsecond/3600;//用整度数表示lL=L0+l;//计算大地经度L}double Fansuan::getB(){return B;}double Fansuan::getL(){return L;}void Fansuan::printLocation(){printf("反算得大地坐标为：大地纬度B= %f°大地经度L=%f°\n",B,L);}Fansuan::~Fansuan(){}。