高斯投影正反算编程

1、高斯投影正反算编程一高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式 以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。二编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C+语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。（2）相关流程图1）正算输入大地坐标B，L和经差L0选择带宽3/6度带计算带号计算弧长计算平面坐标x,y

2、打印x,y计算带号计算弧长计算平面坐标x,y打印x,y开始 6度带3度带选择椭球参数 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 jin double B; double L;double L0;struct jin g100; main(int argc, double *argv)FILE *r

3、=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,&gi.B,&gi.L,&gi.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

4、; 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/

5、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)(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*1

6、00)*100.0; L0miao=(int)(gi.L0)*3600.0+(int)(gi.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;

7、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)(gi.L-Zonewide/2)/Zonewide+1); else if(Zonewide=6) Zonenumb

8、er=(int)gi.L/Zonewide+1; else printf(错误); exit(0); /选择带宽/ double FE=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)/72

9、0; printf(n所选坐标系的转换结果：x=%lf y=%lfn,x,y);fprintf(w,%lf %lfn,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 B

10、B,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 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,&gi.x

11、,&gi.y,&gi.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; do

12、uble 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(

13、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*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

14、*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,Bfe

15、n,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;