高斯投影正反算编程

高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。（2）相关流程图1）正算输入大地坐标B，L和经差L0选择带宽3/6度带计算带号计算弧长计算平面坐标x,y打印x,y计算带号计算弧长计算平面坐标x,y打印x,y开始选择椭球参数3度带6度带2）反算开始输入自然值坐标x,y和经差L0利用迭代算法求解底点纬度利用公式计算B和L打印B和L选择椭球参数三．编程的相关代码（1）正算#includestdio.h#includestdlib.h#includemath.h#includeassert.h#definepi(4*atan(1.0))inti;structjin{doubleB;doubleL;doubleL0;};structjing[100];main(intargc,double*argv[]){FILE*r=fopen(a.txt,r);assert(r!=NULL);FILE*w=fopen(b.txt,w);assert(r!=NULL);inti=0;while(fscanf(r,%lf%lf%lf,&g[i].B,&g[i].L,&g[i].L0)!=EOF){doublea,b;intzuobiao;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;}//选择坐标系//doublef=(a-b)/a;doublee,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一，第二曲率//doublem0,m2,m4,m6,m8;doublea0,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;doubleBmiao,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;doubledb;db=pi/180.0/3600.0;doubleB1,L1,l;B1=Bmiao*db;L1=Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//doubleT=tan(B1)*tan(B1);doublen=e2*e2*cos(B1)*cos(B1);doubleA=l*cos(B1);doubleX,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;//求弧长//doubleN=a/sqrt(1-e*e*sin(B1)*sin(B1));intZonewide;intZonenumber;printf(\n请输入带宽：3度带或6度带Zonewide=);scanf(%d,&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}elseif(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=%lfy=%lf\n,x,y);fprintf(w,%lf%lf\n,x,y);//输出结果到文本文件//}fclose(r);fclose(w);system(pause);return0;}（2）反算#includestdio.h#includestdlib.h#includemath.h#includeassert.h#definepi(4*atan(1.0))doubleX,Y,B1,B2,B3,F,t;doublem0,m2,m4,m6,m8;doublea0,a2,a4,a6,a8,a1,b1;doubleBB,LL,Bf;doublee,e1;intd,m,s,i,zuobiao;doublesort(double,double);structjin{doublex;doubley;doubleL0;};structjing[100];//x,y,L0为输入量：x,y坐标和中央子午线经度//main(intargc,double*argv[]){FILE*r=fopen(c.txt,r);assert(r!=NULL);FILE*w=fopen(d.txt,w);assert(r!=NULL);inti=0;while(fscanf(r,%lf%lf%lf,&g[i].x,&g[i].y,&g[i].L0)!=EOF)//文件为空，无法打开//{doublea1=6378245.0000000000;//克拉索夫斯基椭球参数//doubleb1=6356863.0187730473;doublea75=6378140.0000000000;//1975国际椭球参数//doubleb75=6356755.2881575287;doublea84=6378137.0000000000;//WGS-84系椭球参数//doubleb84=6356752.3142000000;doubleM,N;//mouyou圈曲率半径，子午圈曲率半径//doublet,n;doubleA,B,C;doubleBB,LL,Bf,LL0,BB0;doublea,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);//调用求解底点纬度的函数//doubleq=sqrt(1-e*e*sin(Bf)*sin(Bf));doubleG=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;doubleH,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);//利用公式求解经纬度//intBdu,Bfen,Ldu,Lfen;doubleBmiao,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);return0;}doublesort(doublea,doubleb){doublee,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);doublem0,m2,m4,m6,m8;doublea0,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);//利用迭代算法求解底点纬度//returnB2;}