《高斯投影正反算编程》由会员分享，可在线阅读，更多相关《高斯投影正反算编程（6页珍藏版）》请在金锄头文库上搜索。
1、高斯投影正反算编程一、 题目： 已知部分数据，根据高斯投影正反算思想进行编程，并采用克氏椭球，按3或 6带投影。 正算：已知大地坐标B 、L，二、 已知数据：正算： B=51.38439023 L=111.02131360反算： x=5724004.723 y=19502559.920三、 计算结果：正算结果： x=5724004.723 y=19502559.920反算结果： B=51.38439023 L=111.02131360四、 计算流程图：输入国家统一坐标x、y输入大地坐标B、Li=1,2开始i=1(正算) i=2（反算）选择椭球参数选择椭球参数（不予画出） 计算带号，将国家统一坐
2、标化为高斯平面坐标。 （不予画出）计算带号记算x和y计算B和L打印x 和 y打印B和L五、源代码：#include gaosi.h#include math.h#include stdio.h#include tchar.h#include stdlib.h#define pi 3.141592653589793#define rho 206265void Calculateellipse2plane(double B,double L);void Calculateplane2ellipse(double x,double y);double Dms2Rad(double Dms);doub
3、le D2Dms(double D);double Dms2D(double Dms);int main()double B=0,L=0;double x=0,y=0;int i=0;printf(如使用高斯投影坐标正算，请输入1；反算，请输入2n);scanf_s(%d,&i);if(i=1)printf(请输入大地坐标纬度B(度分秒):n);scanf_s(%lf,&B);printf(请输入大地坐标经度L(度分秒):n);scanf_s(%lf,&L);Calculateellipse2plane(B,L);if(i=2)printf(请输入国家统一坐标x(m):n);scanf_s(%
4、lf,&x);printf(请输入国家统一坐标y(m):n);scanf_s(%lf,&y);Calculateplane2ellipse(x,y);return 0;void Calculateellipse2plane(double B,double L) /高斯投影正算主体double l=0,Lo=0,a0=0,a3=0,a4=0,a5=0,a6=0,n=0,c=0;double x=0,y=0;double m=0,p=0,q=0;int N=0,i=0; /带号printf(如使用6带请输入1，使用3带请输入2n);scanf_s(%d,&i);if(i=1) /已知a点在6带的带
5、号和中央子午线经度N=int(Dms2D(L)/6);Lo=6*N-3;if(i=2) /已知a点在3带的带号和中央子午线经度N=int(Dms2D(L)+1.5)/3);Lo=3*N;l=(Dms2D(L)-Dms2D(Lo)*3600/rho; /单位为 B=Dms2Rad(B);L=Dms2Rad(L);Lo=Dms2Rad(Lo);l=pow(l,2);c=pow(cos(B),2); /c=cos(B)的平方n=6399698.902-(21562.267-(108.973-0.612*c)*c)*c;a0=32140.404-(135.3302-(0.7092-0.0040*c)*
6、c)*c;a3=(0.3333333+0.001123*c)*c-0.1666667;a4=(0.25+0.00252*c)*c-0.04166;a5=0.0083-(0.1667-(0.1968+0.004*c)*c)*c;a6=(0.166*c-0.084)*c;/计算平面坐标，并化为国家统一坐标m=sin(B)*cos(B);p=1+(a3+a5*l)*l;x=6367558.4969*B-(a0-(0.5+(a4+a6*l)*l)*l*n)*m;y=(p)*sqrt(l)*n*cos(B);y=y+500000+N*1000000;printf(国家统一坐标为:n);printf(x=
7、%.4lfn,x);printf(y=%.4lfn,y);system(pause);void Calculateplane2ellipse(double x,double y) /高斯投影反算主体double l=0,Lo=0,b2=0,b3=0,b4=0,b5=0,nf=0,Z=0,beta=0,Bf=0;double B=0,L=0;double c=0,d=0,m=0;int N=0,i=0; /带号N=int(y/1000000);printf(如使用6带请输入1，使用3带请输入2n);scanf_s(%d,&i);if (i=1) /已知a点在6带的带号和中央子午线经度Lo=6*N
8、-3;if (i=2) /已知a点在3带的带号和中央子午线经度Lo=3*N;y=y-N*1000000-500000;beta=x/6367558.4969; /单位为弧度c=cos(beta);c=pow(cos(beta),2); /c=cos(beta)的平方Bf=beta+(50221746+(293622+(2350+22*c)*c)*c)*sin(beta)*cos(beta)*(1e-10);d=pow(cos(Bf),2); /d=cos(Bf)的平方nf=6399698.902-(21562.267-(108.973-0.612*d)*d)*d;Z=y/(nf*cos(Bf)
9、;Z=pow(Z,2);b2=(0.5+0.003369*d)*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*d)*d;b4=0.25+(0.16161+0.00562*d)*d;b5=0.2-(0.1667-0.0088*d)*d;/计算大地坐标B,Lm=(1-(b4-0.12*Z)*Z)*Z*b2;l=(1-(b3-b5*Z)*Z)*sqrt(Z)*rho;B=Bf*rho-m*rho;l=D2Dms(l/3600);B=D2Dms(B/3600);L=Lo+l;printf(大地坐标为:n);printf(B=%.6lfn,B);print
10、f(L=%.6lfn,L);system(pause);double D2Dms(double D) /度化为度分秒double Degree, Miniute;double Second;int Sign;double Dms,Rad;if(D = 0)Sign = 1;elseSign = -1;Rad=D*pi/180;Rad = fabs(Rad * 180.0 / pi);Degree = floor(Rad);Miniute = floor(fmod(Rad * 60.0, 60.0);Second = fmod(Rad * 3600.0, 60.0);Dms = Sign * (
11、Degree + Miniute / 100.0 + Second / 10000.0);return Dms;double Dms2Rad(double Dms) /度分秒化为弧度double Degree, Miniute;double Second;int Sign;double Rad;if(Dms = 0)Sign = 1;elseSign = -1;Dms = fabs(Dms);Degree = floor(Dms);Miniute = floor(fmod(Dms * 100.0, 100.0);Second = fmod(Dms * 10000.0, 100.0);Rad =
12、 Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * pi / 180.0;return Rad;double Dms2D(double Dms) /度分秒化为度double Degree, Miniute;double Second;int Sign;double D;if(Dms = 0)Sign = 1;elseSign = -1;Dms = fabs(Dms);Degree = floor(Dms);Miniute = floor(fmod(Dms * 100.0, 100.0);Second = fmod(Dms * 10000.0, 100.0);D=Sign * (Degree + Miniute / 60.0 + Second / 3600.0);return D;
