1、高斯投影正反算编程一高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。二编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用 C+语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。（2）相关流程图1）正算输入大地坐标B，L和经差 L0选择带宽3/6 度带计算带号计算弧长计算平面坐标
2、x,y打印 x,y计算带号计算弧长计算平面坐标 x,y打印 x,y开始选择椭球参数3 度带 6 度带2）反算开始输入自然值坐标 x,y和经差 L0利用迭代算法求解底点纬度利用公式计算 B和 L打印 B 和 L选择椭球参数三编程的相关代码（1）正算# include “stdio.h“# include “stdlib.h“# include “math.h“# include “assert.h“#define pi (4*atan(1.0)int i;struct jindouble B;double L;double L0;struct jin g100;main(int argc, do
3、uble *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“, int zuobiao;printf(“n 请输入坐标系：北京 54=1，西安80=2，WGS84=3 ：“);scanf(“%d“,getchar();if(zuobiao=1)a=6378245;b=6356863.0187730473; if(zuobiao=2)a=6378140;b=6356755.2881575287
4、;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*
5、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)(gi.B)*3600.0+(int)(gi.B-(int)(gi.B)*100.0)*60.0+(gi.B*100-(int)(gi.B*100)*100.0;Lmiao=(int)(gi.L)*3600.0+(int)(gi.L-(int)(gi.L)*100.0)*60.0+(gi.L*100-(int)(gi.L*100)*100.0;L0miao=(int)(gi.L0)*3600.0+(int)(g
6、i.L0-(int)(gi.L0)*100.0)*60.0+(gi.L0*100-(int)(gi.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+
7、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“,getchar();if(Zonewide=3)Zonenumber=(int)(gi.L-Zonewide/2)/Zonewide+1);else if(Zonewide=6)Zonenumber=(int)gi.L/Zonewide+1;elseprintf(“错误“);exit(0);/选择带宽/doubleFE=Zon
8、enumber*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=%lfn“,x,y);fprintf(w,“%lf %lfn“,x,y);/输出结
9、果到文本文件/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,do
10、uble);struct jindouble x;double y;double L0;struct jin g100;/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“,/克拉索夫斯基椭球参数/double b1=6356863.0187730473;double a75=6378140.
11、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“,getchar();i
12、f(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=gi.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-gi.y*t*A/(2*
13、M)+gi.y*t*H/(24*M)*(5+3*B+C-9*B*C)-gi.y*t*I/(720*M)*(61+90*B+45*B*B);LL0=gi.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
14、(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;doubl
15、e 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=gi.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=(gi.x-F)/a0;B3=B1;B1=B2; while(fabs(B3-B2)10e-10);/利用迭代算法求解底点纬度/return B2;
