高斯投影正反算编程.doc

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

（2）高斯反算基本公式

以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。
二．编程的基本方法和流程图
（1）编程的基本方法
高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。
（2）相关流程图
1）正算

输入大地坐标B，L
和经差L0
选择带宽3/6度带
计算带号
计算弧长

计算平面坐标x,y

打印x,y
计算带号

计算弧长

计算平面坐标x,y

打印x,y

开始

6度带
3度带
选择椭球参数

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))*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;

}