高斯投影正反算编程一高斯投影正反算基本公式文档格式.docx

1、3度带选择椭球参数2）反算输入自然值坐标x,y利用迭代算法求解底点纬度利用公式计算B和L打印B和L三编程的相关代码（1）正算# include stdio.hstdlib.hmath.hassert.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=fopen（a.txt,r）; assert（r!=NULL）; FILE *w=fopen（b.txtw int i=0; while（fscanf（

2、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（%dzuobiao）; 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; /选择坐标系/ double f=（a-b）/a; double

3、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*m8/16; a4=m4/8+3*m6/16+7*m8/32; a6=m6/32+m8/16; a8=m8/128; d

4、ouble 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）（gi.L0-（int）（gi.L0）*100.0）*60.0+（gi.L0*100-（int）（gi.L0*

5、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+a8*sin（8*B1）/8;/求弧长/ double N=a/sqrt（1-e*e

6、*sin（B1）*sin（B1）; int Zonewide; int Zonenumber;n请输入带宽：3度带或6度带Zonewide=Zonewide）; if（Zonewide=3） Zonenumber=（int）（gi.L-Zonewide/2）/Zonewide+1）; else if（Zonewide=6） Zonenumber=（int）gi.L/Zonewide+1; else错误 exit（0）; /选择带宽/ double FE=Zonenumber*1000000+500000;/改写为国家通用坐标/ y=FE+N*A+A*A*A*N*（1-T*T+n*n）/6+A*

7、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;n所选坐标系的转换结果：x=%lf y=%lfn,x,y）; fprintf（w,%lf %lfn/输出结果到文本文件/ fclose（r）; fclose（w）; system（pause return 0;（2）反算double X,Y,B1,B2,B3,F,t;dou

8、ble 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,double）; double x; double y; double L0;/x,y,L0为输入量：x,y坐标和中央子午线经度/c.txtd.txtgi.x,&gi.y,&=EOF）/文件为空，无法打开/ double a1=6378245.0000000000;/克拉索夫斯基椭球参数/ double b1=6356863.0187730473; double a7

9、5=6378140.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;n选择参考椭球：1=克拉索夫斯基椭球，2=1975国际椭球，3=WGS-84系椭球: a=a1; b=b1; a=a75; b=b75; a=a84; b

10、=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*M）+gi.y*t*H/（24*M）*（5+3*B+C-9*B*C）-gi.y*t*I/（720*M）*（61

11、+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（BB0/pi*180）; Bfen=int（BB0/pi*180）*60-Bdu*60）; B

12、miao=BB0/pi*180*3600-Bdu*3600-Bfen*60; /将弧度转化为角度/%d度%d分%lf秒 %d度%d分%lf秒 n,Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao）; fprintf（w,%d%d%lf”%d%d%lf”n/将结果输出到文本文件/double sort（double a,double b） e=sqrt（1-（b/a）*（b/a）; e1=sqrt（a/b）*（a/b）-1）; 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;