高斯投影正反算編程

1、-作者xxxx-日期xxxx高斯投影正反算編程【精品文檔】高斯投影正反算編程一高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式 以上主要通過大地測(cè)量學(xué)基礎(chǔ)課程得到，這不進(jìn)行詳細(xì)的推導(dǎo)，只是列出基本公式指導(dǎo)編程的進(jìn)行。二編程的基本方法和流程圖（1）編程的基本方法高斯投影正反算基本上運(yùn)用了所有的編程基本語句，本文中是利用C+語言進(jìn)行基本的設(shè)計(jì)。高斯正算中對(duì)橢球參數(shù)和帶寬的選擇主要運(yùn)用了選擇語句。而高斯反算中除了選擇語句的應(yīng)用，在利用迭代算法求底點(diǎn)緯度還應(yīng)用了循環(huán)語句。編程中還應(yīng)特別注意相關(guān)的度分秒和弧度之間的相互轉(zhuǎn)換，這是極其重要的。（2）相關(guān)流程圖1）正算輸入大地坐標(biāo)B，L和經(jīng)

2、差L0選擇帶寬3/6度帶計(jì)算帶號(hào)計(jì)算弧長(zhǎng)計(jì)算平面坐標(biāo)x,y打印x,y計(jì)算帶號(hào)計(jì)算弧長(zhǎng)計(jì)算平面坐標(biāo)x,y打印x,y開始 6度帶3度帶選擇橢球參數(shù) 2）反算選擇橢球參數(shù)開始輸入自然值坐標(biāo)x,y和經(jīng)差L0利用迭代算法求解底點(diǎn)緯度利用公式計(jì)算B和L打印B和L 三編程的相關(guān)代碼（1）正算# include ""# include ""# include ""# include ""#define pi (4*atan(1.0)int i;struct jin double B; double L;double L0;stru

3、ct jin g100; 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",&gi.B,&gi.L,&gi.L0)!=EOF) double a,b; int zuobiao; printf("n請(qǐng)輸入坐標(biāo)系：

4、北京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; /選擇坐標(biāo)系/ double f=(a-b)/a; double e,e2; e=sqrt(2*f-f*f); e2=sqrt(a/b)*(a/b)-1);/求橢球的第一，第

5、二曲率/ 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)(gi.B)*3600.0+(i

6、nt)(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)(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,L

7、1,l; B1=Bmiao*db; L1= Lmiao*db; l=L1-L0miao*db;/角度轉(zhuǎn)化為弧度/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;/求弧長(zhǎng)/ double N=a/sqrt(1-e*e*sin(B1)*sin(B1); int Zonewide; int Zonenumber; printf(

8、"n請(qǐng)輸入帶寬：3度帶或6度帶Zonewide="); scanf("%d",&Zonewide); getchar(); if(Zonewide=3) Zonenumber=(int)(gi.L-Zonewide/2)/Zonewide+1); else if(Zonewide=6) Zonenumber=(int)gi.L/Zonewide+1; else printf("錯(cuò)誤"); exit(0); /選擇帶寬/ double FE=Zonenumber*1000000+500000;/改寫為國(guó)家通用坐標(biāo)/ y=FE+N

9、*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所選坐標(biāo)系的轉(zhuǎn)換結(jié)果：x=%lf y=%lfn",x,y);fprintf(w,"%lf %lfn",x,y);/輸出結(jié)果到文本文件/ fclose(r); f

10、close(w);system("pause");return 0; （2）反算# include ""# include ""# include ""# include ""#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(doub

11、le,double);struct jin double x; double y; double L0;struct jin g100;/x,y,L0為輸入量：x,y坐標(biāo)和中央子午線經(jīng)度/ 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",&a

12、mp;gi.x,&gi.y,&gi.L0)!=EOF)/文件為空，無法打開/ double a1=6378245.0000000000;/克拉索夫斯基橢球參數(shù)/ double b1=6356863.0187730473; double a75=6378140.0000000000;/1975國(guó)際橢球參數(shù)/ double b75=6356755.2881575287; double a84=6378137.0000000000;/WGS-84系橢球參數(shù)/ double b84=6356752.3142000000; double M,N;/mouyou圈曲率半徑，子午圈曲率半徑/

13、 double t,n; double A,B,C; double BB,LL,Bf,LL0,BB0; double a,b; printf("n選擇參考橢球：1=克拉索夫斯基橢球，2=1975國(guó)際橢球，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); /調(diào)用

14、求解底點(diǎn)緯度的函數(shù)/ 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*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.

15、0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);/利用公式求解經(jīng)緯度/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; /將弧度轉(zhuǎn)化為角度/ printf(

16、"n所選坐標(biāo)系的轉(zhuǎn)換結(jié)果：%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);/將結(jié)果輸出到文本文件/ 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+