高斯投影正反算编程

高斯投影正反算编程一．高斯投影正反算基本公式（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）正算

# 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))*1 00.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))*10 0.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);