高斯投影正反算编程(可编辑修改word版)

高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式

（2）高斯反算基本公式

以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。

二．编程的基本方法和流程图

（1）编程的基本方法

高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。

（2）相关流程图

1）正算

选择带宽 3/6 度带

计算带号

输入大地坐标 B ，L

和经差 L0

6 度带

3 度带

选择椭球参数

计算带号

计算弧长

计算平面坐标 x,y

打印 x,y

开始

计算平面坐标 x,y

计算弧长

打印 x,y

开始

输入自然值坐标x,y

和经差L0

选择椭球参数

利用迭代算法

求解底点纬度

利用公式计算B

和L

打印B 和L

2）反算

三．编程的相关代码（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 g[100];

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",&g[i].B,&g[i].L,&g[i].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;

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/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)(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;

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*sin(B1)*sin(B1));

int Zonewide;

int Zonenumber;

printf("\n 请输入带宽：3 度带或6 度带Zonewide=");

scanf("%d",&Zonewide);

getchar();

if(Zonewide==3)

{

Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);

}

else if(Zonewide==6)

{

Zonenumber=(int)g[i].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)/720;

printf("\n 所选坐标系的转换结果：x=%lf y=%lf\n",x,y);

fprintf(w,"%lf %lf\n",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 BB,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 g[100];//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",&g[i].x,&g[i].y,&g[i].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;

double 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(Bf));

double G=cos(Bf);

M=a*(1-e*e)/(q*q*q);

N=a/q;

double H,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); //利用公式求解经纬度//

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,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");

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;

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);//利用迭代算法求解底点纬度// return B2;

}